TEXT(2D ATMOSPHERIC BOUNDARY LAYER: T306
TITLE
  DISPLAY
  The case considered is the simulation of a turbulent neutral
  atmospheric boundary layer using the standard k-e model
  and fully-rough logarithmic wall functions. The height of the 
  solution domain is below 150m, and so the flow can be 
  considered as a constant stress layer. In this region the 
  velocity and turbulence variables are given by:
 
      u/u*=ln(z/z0)/K    k=u*^2/cmucd**0.5   enut=K*z*u*
      e=u*^3/(K*z)   or  w=sqrt(k)/([K*z*cmucd]**0.25)
 
  where u*=friction velocity, K=0.41, cmucd=0.09 and z0 is the
  roughness height. A calculation is made using a reference wind
  speed of 4m/s at a reference height of Zr=100m. The wind enters
  from the west boundary over terrain with an aerodynamic
  roughness height z0=0.3m. The domain height is set equal to
  the reference height, and the flow exits the solution domain
  at the east boundary which is located 2.*Zr downstream of the
  inlet. Under these conditions, the friction velocity   
  u*=0.282m/s. The inlet profiles correspond to a fully-developed
  neutral atmospheric boundary. The calculation demonstrates that
  PHOENICS preserves these inlet profiles and returns a uniform
  pressure field. The case is also simulated by using the
  Realisable k-e model and the Wilcox-Kolmogorov k-w model.
  ENDDIS
  PHOTON USE
   p
 
 
 
   up z
   gr y 1
   con p1 y 1 fi;.1;pa
   con u1 y 1 fi;.1;pa
   con ke y 1 fi;.1;pa
   con enut y 1 fi;.1;pa
   con ep y 1 fi;.1;pa
   con udis y 1 fi;.1;pa
   con kdis y 1 fi;.1;pa
  ENDUSE
CHAR(CTURB)
BOOLEAN(KWMOD);KWMOD=F
REAL(TETA,VKC,XLEN,DTF,OMEGIN,OMEGSKY,zdz0)
REAL(USTAR,USTAR2,USTARX,USTARY,Z0,DOMH,UIN,UINX,UINY)
REAL(KCNS,USTAR3,EPCON)
TETA=0.*3.14159265/3.  ! wind angle
UIN =4.0               ! reference velocity
DOMH=100.              ! reference height = domain height 
Z0=0.3                 !  roughness height
VKC=0.41               ! von karman's constant
UINX=UIN*COS(TETA);UINY=UIN*SIN(TETA)
zdz0=domh/z0
USTAR=UIN*VKC/LOG(zdz0)  ! friction velocity
USTAR2=USTAR*USTAR;KCNS=USTAR2/0.3;USTAR3=USTAR2*USTAR
USTARX=USTAR*COS(TETA);USTARY=USTAR*SIN(TETA)
EPCON=USTAR3/VKC
    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
XLEN=2.*DOMH
GRDPWR(X,100,XLEN,1.)
    GROUP 4. Y-direction grid specification
YVLAST=100.0;NY=1
    GROUP 5. Z-direction grid specification
NREGZ=2
IREGZ=1;GRDPWR(Z,1,4.*Z0,1.0)
IREGZ=2;GRDPWR(Z,49,-(DOMH-4.*Z0),1.06)
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,W1,P1);SOLUTN(P1,Y,Y,Y,P,P,P)
SOLUTN(U1,P,P,P,P,P,N);SOLUTN(W1,P,P,P,P,P,N)
SOLVE(C1);SOLUTN(C1,Y,Y,Y,P,P,P)
MESG( Enter the required turbulence model:
MESG(  KE   -  Standard k-e model (default)
MESG(  KO   -  Wilcox   k-o model
MESG(  RKE  -  Realisable k-e model
MESG(
READVDU(CTURB,CHAR,KE)
CASE :CTURB: OF
WHEN KE,2
+ TEXT(2D KE ATMOSPHERIC BOUNDARY LAYER: T306
+ MESG(Standard k-e model
+ TURMOD(KEMODL)
WHEN KO,2
+ TEXT(2D K-F ATMOSPHERIC BOUNDARY LAYER: T306
+ MESG(k-omega model
+ TURMOD(KOMODL);STORE(EP);KWMOD=T
WHEN RKE,3
+ TEXT(2D RK-E ATMOSPHERIC BOUNDARY LAYER: T306
+ MESG(Realisable k-e model
+ TURMOD(KEREAL);STORE(C1E)
ENDCASE
STORE(ENUT,STRS,LEN1,GEN1,UIN,KEIN,EPIN,UDIS,KDIS)
  *** store fully-developed boundary-layer values)
(stored of UIN is (:ustar:*loge(ZG/:Z0:)/:vkc:)!zslstr)
(stored of KEIN is :kcns:!zslstr)
(stored of EPIN is (:epcon:/ZG)!zslstr)
IF(KWMOD) THEN
+ STORE(OMIN)
+(stored of OMIN is :epcon:/(ZG*0.09*:kcns:)!zslstr)
ENDIF
  ** discrepancy with log-law values
(stored of UDIS is 100.*(U1/UIN-1.0) with IF(IX.LT.NX))
(stored of KDIS is 100.*(KE/KEIN-1.0))
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
RHO1=1.18; ENUL=1.E-5 ; PRT(C1)=1.
    GROUP 10. Inter-phase-transfer processes and properties
    GROUP 11. Initialization of variable or porosity fields
  ** Initial guess
REAL(GZM,GZP,GZ,GZN,CONLOG,GUI,GVI,GKI,GEPI,GOMI)
    GROUP 12. Patchwise adjustment of terms
    GROUP 12. Convection and diffusion adjustments
    GROUP 13. Boundary conditions and special sources
      Intialise field to inlet values
PATCH(INITBL,INIVAL,1,NX,1,NY,1,NZ,1,1)
(initial of U1 at INITBL is :USTARX:*LOGE(ZG/:Z0:)/:VKC:)
(initial of KE at INITBL is :KCNS:)
(initial of EP at INITBL is :EPCON:/ZG)
IF(KWMOD) THEN
+(initial of OMEG is :EPCON:/(ZG*0.09*:KCNS:))
ENDIF
  *** Inlet boundary; fully-developed turbulent flow 
PATCH(INBL,WEST,1,1,1,NY,1,NZ,1,1)
(source of P1 at INBL is COVAL(FIXFLU,RHO1*UIN))
(source of U1 at INBL is COVAL(ONLYMS,UIN))
COVAL(INBL,KE,ONLYMS,KCNS)
(source of EP at INBL is COVAL(ONLYMS,EPIN))
IF(KWMOD) THEN
+(source of OMEG at INBL is COVAL(ONLYMS,OMIN))
ENDIF
 
       **OUTLET EAST
PATCH(PEAST,EAST,NX,NX,1,NY,1,NZ,1,1)
COVAL(PEAST,P1,1000,0.)
       **SKY HIGH
REAL(USKY,MUTSKY,COSKY,HDZSKY,PRTKE,PRTEP)
GZN=ZWLAST/Z0;CONLOG=LOG(GZN)/VKC
USKY=USTARX*CONLOG;GEPI=EPCON/ZWLAST
 
   ** entrainment boundary condition
PATCH(SKY,HIGH,1,NX,1,NY,NZ,NZ,1,1)
COVAL(SKY,P1,1E3,0.);COVAL(SKY,U1,ONLYMS,USKY)
COVAL(SKY,KE,ONLYMS,KCNS);COVAL(SKY,EP,ONLYMS,GEPI)
IF(KWMOD) THEN
+ GOMI=GEPI/(0.09*KCNS)
+ COVAL(SKY,OMEG,ONLYMS,GOMI)
ENDIF
      ** diffusive boundary condition
MUTSKY=RHO1*VKC*USTAR*ZWLAST
 
PATCH(SKYD,HIGH,1,NX,1,NY,NZ,NZ,1,1)
(source of U1 at SKYD is COVAL(:MUTSKY:/(0.5*DZW),USKY))
PRTKE=PRT(KE)
(source of KE at SKYD is COVAL(:MUTSKY:/(:PRTKE:*0.5*DZW),KCNS))
IF(KWMOD) THEN
+PRTEP=PRT(OMEG)
+(source of OMEG at SKYD is COVAL(:MUTSKY:/(:PRTEP:*0.5*DZW),GOMI))
ELSE
+PRTEP=PRT(EP)
+(source of EP at SKYD is COVAL(:MUTSKY:/(:PRTEP:*0.5*DZW),GEPI))
ENDIF
    ** use built-in fully-rough wall functions
WALLCO=GRND5;WALLA=Z0
WALL(FRIC,LOW,1,NX,1,NY,1,1,1,1)
COVAL(FRIC,C1,GRND5,1.0)
    GROUP 14. Downstream pressure for PARAB=.TRUE.
    GROUP 15. Termination of sweeps
LSWEEP=500
    GROUP 16. Termination of iterations
LITER(P1)=50
    GROUP 17. Under-relaxation devices
DTF=XLEN/UIN
RELAX(U1,FALSDT,DTF);RELAX(W1,FALSDT,DTF)
RELAX(C1,LINRLX,0.5)
IF(KWMOD) THEN
+ RELAX(KE,FALSDT,DTF);RELAX(OMEG,FALSDT,DTF)
ELSE
+ RELAX(KE,LINRLX,0.5);RELAX(EP,LINRLX,0.5)
ENDIF
    GROUP 18. Limits on variables or increments to them
    GROUP 19. Data communicated by satellite to GROUND
    GROUP 20. Preliminary print-out
ECHO=T
    GROUP 21. Print-out of variables
    GROUP 22. Spot-value print-out
IXMON=11 ; IYMON=1 ; IZMON=11; IPLTL=LSWEEP
    GROUP 23. Field print-out and plot control
    GROUP 24. Dumps for restarts
SELREF=T ; RESFAC=1.E-4
TSTSWP=-1

DISTIL=T
CASE :CTURB: OF
WHEN KE,2
+ EX(P1  )=2.251E-03;EX(U1  )=2.758E+00
+ EX(W1  )=6.918E-04;EX(KE  )=2.663E-01
+ EX(EP  )=8.082E-03;EX(C1  )=8.385E-02
+ EX(KDIS)=3.454E-01;EX(UDIS)=2.488E-01
+ EX(EPIN)=7.922E-03;EX(KEIN)=2.657E-01
+ EX(UIN )=2.755E+00;EX(GEN1)=2.854E-02
+ EX(LEN1)=1.206E+01;EX(STRS)=1.568E-03
+ EX(ENUT)=3.409E+00
WHEN KO,2
+ EX(P1  )=2.763E-03;EX(U1  )=2.758E+00
+ EX(W1  )=7.699E-04;EX(KE  )=2.626E-01
+ EX(EP  )=7.911E-03;EX(C1  )=8.357E-02
+ EX(KDIS)=1.197E+00;EX(UDIS)=2.844E-01
+ EX(OMIN)=3.313E-01;EX(EPIN)=7.922E-03
+ EX(KEIN)=2.657E-01;EX(UIN )=2.755E+00
+ EX(GEN1)=2.846E-02;EX(LEN1)=1.196E+01
+ EX(STRS)=1.511E-03;EX(ENUT)=3.372E+00
+ EX(OMEG)=3.405E-01
WHEN RKE,3
+ EX(P1  )=2.002E-03;EX(U1  )=2.756E+00
+ EX(W1  )=5.624E-04; EX(KE  )=2.627E-01
+ EX(EP  )=8.024E-03; EX(C1  )=8.860E-02
+ EX(KDIS)=1.134E+00; EX(UDIS)=3.146E-01
+ EX(EPIN)=7.922E-03; EX(KEIN)=2.657E-01
+ EX(UIN )=2.755E+00; EX(GEN1)=2.689E-02
+ EX(LEN1)=1.194E+01; EX(STRS)=1.534E-03
+ EX(ENUT)=3.364E+00; EX(C1E )=4.300E-01
+ EX(DWDZ)=4.361E-05; EX(DWDX)=1.241E-05
+ EX(DUDZ)=9.440E-02; EX(DUDX)=3.960E-05
+ EX(EPKE)=3.102E-02;EX(CMU )=9.082E-02
ENDCASE