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 
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
   model.
 
   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 (ASME HTD-60, p75, New York, 1986) and numerically by
   Davidson (Num. Heat Transfer,Vol.18,p129,1990).
 
   The calculation may be performed with or without the
   Boussinesq approximation and with the energy equation solved via
   the temperature TEM1 or the static enthalpy H1.
   Both calculations compute the dynamic laminar viscosity from
   Sutherland's formula.
   When the Boussinesq approximation is not used, variable density 
   and thermal conductivity are also employed .
  ENDDIS
 
   Calculations with BOUSS=T require 1200 sweeps for convergence,
   whereas those with BOUSS=F require 1600 sweeps. For speed of
   computation during demonstration, the case has been set up to
   run for 100 sweeps only.
 
BOOLEAN(BOUSS); CHAR(CHOICE)
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
MESG( Enter required energy variable? (TEM1/H1)
MESG( Enter the required energy variable.
MESG(  TEM1 - Temperature ( default )
MESG(  H1   - Enthalpy
READVDU(CHOICE,CHAR,TEM1)
IF(:CHOICE:.EQ.TEM1) THEN
+ MESG( TEM1 solution selected
ELSE
+ MESG( H1 solution selected
+ CHOICE=H1
ENDIF
REAL(TCOLD,THOT,TDIF,CP,HEIGHT,WIDTH,MU,BETA,TREF,AGRAV)
REAL(HCOLD,HHOT,HREF,DTF,VBUOY,RAY,RHO,TBRUNT,KOND); CP=1.008E3
HEIGHT=2.5;WIDTH=0.5; TCOLD=34.2+273.; TDIF=45.8; THOT=TCOLD+TDIF
HCOLD=CP*TCOLD;HHOT=CP*THOT; TREF=0.5*(THOT+TCOLD);HREF=CP*TREF
AGRAV=-9.81;BETA=1./TREF; RHO=353.505/TREF;MU=2.0383E-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=28
REAL(ALF,GAM,ALFM);ALF=2.0;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=28
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
#solvel
SOLVE(:CHOICE:)
TURMOD(KEMODL-LOWRE)
STORE(ENUT,VISL)
IF(:CHOICE:.EQ.H1) THEN
+ STORE(TMP1); TMP1=LINH
ENDIF
    group 8. terms (in differential equations) & devices
TERMS(:CHOICE:,N,P,P,P,P,P)
    group 9. properties of the medium (or media)
KOND=3.5E-3+7.58E-5*TREF
IF(BOUSS) THEN
+ RHO1=RHO
+ PRNDTL(:CHOICE:)=-KOND
ELSE
+ STORE(DEN1,KOND)
+ FIINIT(DEN1)=RHO; FIINIT(KOND)=KOND
+ RHO1=GRND10; RHO1A=0.; RHO1B=1./353.505
+ PRNDTL(:CHOICE:)=-GRND4;PRLH1A=3.5E-3;PRLH1B=7.5E-5
ENDIF

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

RAY=-AGRAV*BETA*HEIGHT**3*TDIF*MU*CP/KOND/(MU/RHO)**2
RAY
BETA
    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(:CHOICE:.EQ.TEM1) THEN
+ HHOT=THOT; HCOLD=TCOLD; HREF=TREF
ENDIF
FIINIT(:CHOICE:)=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,:CHOICE:,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,:CHOICE:,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
PATCH(KEBUOY,PHASEM,1,NX,1,NY,1,NZ,1,1)
COVAL(KEBUOY,KE,KESOURCE,KESOURCE)
COVAL(KEBUOY,EP,KESOURCE,KESOURCE)
    GROUP 16. Termination of iterations
LITER(V1)=2;LITER(U1)=2;LITER(:CHOICE:)=50
LITER(KE)=2;LITER(EP)=2
    GROUP 17. Under-relaxation devices
VBUOY=(-BETA*TDIF*AGRAV*HEIGHT)**0.5; TBRUNT=HEIGHT/VBUOY
DTF=TBRUNT/NY;KELIN=3
IF(BOUSS) THEN
+ RELAX(V1,FALSDT,DTF); RELAX(U1,FALSDT,DTF)
+ RELAX(KE,FALSDT,DTF); RELAX(EP,FALSDT,DTF)
+ RELAX(:CHOICE:,FALSDT,50.*DTF);LSWEEP=1200
ELSE
+ RELAX(V1,FALSDT,DTF); RELAX(U1,FALSDT,DTF)
+ RELAX(KE,FALSDT,DTF); RELAX(EP,FALSDT,DTF)
+ RELAX(:CHOICE:,FALSDT,50.*DTF);LSWEEP=1600;DENPCO=T
ENDIF
        LSWEEP=100
    GROUP 18. Limits on variables or increments to them
VARMIN(:CHOICE:)=HCOLD; VARMAX(:CHOICE:)=HHOT
    GROUP 22. Spot-value print-out
IYMON=NY/2; IXMON=5; NYPRIN=NY/5; NXPRIN=NX/5; NPLT=20
TSTSWP=-1;YPLS=T;WALPRN=T; ITABL=3
    GROUP 23. Field print-out and plot control