TALK=T;RUN( 1, 1)
** LOAD(x309) from the x Input Library
TEXT(Realisable KE_2D Elliptic Round Free Jet: T309
TITLE
DISPLAY
The problem considered is the submerged free turbulent round
jet in stagnant surroundings. The jet issues from a pipe of
diameter D at a Reynolds number of 38,500. The calculation is
carried out with the elliptic solver for a domain which is 30D
long by 7.5D wide.
Calculations are performed with the standard, 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.086 in the self-similar region of the jet.
The present calculations predict the following spreading rates:
data KE CK RNG RKE KW KWR KWM KWS
dyh/dz 0.086 0.111 0.107 0.167 0.10 0.077 0.082 0.119 0.119
Apart from the Wilcox 2008 k-w model, all other models seriously
overestimate the measured spreading rate, and especially the
standard k-w model and the RNG k-e model. The realisable k-e model
is reported to predict the correct spreading rate, but here the
improvement over the standard k-e model isn't remarkable. This
needs further investigation. The value reported above for the
Wilcox 1988 model was obtained with a 40% 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,DIAM,PRADO,PRADI,CD)
REAL(OM_FS,OMINJ,LAMVIS,ENUT_FS,TINJ,TIN_FS,REYNO)
BOOLEAN(KWMOD);KWMOD=F
REYNO=3.85E4
DIAM=0.058;PRADI=0.5*DIAM;PRADO=15.*PRADI
CD=0.1643
** jet-inflow conditions
WINJ=20.;TINJ=0.05
KEINJ=(TINJ*WINJ)**2; EPINJ=CD*KEINJ**1.5/(0.1*PRADI)
** laminar kinematic viscosity
LAMVIS=WINJ*DIAM/REYNO
** free-stream conditions
WIN_FS=WINJ/100.;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
CARTES=F;XULAST=0.1
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.*DIAM),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 Round Free Jet
+ MESG(Standard k-e model
+ TURMOD(KEMODL)
WHEN CK,2
+ TEXT(Chen-Kim KE_2D Elliptic Round Free Jet
+ MESG(Chen Kim k-e model
+ TURMOD(KECHEN)
WHEN RNG,3
+ TEXT(RNG KE_2D Elliptic Round Free Jet
+ MESG(RNG k-e model
+ TURMOD(KERNG)
WHEN RKE,3
+ TEXT(Realisable KE_2D Elliptic Round Free 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 Elliptic Round Free Jet
+ MESG(Wilcox 1988 k-w model
+ TURMOD(KOMODL);KWMOD=T
+ STORE(EP);OMINJ=EPINJ/(0.09*KEINJ)
+ OM_FS=EP_FS/(0.09*KE_FS)
** use enut_fs=0.1*enu_lam
+ OM_FS=KE_FS/(0.01*LAMVIS)
** lower C2f to reduce excessive jet spreading
+ SPEDAT(KECONST,C2E,R,0.6*0.075)
WHEN KWR,3
+ TEXT(Wilcox 2008 k-w Elliptic Round Free 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)
WHEN KWM,3
TEXT(Menter k-w 2D Elliptic Round Free Jet
+ MESG(Menter 1992 k-w model
+ TURMOD(KWMENTER);KWMOD=T
+ OMINJ=EPINJ/(0.09*KEINJ)
+ OM_FS=EP_FS/(0.09*KE_FS)
+ STORE(CDWS,SIGK,SIGW)
WHEN KWS,3
TEXT(SST k-w 2D Elliptic Round Free Jet
+ MESG(Menter 1992 k-w SST model
+ TURMOD(KWSST);KWMOD=T
+ OMINJ=EPINJ/(0.09*KEINJ)
+ OM_FS=EP_FS/(0.09*KE_FS)
+ STORE(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)
+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)
+LSWEEP=1800
+OUTPUT(EP,P,P,P,P,Y,Y)
WHEN KWM,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)
+LSWEEP=1800
+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,0.25)
+RELAX(KE,FALSDT,4.*RLXFAC);RELAX(OMEG,FALSDT,4.*RLXFAC)
+LSWEEP=1800
+OUTPUT(EP,P,P,P,P,Y,Y)
ENDCASE
GROUP 18. Limits on variables or increments to them
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
SPEDAT(SET,GXMONI,PLOTALL,L,T)
GROUP 23. Field print-out and plot control
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)
** check using ground coding in public ground.for
LG(1)=T
STORE(YHLF)
(make1 yhjet)
(store1 yhjet is (YHLF[,1,NZ-1]-YHLF[,1,:IZ1:])/(ZGM[,1,NZ-1]-ZGM[,1,:IZ1:]))
(print yhjet is yhjet)
** generate axial .csv file named IY1.csv
** compute jet half width & store at jet axis for profile plot
(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)
RESFAC=1.E-4
DISTIL=T
CASE :CTURB: OF
WHEN KE,2
+EX(P1 )=9.599E-02;EX(V1 )=1.382E-01
+EX(W1 )=3.588E+00;EX(KE )=2.341E+00
+EX(EP )=1.389E+02;EX(ZGM )=7.009E-01
+EX(YH )=1.195E-03;EX(YGP )=1.437E-01
+EX(WH )=6.250E+00;EX(EPKE)=1.994E+01
+EX(ENUT)=9.597E-03;EX(YH3D)=7.315E-02
+EX(LEN1)=1.389E-02
WHEN CK,2
+EX(P1 )=7.320E-02;EX(V1 )=1.227E-01
+EX(W1 )=3.723E+00;EX(KE )=2.074E+00
+EX(EP )=1.291E+02;EX(ZGM )=7.009E-01
+EX(YH )=9.887E-04;EX(YGP )=1.437E-01
+EX(WH )=6.880E+00;EX(EPKE)=2.261E+01
+EX(ENUT)=8.212E-03;EX(YH3D)=6.058E-02
+EX(LEN1)=1.448E-02
WHEN RNG,3
+EX(P1 )=1.099E-01;EX(V1 )=1.431E-01
+EX(W1 )=3.536E+00;EX(KE )=2.441E+00
+EX(EP )=1.281E+02;EX(ZGM )=7.009E-01
+EX(YH )=1.303E-03;EX(YGP )=1.437E-01
+EX(WH )=6.298E+00;EX(EPKE)=2.038E+01
+EX(ENUT)=1.253E-02;EX(YH3D)=8.003E-02
+EX(LEN1)=1.620E-02
WHEN RKE,3
+EX(P1 )=1.305E-01;EX(V1 )=1.469E-01
+EX(W1 )=3.445E+00;EX(KE )=2.526E+00
+EX(EP )=1.544E+02;EX(ZGM )=7.009E-01
+EX(YH )=1.415E-03;EX(YGP )=1.437E-01
+EX(WH )=5.621E+00;EX(C1E )=4.774E-01
+EX(DWDZ)=4.708E+00;EX(DWDY)=5.334E+01
+EX(DVDZ)=7.864E-01;EX(DVDY)=4.056E+00
+EX(DUDX)=2.478E+00;EX(EPKE)=1.887E+01
+EX(CMU )=1.103E-01;EX(ENUT)=1.300E-02
+EX(YH3D)=8.649E-02;EX(LEN1)=1.269E-02
WHEN KW,2
+EX(P1 )=9.930E-01;EX(V1 )=1.183E-01
+EX(W1 )=1.265E+00;EX(KE )=9.985E+00
+EX(EP )=5.880E+02;EX(YH3D)=2.467E-01
+EX(ZGM )=7.009E-01;EX(YH )=4.057E-03
+EX(YGP )=1.437E-01;EX(WH )=1.289E+00
+EX(LEN1)=3.344E-02;EX(OMEG)=1.116E+02
+EX(ENUT)=3.202E-02
WHEN KWR,3
+EX(P1 )=1.653E-01;EX(V1 )=1.550E-01
+EX(W1 )=3.359E+00;EX(KE )=3.095E+00
+EX(YH3D)=8.553E-02
+EX(ZGM )=7.009E-01;EX(LEN1)=2.541E-02
+EX(DWDZ)=4.919E+00;EX(DWDY)=4.797E+01
+EX(DVDZ)=1.216E+00;EX(DVDY)=4.177E+00
+EX(DUDX)=2.592E+00;EX(GEN1)=1.750E+04
+EX(FBP )=8.659E-01;EX(XWP )=4.811E+02
+EX(OMEG)=1.757E+02;EX(ENUT)=1.250E-02
+EX(LEN1)=1.945E-02;EX(WH )=5.212E+00
+EX(EP )=1.703E+02;EX(YH3D)=8.322E-02
+EX(YH )=1.363E-03;EX(YGP )=1.437E-01
WHEN KWM,3
+EX(P1 )=1.519E-01;EX(V1 )=1.564E-01
+EX(W1 )=3.351E+00;EX(KE )=2.869E+00
+EX(EP )=1.694E+02;EX(YH3D)=9.558E-02
+EX(ZGM )=7.009E-01;EX(YH )=1.518E-03
+EX(YGP )=1.437E-01;EX(WH )=5.219E+00
+EX(OMEG)=1.898E+02;EX(ENUT)=1.240E-02
+EX(LEN1)=1.651E-02;EX(YH3D)=9.256E-02
+EX(SIGW)=1.168E+00;EX(SIGK)=1.000E+00
+EX(CDWS)=5.716E+03
WHEN KWS,3
+EX(P1 )=1.462E-01;EX(V1 )=1.551E-01
+EX(W1 )=3.380E+00;EX(KE )=2.804E+00
+EX(EP )=1.654E+02;EX(YH3D)=9.369E-02
+EX(ZGM )=7.009E-01;EX(YH )=1.490E-03
+EX(YGP )=1.437E-01;EX(WH )=5.331E+00
+EX(GEN1)=1.638E+04;EX(CDWS)=5.849E+03
+EX(OMEG)=1.928E+02;EX(ENUT)=1.197E-02
+EX(LEN1)=1.605E-02;EX(YH3D)=9.113E-02
+EX(SIGW)=1.168E+00;EX(SIGK)=1.000E+00
ENDCASE
LIBREF = 309
STOP