** PHOENICS VALIDATION CASE
TEXT(CHEN-KIM_2D ABRUPT PIPE EXPANS  :T213
TITLE
  DISPLAY
  The problem considered is the calculation of turbulent flow and
  heat transfer in a sudden pipe expansion. The geometry consists of
  an axisymmetric sudden pipe expansion of diameter ratio Do/Di=
  1.945. The Reynolds number based upon the bulk velocity and
  diameter of the larger pipe is 2.0E5. The thermal field is
  calculated for a constant heat flux through the side wall and an
  adiabatic back step. The laminar Prandtl number is taken as 1.0
  and the turbulent Prandtl number as 0.86.
 
  The calculation may be performed with one of the following low-Re
  turbulence models: the Lam-Bremhorst k-e model; the Chen-Kim k-e
  model; the 2-Layer k-e model; or the Wilcox k-omega model. The
  numerical integration is taken right down to the wall, although
  the mesh employed is still relatively coarse and so the solution
  cannot be considered as grid independent.
  ENDDIS
 
  The calculations are started at z= -4h and terminated at z=16h,
  where a fixed-pressure boundary condition is applied. Here, h is
  the step height given by h=0.5*(Do-Di). Uniform profiles of w, T,
  k and e are specified at the inlet to the calculation domain.
 
  The fine-mesh calculations employ a relatively coarse non-uniform
  mesh of NX=80 by NY=80 which concentrates grid cells close to the
  wall. About 2000 sweeps are required for a converged solution with
  this mesh. However, for speed of computation the calculation is
  performed on a coarse mesh of NX=25 by NY=20. Consequently, the
  solution is unlikely to be of sufficient quality because most of
  the near-wall cells are located too far from the wall for laminar
  conditions to apply on KE and EP. This solution accuracy will be
  particularly poor for the k-omega model, which requires fineg=t
  for an acceptable solution. PHOENICS will actually activate
  log-law wall functions for velocities if the near-wall point lies
  outside the viscous sublayer. For demonstration purposes the
  calculation is set up for 20 sweeps, but 250 sweeps are required
  for complete convergence.
 
    This AUTOPLOT sequence provides a plot of the axial
    velocity W1 along the outer surface of the solution
    domain versus normalised axial distance Z. The axial
    coordinate 0.0 corresponds to the step location. The
    reattachment point corresponds the axial location where
    W1 changes from negative to positive.
 
   AUTOPLOT USE
   file
   phi 5
 
 
   da 1 w1 y m
   shift z -.972 1
   divide z .243 1
   col9 1
   level y 0;level z 0
   msg Velocity (W1) profile
   msg Press e to END
   ENDUSE
 
CHAR(CTURB,TLSC);BOOLEAN(FINEG);INTEGER(TMODEL,NYS,NZS,NZS2)
REAL(HGHT,DIAMI,DIAMO,SLEN,SLEN2,REYNO,GRADI,RADO)
REAL(PLEN,WIN,WOUT,TKEIN,EPSIN,FLOW,AIN,QIN,AREAO)
FINEG=F
     ** Calculation of domain specifications
REYNO=2.E5;DIAMO=1.0;DIAMI=DIAMO/1.945;WIN=1.0
WOUT=WIN*(DIAMI/DIAMO)**2;ENUL=(WOUT*DIAMO)/REYNO
     ** Calculation of KE (where fric=0.018)...
WOUT
DIAMO
RADO=0.5*DIAMO;GRADI=0.5*DIAMI;HGHT=RADO-GRADI
SLEN=4.*HGHT;SLEN2=5.5*HGHT;PLEN=20.*HGHT
TKEIN=0.25*WIN*WIN*0.018;EPSIN=0.1643*TKEIN**1.5/(0.09*HGHT)
    GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1;AIN=DIAMI*DIAMI*XULAST/8.;FLOW=WIN*AIN
    GROUP 4. Y-direction grid specification
    GROUP 5. Z-direction grid specification
  ** Radius of small pipe = 0.257 Radius of large pipe  = 0.5
  ** Full length of channel = 4.86
IF(FINEG) THEN
+ NZS=14;NZS2=61;NYS=40
+ SUBGRD(Y,1,NYS,GRADI,-2.0);SUBGRD(Y,NY+1,-80,HGHT,1.7)
+ SUBGRD(Z,1,NZS,SLEN,-2.0);SUBGRD(Z,NZ+1,NZS2,SLEN2,2.0)
+ SUBGRD(Z,NZ+1,80,PLEN-SLEN2-SLEN,1.2)
ELSE
+ NZS=7;NZS2=17;NYS=10
+ SUBGRD(Y,1,NYS,GRADI,-2.0);SUBGRD(Y,NY+1,-20,HGHT,1.7)
+ SUBGRD(Z,1,NZS,SLEN,-2.0);SUBGRD(Z,NZ+1,NZS2,SLEN2,2.0)
+ SUBGRD(Z,NZ+1,25,PLEN-SLEN2-SLEN,1.2)
ENDIF
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
SOLVE(P1,W1,V1,H1);SOLUTN(P1,Y,Y,Y,N,N,N)
SOLUTN(W1,P,P,P,P,P,N);SOLUTN(V1,P,P,P,P,P,N)
SOLUTN(H1,Y,Y,Y,N,N,N);STORE(ENUT)
MESG( Enter the required turbulence model:
MESG(  LAMB - Low-Re Lam-Brem. k-e model
MESG(  CHEN - Low-Re Lam-Brem. k-e model (default)
MESG(  TWOL - Low-Re 2-layer   k-e model
MESG(  KO   - Low-Re k-omega model
MESG(
READVDU(CTURB,CHAR,CHEN)
CASE :CTURB: OF
WHEN LAMB,4
+ TEXT(LAM-BRE K-E_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re Lam-Bem. k-e turbulence model
+ TMODEL=1;KELIN=1;TLSC=EP
+ TURMOD(KEMODL-LOWRE);STORE(REYN)
WHEN CHEN,4
+ TEXT(CHEN-KIM_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re Chen-Kim k-e turbulence model
+ TMODEL=2;KELIN=1;TLSC=EP
+ TURMOD(KECHEN-LOWRE);STORE(REYN)
WHEN TWOL,4
+ TEXT(2-LAYER K-E_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re 2-layer k-e turbulence model
+ TMODEL=3;KELIN=1;TLSC=EP
+ TURMOD(KEMODL-2L);STORE(REYN)
WHEN KO,2
+ TEXT(K-OMEGA_2D ABRUPT PIPE EXPANS  :T213
+ MESG(Low-Re k-omega turbulence model
+ TMODEL=4;TLSC=OMEG
+ TURMOD(KOMODL-LOWRE);STORE(REYT)
+ STORE(EP);EPSIN=EPSIN/(0.09*TKEIN)
ENDCASE
    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,P,P,P,P,P)
    GROUP 9. Properties of the medium (or media)
RHO1=1.0;PRT(H1)=0.86
    GROUP 11. Initialization of variable or porosity fields
CONPOR(0.0,VOLUME,1,NX,-(NYS+1),-NY,-1,-NZS)
     ** Initial values
FIINIT(W1)=WIN;FIINIT(V1)=0.0;FIINIT(P1)=1.3E-4
FIINIT(KE)=TKEIN
FIINIT(EP)=(0.09*TKEIN*TKEIN)/(0.01*WOUT*0.5*DIAMO)
IF(TMODEL.EQ.4) THEN
+ FIINIT(OMEG)=TKEIN/(0.01*WOUT*0.5*DIAMO)
ELSE
+ FIINIT(EP)=(0.09*TKEIN*TKEIN)/(0.01*WOUT*0.5*DIAMO)
ENDIF
    GROUP 13. Boundary conditions and special sources
     ** Inlet
INLET(INLET,LOW,1,NX,1,NYS,1,1,1,1);VALUE(INLET,P1,RHO1*WIN)
VALUE(INLET,W1,WIN);VALUE(INLET,KE,TKEIN)
VALUE(INLET,:TLSC:,EPSIN);VALUE(INLET,H1,0.0)
     ** Exit
PATCH(OUTLET,HIGH,1,NX,1,NY,NZ,NZ,1,1);COVAL(OUTLET,P1,1.0E5,0.0)
COVAL(OUTLET,H1,ONLYMS,SAME)
     ** N-wall
WALL(WFUNNORT,NORTH,1,NX,NY,NY,NZS+1,NZ,1,1)
   ** Heat input; take cp=1 & Qin to raise temperature from 0 to 1.
QIN=RHO1*FLOW;AREAO=XULAST*RADO*(PLEN-SLEN)
PATCH(HEATIN,NORTH,1,NX,NY,NY,NZS+1,NZ,1,1)
COVAL(HEATIN,H1,FIXFLU,QIN/AREAO)
    GROUP 15. Termination of sweeps
LSWEEP=250;TSTSWP=-1
    GROUP 17. Under-relaxation devices
REAL(DTF,DTFKE);DTF=ZWLAST/NZ/WIN
RELAX(P1,LINRLX,1.0); RELAX(W1,FALSDT,DTF); RELAX(V1,FALSDT,DTF)
IF(TMODEL.EQ.4) THEN
+ RELAX(KE,FALSDT,DTF/5); RELAX(:TLSC:,FALSDT,DTF/5)
ELSE
+ RELAX(KE,FALSDT,DTF); RELAX(:TLSC:,FALSDT,DTF)
ENDIF
RELAX(H1,FALSDT,1.0)
    GROUP 18. Limits on variables or increments to them
   ** 2-layer k-e model does not work on PHOENICS V2.2
      when SELREF=T (default setting in PHOENICS)
IF(TMODEL.EQ.3) THEN
+ SELREF=F; RESREF(P1)=FLOW*1.E-12
+ RESREF(W1)=RHO1*RESREF(P1)*WIN; RESREF(V1)=RESREF(W1)
+ RESREF(KE)=RHO1*RESREF(P1)*TKEIN
+ RESREF(EP)=RHO1*RESREF(P1)*EPSIN; RESREF(H1)=QIN*1.E-12
ENDIF
+ VARMAX(H1)=1.E9;VARMIN(H1)=0.0;VARMIN(ENUT)=1.E-10
    GROUP 21. Print-out of variables
    GROUP 22. Monitor print-out
IZMON=NZS+2;IPLTL=2000;IYMON=NYS-2
ITABL=3;NPLT=10;NPRMON=10000;WALPRN=T