TALK=T;RUN( 1, 1)
** LOAD(x306) from the x Input Library
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 leaves the solution domain
at the east boundary 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 to within a few percent and returns a
uniform pressure field. The case can also be simulated by using
the Realisable k-e model, the Wilcox 1988 & 2008 k-w models, and
the Menter k-w & k-w SST models.
ENDDIS
Notes:
1) The inlet log-law profiles aren't preserved exactly because the
analytical profiles aren't represented exactly in the discrete
form of the finite-volume equations. The solutions show a spike
in the k profile at the second cell from the wall. This originates
from the discretisation practices employed for the turbulence
production rate Pk in the turbulence transport equations and the
vertical diffusion terms in the axial momentum equation.
2) Specifically, the momentum equation uses face-centred differencing
for the vertical diffusion terms, which means the shearing stresses
in the momentum equation also use face-centred differencing. However,
the production rate Pk uses cell-centred differencing for the
velocity gradient associated with mean shear (together with a
cell-centred value for the turbulent viscosity ). This is a standard
approach used in CFD codes for computational convenience when dealing
with generic applications, but as shown by Norris & Richards 2010) it
fails to produce a uniform k profile for the limiting case of the
constant-stress layer.
3) The discretisation practices of Norris and reynolds (2010) can be
adopted to eliminate the anomoly, but given the relative complexity
of this practice coupled with increased computer time, it isn't worthwhile
because the velocity field is well predicted to within a few percent of
the analytical values; and also the resulting disturbances to what
should be a uniform pressure field are insignificant compared to the
axial momentum flux.
4)Norris, S.E. & Richards, P.J.,"Appropriate boundary conditions for
computational wind engineering models revisited", 5th Int. Symp. on
Computational Wind Engineering (CWE2010), Chapel Hill, North Carolina,
USA May 23-27, (2010).
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( RKE - Realisable k-e model
MESG( KW - Wilcox 1988 k-w model
MESG( KWR - Wilcox 2008 k-w model
MESG( KWM - Menter 1992 k-w model
MESG( KWS - k-w SST 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 RKE,3
+ TEXT(2D RK-E ATMOSPHERIC BOUNDARY LAYER: T306
+ MESG(Realisable k-e model
+ TURMOD(KEREAL);STORE(C1E)
WHEN KW,2
+ TEXT(2D Wilcox 1988 k-w ATMOSPHERIC BOUNDARY L: T306
+ MESG(Wilcox 1988 k-w model
+ TURMOD(KWMODL);KWMOD=T
WHEN KWR,3
+ TEXT(2D Wilcox 2008 k-w ATMOSPHERIC BOUNDARY L: T306
+ MESG(Wilcox 2008 k-w model
+ TURMOD(KWMODL);KWMOD=T
WHEN KWM,3
TEXT(Menter k-w ATMOSPHERIC BOUNDARY LAYER: T306
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER)
+ STORE(EP)
+ STORE(BF1);FIINIT(BF1)=1.0
+ KWMOD=T
WHEN KWS,3
TEXT(SST k-w ATMOSPHERIC BOUNDARY LAYER: T306
+ MESG(Menter 1992 k-w SST model
+ TURMOD(KWSST)
+ STORE(EP)
+ STORE(BF1,BF2,GEN1);FIINIT(BF1)=1.0
+ FIINIT(BF2)=1.0
+ KWMOD=T
+(stored of S is GEN1^0.5)
+(stored of SANY is :ustar:/(:vkc:*ZG) )
+ STORE(DUDZ)
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))
CASE :CTURB: OF
WHEN KW,2
+PRTKE=PRT(KE);PRTEP=PRT(OMEG)
WHEN KWR,3
+PRTKE=PRT(KE);PRTEP=PRT(OMEG)
WHEN KWM,3
+PRTKE=2.0;PRTEP=2.0
WHEN KWS,3
+PRTKE=2.0;PRTEP=2.0
ENDCASE
IF(KWMOD) THEN
+(source of KE at SKYD is COVAL(:MUTSKY:/(:PRTKE:*0.5*DZW),KCNS))
+(source of OMEG at SKYD is COVAL(:MUTSKY:/(:PRTEP:*0.5*DZW),GOMI))
ELSE
+PRTKE=PRT(KE)
+(source of KE at SKYD is COVAL(:MUTSKY:/(:PRTKE:*0.5*DZW),KCNS))
+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)
IF(KWMOD) THEN
+ COVAL(FRIC,OMEG,GRND5,GRND5)
ENDIF
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 RKE,3
+ EX(P1 )=1.822E-03;EX(U1 )=2.756E+00
+ EX(W1 )=5.003E-04;EX(KE )=2.627E-01
+ EX(EP )=8.024E-03;EX(C1 )=8.860E-02
+ EX(KDIS)=1.111E+00;EX(UDIS)=3.414E-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.556E-05;EX(DWDX)=1.324E-05
+ EX(DUDZ)=9.440E-02;EX(DUDX)=4.080E-05
+ EX(EPKE)=3.102E-02;EX(CMU )=9.082E-02
WHEN KW,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 KWR,3
+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(OMIN)=3.313E-01;EX(KDIS)=1.197E+00
+EX(UDIS)=2.844E-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 KWM,3
+EX(P1 )=2.708E-03;EX(W1 )=7.539E-04
+EX(UDIS)=2.782E-01
+EX(U1 )=2.758E+00;EX(KE )=2.627E-01
+EX(EP )=7.927E-03;EX(C1 )=8.360E-02
+EX(OMIN)=3.313E-01;EX(KDIS)=1.129E+00
+EX(UDIS)=2.782E-01;EX(EPIN)=7.922E-03
+EX(KEIN)=2.657E-01;EX(UIN )=2.755E+00
+EX(GEN1)=2.846E-02;EX(LEN1)=1.180E+01
+EX(STRS)=1.515E-03;EX(ENUT)=3.327E+00
+EX(BF1 )=1.000E+00;EX(OMEG)=3.405E-01
+EX(LTLS)=2.133E+03;EX(WDIS)=2.921E+01
WHEN KWS,3
+EX(P1 )=4.209E-03;EX(U1 )=2.760E+00
+EX(W1 )=1.463E-03;EX(EP )=7.847E-03
+EX(C1 )=8.286E-02;EX(KDIS)=1.601E+00
+EX(UDIS)=5.000E-01;EX(STRS)=1.494E-03
+EX(ENUT)=3.356E+00;EX(GEN1)=2.837E-02
+EX(OMEG)=3.381E-01;EX(KE )=2.624E-01
+EX(OMIN)=3.313E-01;EX(EPIN)=7.922E-03
+EX(KEIN)=2.657E-01;EX(UIN )=2.755E+00
+EX(LEN1)=1.197E+01;EX(BF2 )=1.000E+00
+EX(BF1 )=1.000E+00;EX(LTLS)=2.133E+03
+EX(WDIS)=2.921E+01;EX(DUDZ)=9.441E-02
+EX(SANY)=9.940E-02;EX(S )=9.243E-02
ENDCASE
LIBREF = 306
STOP