TALK=T;RUN( 1, 1)
 ** LOAD(x310) from the x Input Library
TEXT(Realisable KE_2D Elliptic Plane Free Jet: T310
TITLE
  DISPLAY
  The problem considered is the submerged free turbulent plane
  jet in stagnant surroundings. The jet issues from a slot H
  at a Reynolds number of 14,000. The calculation exploits
  symmetry, and is carried out with the elliptic solver for a
  domain which is 30H long by 7.5H wide.
 
  Calculations are performed with the standard k-e & the Chen-Kim, 
  RNG and Realisable k-e variants. Simulations are also made with
  the 1988 & 2008 Wilcox k-w models, the Menter k-w model and the 
  k-w SST variant. The experimental data indicate a velocity half-
  width spreading rate of 0.11 in the self-similar region of
  the jet. The present calculations predict the following
  half-width spreading rates:
          data    KE    CK     RNG    RKE    KW    KWR   KWM   SST
   dy/dz  0.11  0.099  0.076  0.119  0.099  0.115 0.106 0.106 0.102
 
  The Realisable, RNG & standard k-e models produce fairly good 
  agreement with the data, and likewise with the Wilcox 2008, 
  Menter and k-w-SST models. However, the Chen-Kim k-e model 
  significantly underestimates the jet spreading, and the Wilcox 
  1988 k-w model with the default model constants produces an 
  excessive spreading rate. The value reported above for the 
  Wilcox 1988 model was obtained with a 30% reduction in the
  value of C2w in the w-equation. It should be mentioned that the 
  results reported here have not been subjected to a grid-sensitivity
  test. The half-width spreading rate is calculated approximately 
  using In-Form commands and then written to the text-output file 
  named inforout.
  ENDDIS
    GROUP 1. Run title and other preliminaries
REAL(WINJ,WIN_FS,KE_FS,EP_FS,KEINJ,EPINJ,HSLOT,PRADO,PRADI,CD)
REAL(OM_FS,OMINJ,LAMVIS,ENUT_FS,TINJ,TIN_FS,REYNO)
BOOLEAN(KWMOD);KWMOD=F
REYNO=1.4E4
HSLOT=0.025;PRADI=0.5*HSLOT;PRADO=15.*PRADI
CD=0.1643
 
  ** jet-inflow conditions
WINJ=10.;TINJ=0.05
KEINJ=(TINJ*WINJ)**2; EPINJ=CD*KEINJ**1.5/(0.1*PRADI)
   ** laminar kinematic viscosity
LAMVIS=WINJ*2.*HSLOT/REYNO
LAMVIS
   ** free-stream conditions
WIN_FS=WINJ/1000.;TIN_FS=0.05
ENUT_FS=LAMVIS
KE_FS=(TIN_FS*WIN_FS)**2;EP_FS=0.09*KE_FS**2/ENUT_FS
    GROUP 3. X-direction grid specification
    GROUP 4. Y-direction grid specification
INTEGER(NYF,NYO,NYG)
NYF=10;NYO=50;NYG=NYF+NYO
NREGY=2;NY=46
IREGY=1;GRDPWR(Y,NYF,PRADI,1.0)
IREGY=2;GRDPWR(Y,NYO,-(PRADO-PRADI),1.04)
 
    GROUP 5. Z-direction grid specification
NZ=120;GRDPWR(Z,NZ,-(30.*HSLOT),1.01)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);STORE(ENUT)
SOLUTN(P1,P,P,Y,P,P,P);SOLUTN(V1,P,P,P,P,P,N)
SOLUTN(W1,P,P,P,P,P,N)
CHAR(CTURB)
MESG( Enter the required turbulence model:
MESG(  KE   - Standard k-e model
MESG(  CK   - Chen Kim k-e model
MESG(  RNG  - RNG  k-e model
MESG(  RKE  - Realisable k-e model (default)
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,RKE)
CASE :CTURB: OF
WHEN KE,2
+ TEXT(Standard KE_2D Elliptic Plane Jet
+ MESG(Standard k-e model
+ TURMOD(KEMODL)
WHEN CK,2
+ TEXT(Chen-Kim KE_2D Elliptic Plane Jet
+ MESG(Chen Kim k-e model
+ TURMOD(KECHEN)
WHEN RNG,3
+ TEXT(RNG KE_2D Elliptic Plane Jet
+ MESG(RNG k-e model
+ TURMOD(KERNG)
WHEN RKE,3
+ TEXT(Realisable KE_2D Elliptic Plane Jet
+ MESG(RK k-e model
+ TURMOD(KEREAL)
+ STORE(CMU,C1E)
+ OUTPUT(CMU,P,P,P,P,Y,Y);OUTPUT(C1E,P,P,P,P,Y,Y)
WHEN KW,2
+ TEXT(Wilcox 1988 k-w 2D Elliptic Plane Jet
+ MESG(Wilcox 1988 k-w Model
+ TURMOD(KWMODL);KWMOD=T
+ OMINJ=EPINJ/(0.09*KEINJ)
+ ENUT_FS=0.001*LAMVIS
+ EP_FS=0.09*KE_FS**2/ENUT_FS
+ OM_FS=EP_FS/(0.09*KE_FS)
  ** lower C2f from 0.075 to match measured spreading rate
+ SPEDAT(KECONST,C2E,R,0.7*0.075)
WHEN KWR,3
+ TEXT(Wilcox 2008 k-w 2D Elliptic Plane Jet
+ MESG(Wilcox 2008 k-w Model
+ TURMOD(KWMODLR);KWMOD=T;FIINIT(FBP)=1.0
+ OMINJ=EPINJ/(0.09*KEINJ)
+ OM_FS=EP_FS/(0.09*KE_FS)
+ STORE(CDWS)
  ** deactivate vortex stretching term, as this is zero
     in 2d incompressible planar flow
+ SPEDAT(KWMOD,VORTSP,L,F)
WHEN KWM,3
TEXT(Menter k-w 2D Elliptic Plane Jet
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER);KWMOD=T
+ STORE(EP);OMINJ=EPINJ/(0.09*KEINJ)
+ OM_FS=EP_FS/(0.09*KE_FS)
+ STORE(BF1,CDWS,SIGK,SIGW)
WHEN KWS,3
TEXT(SST k-w 2D Elliptic Plane Jet
+ MESG(Menter 1992 k-w SST model
+ TURMOD(KWSST);KWMOD=T
+ STORE(EP);OMINJ=EPINJ/(0.09*KEINJ)
+ OM_FS=EP_FS/(0.09*KE_FS)
+ STORE(BF1,BF2,GEN1,CDWS,SIGK,SIGW)
ENDCASE
STORE(LEN1)
    GROUP 8. Terms (in differential equations) & devices
    GROUP 9. Properties of the medium (or media)
ENUL=LAMVIS
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WIN_FS
PATCH(INIT,INIVAL,1,NX,1,NYF,1,NZ,1,LSTEP)
INIT(INIT,W1,0.0,WINJ); FIINIT(KE)=KEINJ; FIINIT(EP)=EPINJ
    GROUP 13. Boundary conditions and special sources
  ** Jet Inlet Conditions
INLET(IN1,LOW,1,NX,1,NYF,1,1,1,LSTEP)
VALUE(IN1,P1,RHO1*WINJ); VALUE(IN1,W1,WINJ)
VALUE(IN1,KE,KEINJ);VALUE(IN1,EP,EPINJ)
IF(KWMOD) THEN
+VALUE(IN1,OMEG,OMINJ);FIINIT(OMEG)=OMINJ
ENDIF
  ** Free Boundary Conditions
INLET(IN2,LOW,1,NX,#2,#2,1,1,1,LSTEP)
VALUE(IN2,P1,RHO1*WIN_FS);VALUE(IN2,W1,WIN_FS)
VALUE(IN2,KE,KE_FS);VALUE(IN2,EP,EP_FS)
IF(KWMOD) THEN
+VALUE(IN2,OMEG,OM_FS)
ENDIF
 
PATCH(FREEB,NORTH,1,NX,NYG,NYG,1,NZ,1,LSTEP)
COVAL(FREEB,W1,ONLYMS,WIN_FS)
COVAL(FREEB,KE,ONLYMS,KE_FS);COVAL(FREEB,EP,ONLYMS,EP_FS)
COVAL(FREEB,P1,1.E4,0.0)
IF(KWMOD) THEN
+VALUE(FREEB,OMEG,OM_FS)
ENDIF
 
OUTLET(OUT,HIGH,1,NX,1,NYG,NZ,NZ,1,LSTEP)
COVAL(OUT,P1,1.E4,0.0)
VALUE(OUT,V1,0.0); VALUE(OUT,W1,0.0)
VALUE(OUT,KE,0.0);VALUE(OUT,EP,0.0)
IF(KWMOD) THEN
+VALUE(OUT,OMEG,OM_FS)
ENDIF
 
    GROUP 15. Termination of sweeps
LSWEEP=1000
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
IF(.NOT.KWMOD) THEN
+ KELIN=3
ENDIF	
REAL(RLXFAC); RLXFAC=8.*ZWLAST/WINJ/NZ
RELAX(V1,FALSDT,RLXFAC); RELAX(W1,FALSDT,RLXFAC)
RELAX(KE,LINRLX,0.4); RELAX(EP,LINRLX,0.4)
CASE :CTURB: OF
WHEN RKE,3
+ RELAX(KE,FALSDT,RLXFAC); RELAX(EP,FALSDT,RLXFAC)
+ RELAX(ENUT,LINRLX,0.2);VARMAX(ENUT)=0.1
WHEN KW,2
+RLXFAC=ZWLAST/WINJ/NZ
+RELAX(V1,FALSDT,4.*RLXFAC);RELAX(W1,FALSDT,4.*RLXFAC)
+RELAX(EP,LINRLX,1.0);RELAX(ENUT,LINRLX,0.25)
+RELAX(KE,FALSDT,4.*RLXFAC);RELAX(OMEG,FALSDT,4.*RLXFAC)
+VARMAX(ENUT)=0.2
+LSWEEP=1800
+OUTPUT(EP,P,P,P,P,Y,Y)
WHEN KWR,3
+RLXFAC=ZWLAST/WINJ/NZ
+RELAX(V1,FALSDT,4.*RLXFAC);RELAX(W1,FALSDT,4.*RLXFAC)
+RELAX(EP,LINRLX,1.0);RELAX(ENUT,LINRLX,0.25)
+RELAX(KE,FALSDT,4.*RLXFAC);RELAX(OMEG,FALSDT,4.*RLXFAC)
+RELAX(FBP,LINRLX,1.E-10)
+LSWEEP=1800
WHEN KWM,3
+RLXFAC=ZWLAST/WINJ/NZ
+RELAX(V1,FALSDT,RLXFAC);RELAX(W1,FALSDT,RLXFAC)
+RELAX(EP,LINRLX,1.0);RELAX(ENUT,LINRLX,1.0)
+RELAX(KE,FALSDT,RLXFAC);RELAX(OMEG,FALSDT,RLXFAC)
+LSWEEP=2500
+OUTPUT(EP,P,P,P,P,Y,Y)
WHEN KWS,3
+RLXFAC=ZWLAST/WINJ/NZ
+RELAX(V1,FALSDT,4.*RLXFAC);RELAX(W1,FALSDT,4.*RLXFAC)
+RELAX(EP,LINRLX,1.0);RELAX(ENUT,LINRLX,1.0)
+RELAX(KE,FALSDT,0.5*RLXFAC);RELAX(OMEG,FALSDT,0.5*RLXFAC)
+LSWEEP=2500
+OUTPUT(EP,P,P,P,P,Y,Y)
ENDCASE
    GROUP 18. Limits on variables or increments to them
RESFAC=1.E-4
    GROUP 20. Preliminary print-out
ECHO=T
    GROUP 21. Print-out of variables
NYPRIN=1
    GROUP 22. Spot-value print-out
TSTSWP=-1
IYMON=NYF+2;IZMON=NZ-1;NPLT=10;ITABL=3;NZPRIN=1
    GROUP 23. Field print-out and plot control
SPEDAT(SET,GXMONI,PLOTALL,L,T)
    GROUP 24. Dumps for restarts

	** compute jet half-width at each axial station
      --------------------------------------------
(stored of WH is 0.5*W1[&1&] )
(stored of YGP is YG)
PATCH(HWIDTH,CELL,1,NX,1,NY,1,NZ-1,1,LSTEP)
(stored of YH is 0.0)
(stored of YH at HWIDTH is YGP with IF(W1.GE.WH.AND.W1[,+1,].LE.WH))

   ** compute jet half-width spreading rate by sampling two axial stations
      --------------------------------------------------------------------
    ** first axial station
INTEGER(IZ1);IZ1=3*NZ/4
PATCH(HWIDTH1,CELL,1,NX,2,NY-1,IZ1,IZ1,1,LSTEP)
(make1 YH1)
(store1 of YH1 at HWIDTH1 is MAX(YH,1.1E-10))
(print YH1 is YH1)
    ** second axial station
PATCH(HWIDTH2,CELL,1,NX,2,NY-1,NZ-1,NZ-1,1,LSTEP)
(make1 YH2)
(store1 of YH2 at HWIDTH2 is MAX(YH,1.1E-10))
(print YH2 is YH2)
    ** jet half-width spreading rate
(make1 DYHDZ)
(stored of ZGM is ZG)
(store1 of DYHDZ is (YH2-YH1)/(ZGM[,1,NZ-1]-ZGM[,1,:IZ1:]))
(print DYHDZ IS DYHDZ)

   ** compute jet half width & store at jet axis for axial profile 
      plot in .csv file named IY1.CSV 
(make1 YH0D is 0)
(store1 of YH0D at HWIDTH is YGP with IF(W1.GE.WH.AND.W1[,+1,].LE.WH))
(stored of YH3D is YH0D)   
   
PATCH(IY1,PROFIL,1,1,1,1,1,NZ-1,1,1)
PLOT(IY1,W1,0.0,0.0);PLOT(IY1,KE,0.0,0.0);PLOT(IY1,YH3D,0.0,0.0)

DISTIL=T
CASE :CTURB: OF
WHEN KE,2
+EX(P1  )=7.548E-02;EX(V1  )=1.934E-01
+EX(W1  )=2.513E+00;EX(KE  )=1.000E+00
+EX(EP  )=5.415E+01;EX(YH3D)=2.818E-02
+EX(ZGM )=3.021E-01;EX(YH  )=4.607E-04
+EX(YGP )=6.195E-02;EX(WH  )=4.097E+00
+EX(EPKE)=2.713E+01;EX(LEN1)=4.238E-02
+EX(ENUT)=3.215E-03
WHEN CK,2
+EX(P1  )=5.890E-02;EX(V1  )=1.634E-01 
+EX(W1  )=2.472E+00;EX(KE  )=7.798E-01 
+EX(EP  )=4.749E+01;EX(ZGM )=3.021E-01 
+EX(YH  )=3.794E-04;EX(YGP )=6.195E-02 
+EX(WH  )=4.311E+00;EX(EPKE)=3.009E+01 
+EX(LEN1)=1.164E-02;EX(ENUT)=2.379E-03 
+EX(P1  )=5.636E-02;EX(V1  )=1.706E-01 
+EX(YH3D)=2.318E-02
WHEN RNG,3
+EX(W1  )=2.504E+00;EX(KE  )=9.638E-01
+EX(EP  )=4.911E+01;EX(ZGM )=3.021E-01
+EX(YH  )=4.484E-04;EX(YGP )=6.195E-02
+EX(WH  )=4.151E+00;EX(EPKE)=2.829E+01
+EX(LEN1)=4.415E-02;EX(ENUT)=3.016E-03
+EX(P1  )=7.015E-02;EX(V1  )=1.882E-01 
+EX(YH3D)=2.740E-02
WHEN RKE,3
+EX(P1  )=8.042E-02;EX(V1  )=1.931E-01
+EX(W1  )=2.548E+00;EX(KE  )=1.038E+00 
+EX(EP  )=5.927E+01;EX(ZGM )=3.021E-01 
+EX(YH  )=4.690E-04;EX(YGP )=6.195E-02 
+EX(WH  )=4.051E+00;EX(LEN1)=4.810E-02 
+EX(C1E )=6.296E-01;EX(DWDZ)=5.430E+00 
+EX(DWDY)=8.830E+01;EX(DVDZ)=1.440E+00 
+EX(DVDY)=5.612E+00;EX(EPKE)=2.701E+01 
+EX(CMU )=6.639E-02;EX(ENUT)=3.312E-03 
+EX(YH3D)=2.861E-02
WHEN KW,2
+EX(P1  )=3.183E-01;EX(V1  )=1.864E-01 
+EX(W1  )=2.012E+00;EX(KE  )=4.317E+00 
+EX(EP  )=2.410E+02;EX(YH3D)=7.309E-02 
+EX(ZGM )=3.021E-01;EX(YH  )=1.200E-03 
+EX(YGP )=6.195E-02;EX(WH  )=1.889E+00 
+EX(LEN1)=1.939E-02;EX(OMEG)=2.035E+02 
+EX(ENUT)=1.444E-02 
WHEN KWR,3
+EX(P1  )=1.336E-01;EX(V1  )=2.266E-01
+EX(W1  )=2.679E+00;EX(KE  )=1.494E+00
+EX(EP  )=7.109E+01;EX(XWP )=1.000E-10
+EX(ZGM )=3.021E-01;EX(YH  )=5.939E-04
+EX(YGP )=6.195E-02;EX(WH  )=3.669E+00
+EX(LEN1)=4.484E-02;EX(YH3D)=3.629E-02
+EX(CDWS)=1.069E+04;EX(DWDZ)=6.647E+00
+EX(DWDY)=7.429E+01;EX(DVDZ)=2.022E+00
+EX(DVDY)=6.713E+00;EX(GEN1)=2.560E+04
+EX(OMEG)=2.395E+02
+EX(ENUT)=5.027E-03;EX(FBP )=1.000E+00
WHEN KWM,3
+EX(ZGM )=3.021E-01
+EX(YH  )=5.835E-04;EX(YGP )=6.195E-02
+EX(EP  )=6.861E+01;EX(WH  )=3.764E+00 
+EX(SIGW)=1.168E+00;EX(SIGK)=1.000E+00
+EX(P1  )=1.045E-01;EX(V1  )=2.149E-01 
+EX(W1  )=2.613E+00;EX(KE  )=1.265E+00 
+EX(LEN1)=6.156E-02;EX(CDWS)=8.810E+03 
+EX(OMEG)=2.739E+02;EX(ENUT)=4.786E-03 
+EX(YH3D)=3.469E-02;EX(YH  )=5.673E-04
WHEN KWS,3
+EX(EP  )=7.554E+01;EX(ZGM )=3.021E-01
+EX(YGP )=6.195E-02;EX(ENUT)=2.507E-02 
+EX(SIGW)=1.168E+00;EX(SIGK)=1.000E+00 
+EX(P1  )=1.027E-01;EX(V1  )=2.140E-01
+EX(W1  )=2.609E+00;EX(KE  )=1.236E+00
+EX(EP  )=6.726E+01;EX(YH  )=5.669E-04
+EX(WH  )=3.795E+00;EX(LEN1)=1.519E-02
+EX(CDWS)=8.912E+03;EX(GEN1)=2.684E+04
+EX(OMEG)=2.773E+02;EX(ENUT)=4.012E-03
+EX(LEN1)=1.041E-02;EX(YH3D)=3.425E-02
ENDCASE
  LIBREF = 310
STOP