TALK=T;RUN( 1, 1)
 PHOTON USE
   p;;;;
 
   view z; Gr ou z 1
   msg velocity vectors
   vec z 1 sh; pause
   msg wall-distance contours
   con wdis z 1 fi;0.001; Pause
   msg effective-viscosity contours
   con enut z 1 fi;0.001; Pause
 
   ENDUSE
  GROUP 1. RUN TITLE AND OTHER PRELIMINARIES
TEXT(2D Turbulent Free Convection In A Cavity: T219
TITLE
  DISPLAY
   Turbulent buoyancy-driven air flow in a tall rectangular
   cavity of 5:1 aspect ratio with a Rayleigh number of 4.E10 is
   calculated using a buoyancy-extended low-Reynolds-number k-e
   and k-w models. The vertical walls are isothermal with a 
   temperature difference of 45.8 degC; the horizontal walls are 
   adiabatic.
 
   This problem has been studied experimentally by Cheesewright
   et al (1986) and numerically by Davidson (1990) and Peng &
   Davidson (1999). The primary objective is to compute the
   Nusselt number distribution along the hot and cold walls, and
   the overall Nusselt number. The calculation may be performed 
   with or without the Boussinesq approximation, and the dynamic 
   laminar viscosity is computed from Sutherland's formula. For 
   the Boussinesq approximation, a uniform density and thermal 
   conductivity are employed in the simulations.
   
   The main result is the average Nusselt number, and the table
   below compares the predicted values with the value deduced
   from the measurements:
                Data      LB     KW-R    KW-SST
       Nu,av     160      154	  160     153
	   
   All of the models produce good agreement with the measured 
   value.	   

   ENDDIS
  
   Davidson, L., "Calculation of the turbulent buoyancy-driven
   flow in a rectangular cavity using an efficient solver and two
   different low-Reynolds number k-e models" Numerical Heat  
   Transfer, Vol.18, p129-147, (1990).
  
   Peng,S.H. & Davidson,L.,"Computation of turbulent buoyant flows 
   in enclosures with low-Reynolds-number k-w models", Int.J. Heat 
   & Fluid Flow, Vol.20, p172-184,(1999).
   
   Cheesewright,R., King,K.J. & Ziai,S., "Experimental data for the
   validation of computer codes for the prediction of two-dimensional
   buoyant cavity Flows", In: J.A.C. Humphrey, C.T. Avedisian, B.W.
   Le Tourneau, M.M. Chen (Eds.), Significant Questions in Buoyancy
   Affected Enclosure or Cavity Flows, HTD 60, ASME, p75-81,(1986).
   
BOOLEAN(BOUSS,KO);KO=F
MESG( Do you want the Boussinesq approximation? (y/N)
READVDU(ANS,CHAR,N)
IF(:ANS:.EQ.Y) THEN
+ BOUSS=T
+ MESG( Boussinesq approximation
ELSE
+ BOUSS=F
+ MESG( Variable density option
ENDIF

REAL(TCOLD,THOT,TDIF,CP,HEIGHT,WIDTH,MU,BETA,TREF,AGRAV)
REAL(HCOLD,HHOT,HREF,DTF,VBUOY,RAY,RHO,TBRUNT,COND,GRAS)
REAL(PRLAM,ENULAM,GRSHOF,VREF) 

   ** Cavity dimensions
HEIGHT=2.5;WIDTH=0.5
   ** Temperature levels
TCOLD=34.2+273.; TDIF=45.8; THOT=TCOLD+TDIF
TREF=0.5*(THOT+TCOLD)
   ** Physical properties
CP=1.008E3;BETA=1./TREF; RHO=353.505/TREF;MU=2.0383E-5
ENULAM=MU/RHO
   ** Enthalpies
HCOLD=CP*TCOLD;HHOT=CP*THOT;HREF=CP*TREF

AGRAV=-9.81

VREF=(-AGRAV*BETA*TDIF*HEIGHT)**0.5

    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
   ** for a fine mesh use NX=NY=56 and ALF=3.5
NX=56
REAL(ALF,GAM,ALFM);ALF=3.5;ALFM=-ALF
DO JJ=1,NX
+ GAM=ALF*(2.*(JJ/NX)-1.)
+ XFRAC(JJ)=-0.5*TANH(GAM)/TANH(ALFM)+0.5
+ XFRAC(JJ)=XFRAC(JJ)*WIDTH
ENDDO
    group 4. y-direction grid specification
NY=56
DO JJ=1,NY
+ GAM=ALF*(2.*(JJ/NY)-1.)
+ YFRAC(JJ)=-0.5*TANH(GAM)/TANH(ALFM)+0.5
+ YFRAC(JJ)=YFRAC(JJ)*HEIGHT
ENDDO
    group 7. variables stored, solved & named
SOLVE(P1,U1,V1,TEM1)

REAL(FRIC,TKEIN,UMAX,MIXL,EPSIN)
UMAX=0.4
FRIC=0.001;TKEIN=0.25*UMAX*UMAX*FRIC;MIXL=0.09*WIDTH
CHAR(CTURB,TLSC)
MESG( Enter the required turbulence model:
MESG(  LB - Low-Re Lam-Brem. k-e model (default)
MESG(  KWR  - Wilcox 2008 Low-Re k-w model
MESG(  KWS  - Low-Re k-w SST model
MESG(
READVDU(CTURB,CHAR,LB)
CASE :CTURB: OF
WHEN LB,2
+ TEXT(2D Turb. Buoyant Cavity Flow K-E:T219 
+ MESG(Low-Re Lam-Bem. k-e model
+ KELIN=1;TLSC=EP;KELIN=3
+ TURMOD(KEMODL-LOWRE);STORE(REYN)
WHEN KWR,3
+ TEXT(2D Turb. Buoyant Cavity Flow K-W R:T219 
+ KO=T
+ MESG(Wilcox 2008 low-Re k-w model
+ TLSC=OMEG
+ TURMOD(KWMODLR-LOWRE)
+ EPSIN=EPSIN/(0.09*TKEIN)
WHEN KWS,3
+ TEXT(2D Turb. Buoyant Cavity Flow K-W SST:T219 
+ KO=T
+ MESG(Menter k-w SST model
+ TURMOD(KWSST-LOWRE)
+ TLSC=OMEG
+ EPSIN=EPSIN/(0.09*TKEIN)
+ FIINIT(BF1)=1.0;FIINIT(BF2)=1.0
+ RELAX(BF1,LINRLX,0.05);RELAX(BF2,LINRLX,1.0)
ENDCASE
STORE(TREY) ! local turbulent Reynolds number
    group 8. terms (in differential equations) & devices
TERMS(TEM1,N,P,P,P,P,P)
    group 9. properties of the medium (or media)
COND=3.5E-3+7.58E-5*TREF
IF(BOUSS) THEN
+ RHO1=RHO
+ PRNDTL(TEM1)=-KOND
ELSE
+ STORE(DEN1,KOND)
+ FIINIT(DEN1)=RHO; FIINIT(KOND)=COND
+ RHO1=GRND10; RHO1A=0.; RHO1B=1./353.505
+ PRNDTL(TEM1)=-GRND4;PRLH1A=3.5E-3;PRLH1B=7.5E-5
ENDIF

CP1 =CP
ENUL=SUTHLAND;ENULA=1.458E-6;ENULB=110.4
PRT(TEM1)=0.9

PRLAM=MU*CP/COND
PRLAM
GRSHOF=-AGRAV*BETA*HEIGHT**3*TDIF/(ENULAM)**2
GRSHOF
RAY=GRSHOF*PRLAM
RAY
BETA
COND
STORE(ENUT,ENUL,YPLS,LEN1)
STORE(DWAL,QH,QC,HTH,HTC,NUH,NUC)
  __________________________________________________________________
   SAVE7BEGIN
(STORED of DWAL is 0.5*DXU)
   ** wall heat flux in W/m^2
(STORED QH at HOT is KOND*(THOT-TEM1)/DWAL)
(STORED QC at COLD is KOND*(TCOLD-TEM1)/DWAL)
   ** local heat transfer coefficient 
(STORED HTH at HOT is QH/TDIF)
(STORED HTC at COLD is QC/TDIF)
   ** Local Nusselt number 
      Nu = h_loc*H/k = (q_loc/DT)*H/k = k*((Tw-Tp)/dx)*H/(k*DT)
	     = ((Tw-Tp)/dx)*H/DT	  
(STORED NUH at HOT is ((THOT-TEM1)/DWAL)*:HEIGHT:/:TDIF:)
(STORED NUC at COLD is ABS((TCOLD-TEM1)/DWAL)*:HEIGHT:/:TDIF:)

   ** compute overall heat transfer rates to compare
      with NETT source printout in RESULT file
(MAKE1 QHOT is 0.0)
(STORE1 QHOT at HOT is SUM(QH*AEAST))
(PRINT Q_HOT is QHOT)

(MAKE1 QCOLD is 0.0)
(STORE1 QCOLD at COLD is SUM(QC*AEAST))
(PRINT of Q_Cold is QCOLD)

   ** Compute average Nusselt number
(MAKE1 NUAHOT is 0.0)
(STORE1 of NUAHOT at HOT is QHOT/(:TDIF:*:COND:))     
(PRINT1 of Nu_av_Hot is NUAHOT)
 
(MAKE1 NUACOLD is 0.0)
(STORE1 of NUACOLD at COLD is ABS(QCOLD)/(:TDIF:*:COND:))     
(PRINT1 of Nu_av_Cold is NUACOLD)

   ** dimensionless velocities
STORE(TNO,U1N,V1N)
(STORED OF U1N is U1/:VREF:)  
(STORED OF V1N is V1/:VREF:)  
   ** dimensionless temperature
(STORED of TNO IS (TEM1-:TCOLD:)/:TDIF:)
  SAVE7END
    _____________________________________________________________________

    GROUP 11. Initialization of variable or porosity fields
FIINIT(ENUT)=40.*MU/RHO
VBUOY=(-2.*BETA*0.5*TDIF*AGRAV*0.1*HEIGHT)**0.5
VBUOY
FIINIT(KE)=(0.1*VBUOY)**2
FIINIT(EP)=0.09*FIINIT(KE)**2/FIINIT(ENUT)
IF(TEM1.EQ.TEM1) THEN
+ HHOT=THOT; HCOLD=TCOLD; HREF=TREF
ENDIF
FIINIT(TEM1)=0.5*(HHOT+HCOLD)
    group 13. boundary conditions and special sources
   1. Hot wall boundary: Constant Temperature of 65.8 deg C.
WALL (HOT,WEST,1,1,1,NY,1,NZ,1,1)
COVAL(HOT,TEM1,LOGLAW,HHOT)
   2. Cold wall boundary: Constant Temperature of 20 deg C.
WALL (COLD,EAST,NX,NX,1,NY,1,NZ,1,1)
COVAL(COLD,TEM1,LOGLAW,HCOLD)
   3. Top wall boundary: adiabatic but with friction
WALL(TOPWAL,NORTH,1,NX,NY,NY,1,NZ,1,1)
   4. Bottom wall boundary: adiabatic but with friction
WALL(BOTWAL,SOUTH,1,NX,1,1,1,NZ,1,1)
   5. Buoyancy force
BUOYA=0.;BUOYB=AGRAV;BUOYC=0.
PATCH(BUOY,PHASEM,1,NX,1,NY,1,NZ,1,1)
IF(BOUSS) THEN
+ COVAL(BUOY,V1,FIXFLU,BOUSS)
+ DVO1DT=BETA; BUOYE=HREF
ELSE
+ COVAL(BUOY,V1,FIXFLU,DENSDIFF);BUOYD=RHO
ENDIF
   6. Reference pressure at the centre of the cavity
PATCH(REFP,CELL,NX/2,NX/2,NY/2,NY/2,1,NZ,1,1)
COVAL(REFP,P1,FIXP,0.0); COVAL(REFP,V1,ONLYMS,0.0)
COVAL(REFP,H1,ONLYMS,SAME)

   7. Buoyancy terms in the k-e model 
SPEDAT(SET,KEBUOY,GCEB,R,1.0)  ! C3EB = 1.0 by default   
    ** c3eb<0 means C3EB = tanh(V/U)
SPEDAT(SET,KEBUOY,GCEB,R,-1.0) 
PATCH(KEBUOY,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(KEBUOY,KE,KESOURCE,KESOURCE)
COVAL(KEBUOY,:TLSC:,KESOURCE,KESOURCE)
STORE(GENB,C3EB) ! Output of buoyancy production & C3EB
    GROUP 16. Termination of iterations
LITER(P1)=40	
LITER(V1)=2;LITER(U1)=2;LITER(TEM1)=50
LITER(KE)=2;LITER(:TLSC:)=2
    GROUP 17. Under-relaxation devices
VBUOY=(-BETA*TDIF*AGRAV*HEIGHT)**0.5; TBRUNT=HEIGHT/VBUOY
DTF=TBRUNT/NY
IF(BOUSS) THEN
+ RELAX(V1,FALSDT,DTF); RELAX(U1,FALSDT,DTF)
+ RELAX(KE,FALSDT,DTF*5); RELAX(:TLSC:,FALSDT,DTF)
ELSE
+ RELAX(V1,FALSDT,DTF); RELAX(U1,FALSDT,DTF)
+ RELAX(KE,FALSDT,DTF); RELAX(:TLSC:,FALSDT,DTF)
+ DENPCO=T
ENDIF
RELAX(TEM1,FALSDT,50.*DTF)
   RELAX(KE,FALSDT,DTF); RELAX(:TLSC:,FALSDT,DTF)
LSWEEP=8000
    GROUP 18. Limits on variables or increments to them
VARMIN(TEM1)=HCOLD; VARMAX(TEM1)=HHOT
    GROUP 22. Spot-value print-out
IYMON=NY/2; IXMON=5; NYPRIN=NY/5; NXPRIN=NX/5; NPLT=LSWEEP/20
TSTSWP=-1; ITABL=3
    GROUP 23. Field print-out and plot control
   ** Output Hot & Cold wall Nusselt numbers to .csv file	
PATCH(NUHOT,PROFIL,1,1,1,NY,1,NZ,1,LSTEP)
COVAL(NUHOT,NUH,0.0,0.0)	
PATCH(NUCOLD,PROFIL,NX,NX,1,NY,1,NZ,1,LSTEP)
COVAL(NUCOLD,NUC,0.0,0.0)	

SPEDAT(SET,GXMONI,PLOTALL,L,T)
DISTIL=T 

CASE :CTURB: OF
WHEN LB,2
+EX(P1  )=2.086E-01;EX(U1  )=8.015E-03 
+EX(V1  )=5.177E-02;EX(KE  )=6.196E-04 
+EX(EP  )=8.658E-04;EX(EPKE)=1.422E+06 
+EX(NUC )=8.657E-01;EX(NUH )=7.796E-01 
+EX(HTC )=4.602E-02;EX(HTH )=4.668E-02 
+EX(QC  )=2.108E+00;EX(QH  )=2.138E+00 
+EX(DWAL)=4.464E-03;EX(ENUL)=1.858E-05 
+EX(ENUT)=6.212E-05;EX(KOND)=2.827E-02 
+EX(DEN1)=1.073E+00;EX(REYN)=2.763E+01 
+EX(LTLS)=4.122E-03;EX(WDIS)=2.636E-02 
+EX(TNO )=5.034E-01;EX(TEM1)=3.303E+02 
+EX(YPLS)=1.008E-02;EX(V1N )=2.757E-02
+EX(U1N )=4.267E-03;EX(NUC )=4.328E+00 
+EX(NUH )=3.898E+00;EX(TREY)=3.357E+00 
+EX(C3EB)=6.060E-01;EX(GENB)=2.424E-05 
+EX(ENUT)=5.987E-05;EX(TREY)=3.235E+00 
 EX(LEN1)=4.112E-03
WHEN KWR,3
+EX(P1  )=1.636E-01;EX(U1  )=1.111E-02
+EX(V1  )=6.188E-02;EX(KE  )=5.287E-04
+EX(EP  )=9.268E-04;EX(NUC )=8.318E-01
+EX(NUH )=7.490E-01;EX(HTC )=4.422E-02
+EX(HTH )=4.485E-02;EX(QC  )=2.025E+00
+EX(QH  )=2.054E+00;EX(DWAL)=4.464E-03
+EX(YPLS)=1.075E-02;EX(ENUL)=1.857E-05
+EX(ENUT)=4.976E-05;EX(KOND)=2.827E-02
+EX(DEN1)=1.073E+00
+EX(DVDY)=6.471E-01;EX(DVDX)=1.884E+01
+EX(DUDY)=1.084E+00;EX(DUDX)=6.311E-01
+EX(GEN1)=1.566E+03;EX(FBP )=9.833E-01
+EX(XWP )=1.336E-02;EX(OMEG)=1.638E+04
+EX(TNO )=5.023E-01;EX(V1N )=3.295E-02
+EX(U1N )=5.916E-03;EX(TEM1)=3.302E+02
+EX(C3EB)=6.490E-01;EX(GENB)=1.958E-05
+EX(NUC )=4.159E+00;EX(NUH )=3.745E+00
+EX(LEN1)=2.672E-03;EX(TREY)=2.688E+00
+EX(P1  )=1.600E-01;EX(KE  )=5.638E-04
+EX(GENB)=2.435E-05;EX(LEN1)=4.097E-03
+EX(ENUT)=6.532E-05;EX(TREY)=3.528E+00
+EX(XWP )=2.916E-02
WHEN KWS,3
+EX(P1  )=2.145E-01;EX(U1  )=9.293E-03
+EX(V1  )=5.816E-02;EX(KE  )=3.490E-04
+EX(EP  )=7.322E-04;EX(NUC )=8.481E-01
+EX(NUH )=7.638E-01;EX(HTC )=4.508E-02
+EX(HTH )=4.573E-02;EX(QC  )=2.065E+00
+EX(QH  )=2.095E+00;EX(DWAL)=4.464E-03
+EX(YPLS)=1.020E-02;EX(ENUL)=1.858E-05
+EX(ENUT)=2.441E-05;EX(KOND)=2.827E-02
+EX(DEN1)=1.073E+00;EX(LTLS)=4.122E-03
+EX(WDIS)=2.636E-02;EX(GEN1)=1.706E+03
+EX(BF2 )=9.526E-01;EX(BF1 )=8.270E-01
+EX(OMEG)=1.542E+04;EX(TNO )=5.033E-01
+EX(V1N )=3.097E-02;EX(U1N )=4.948E-03
+EX(TEM1)=3.302E+02;EX(C3EB)=6.160E-01
+EX(GENB)=1.153E-05;EX(NUC )=4.240E+00
+EX(NUH )=3.819E+00;EX(LEN1)=1.616E-03
+EX(TREY)=1.319E+00
ENDCASE
STOP