TALK=T;RUN( 1, 1)
GROUP 1. Run title and other preliminaries
TEXT(2D TRANSONIC UNDEREXPANDED FREE JET: T107
TITLE
DISPLAY
The problem considered is an axisymmetric sonic jet discharging
into stagnant surroundings from a nozzle at a pressure 3.56
times higher than the ambient pressure. The stagnation enthalpy of
the nozzle fluid is equal to that of the free stream, so that with
the assumption of unit Prandtl numbers, the energy equation need
not be solved. The turbulence can be represented by means of either
the k-e model, the Wilcox 2008 k-w model, or the k-w-SST model. A
relatively coarse mesh of 30 radial by 90 axial cells is used in
the computations. This case is the same as library N113 case, which
uses the standard k-e model togther with the higher-order van-Leer
convection scheme. Here, the default hybrid differencing scheme is
used in the simulations.
ENDDIS
This case was studied experimentally by Donaldson and Snedeker
(J.Fluid Mech, 45, p281, [1971]) and numerically by Palacio et al
(Int.J.Heat Mass Transfer, 33, p1193, [1990] ). The flow shows a
rapid initial expansion of the nozzle fluid, with the experiments
indicating Mach-disc formation at z/D=1.54 with a Mach number of
3.5 just upstream of the disc and a Mach number of 0.5 just
downstream of the disc.
PHOTON USE
P
0.20443E+04 0.15633E+04 CR
CON MACH X 1 FI;.5
msg Mach number contours
msg Press to continue
pause
ENDUSE
AUTOPLOT USE
file
phida 3
D 1 MACH Y 1;PLOT;REDR;LEVEL Y 1;level x 1.54
msg Mach number distribution along flow axis
msg Press to continue
pause
ENDUSE
REAL(AIN,CP,GAM,GM1,PTOT,HTOT,RHTOT,PEXIT,HAMB,MIN,PIN,HIN,RHOIN)
REAL(WIN,EPSIN,TKEIN,DTF,DN,RAD,PRAT,RCON,RGAM,TTOT,TIN,TKEFRE)
REAL(EPSFRE)
CHAR(CTURB);BOOLEAN(HSTAG,KWMOD)
HSTAG=T
** Gas properties
GAM=1.4;GM1=GAM-1.;RCON=1.0;CP=RCON*GAM/GM1;RGAM=1./GAM
** Inlet conditions
PTOT=1.0;RHTOT=1.0;TTOT=1.0;MIN=1.0;DN=1.0;RAD=0.5*DN
HTOT=CP*TTOT
PIN=PTOT/(1.+GM1*MIN*MIN/2.)**(GAM/GM1)
RHOIN=RHTOT*(PIN/PTOT)**RGAM;WIN=MIN*(GAM*PIN/RHOIN)**0.5
TIN=PIN/(RCON*RHOIN);HIN=CP*TIN
** Exit pressure
PRAT=3.563;PEXIT=PIN/PRAT;HAMB=HTOT
GROUP 3. X-direction grid specification
CARTES=F;XULAST=0.1;AIN=0.5*XULAST*RAD*RAD
GROUP 4. Y-direction grid specification
NY=30;YVLAST=DN;NREGY=2
IREGY=1;GRDPWR(Y,20,RAD,1.0);IREGY=2;GRDPWR(Y,10,RAD,1.3)
GROUP 5. Z-direction grid specification
NREGZ=3
IREGZ=1;GRDPWR(Z,65,2.5,1.0);IREGZ=2;GRDPWR(Z,15,1.5,1.2)
IREGZ=3;GRDPWR(Z,10,1.0,1.2)
GROUP 7. Variables stored, solved & named
SOLVE(P1);SOLUTN(P1,Y,Y,Y,N,N,N)
SOLVE(V1,W1);STORE(RHO1,MACH,ENUT,LEN1,VABS)
IF(HSTAG) THEN
+ STORE(H1);STORE(TMP1)
ELSE
+ STORE(HTOT);SOLVE(H1);SOLUTN(H1,Y,Y,Y,N,N,N)
ENDIF
SOLUTN(V1,P,P,P,P,P,N);SOLUTN(W1,P,P,P,P,P,N)
SOLUTN(KE,P,P,P,P,P,N);SOLUTN(EP,P,P,P,P,P,N)
MESG( Enter the required turbulence model:
MESG( KE - Standard k-e model (default)
MESG( KWR - Wilcox 2008 k-w model
MESG( KWS - k-w SST model
MESG(
READVDU(CTURB,CHAR,KE)
CASE :CTURB: OF
WHEN KE,2
TEXT(T107:KE-2D TRANSONIC UNDEREXPANDED JET
+ MESG(Standard k-e model (default)
+ TURMOD(KEMODL);KWMOD=F
+ SOLUTN(KE,P,P,P,P,P,N);SOLUTN(EP,P,P,P,P,P,N)
WHEN KWR,3
+ TEXT(T107:KW 2008-2D UNDEREXPANDED JET
+ MESG(Wilcox 2008 k-w model
+ TURMOD(KWMODLR);KWMOD=T;STORE(FBP);FIINIT(FBP)=1.0
WHEN KWS,3
TEXT(T107:KW SST-2D UNDEREXPANDED JET
+ MESG( k-w SST model
+ TURMOD(KWSST);KWMOD=T;STORE(BF1,BF2,GEN1)
+ STORE(SIGK,SIGW,CDWS)
+ FIINIT(BF1)=0.0;FIINIT(BF2)=0.0
ENDCASE
GROUP 8. Terms (in differential equations) & devices
IF (.NOT.HSTAG) THEN
+ TERMS(H1,Y,Y,Y,N,Y,N)
ENDIF
** Activate compressibility corrections of Malin & Sanchez
UCONV=T;NAMGRD=CONV
GROUP 9. Properties of the medium (or media)
RHO1=IDEALGAS;RHO1B=1./RCON;PRESS0=0.0;DRH1DP=IDEALGAS;RHO1C=RGAM
** Enthalpy-temperature relationship
IF(HSTAG) THEN
+ TMP1=VARSTAGH;CP1=CP
ELSE
+ TMP1=LINH;TMP1A=0.0;CP1=CP
ENDIF
ENUL=2.E-5
GROUP 11. Initialization of variable or porosity fields
TKEIN=(0.05*WIN)**2;EPSIN=0.1643*TKEIN**1.5/(0.1*RAD)
IF(KWMOD) THEN
+ EPSIN=EPSIN/(0.09*TKEIN)
+ FIINIT(OMEG)=EPSIN
ELSE
+ FIINIT(EP)=EPSIN
ENDIF
FIINIT(KE)=TKEIN;FIINIT(W1)=WIN;FIINIT(V1)=0.0
FIINIT(RHO1)=0.5*RHOIN;FIINIT(P1)=0.5*(PIN+PEXIT)
IF(HSTAG) THEN
+ FIINIT(H1)=3.5
ELSE
+ FIINIT(H1)=0.5*HIN
ENDIF
GROUP 13. Boundary conditions and special sources
** Nozzle discharge
PATCH(IN,LOW,1,1,1,NY/2,1,1,1,1)
COVAL(IN,P1,FIXFLU,RHOIN*WIN)
COVAL(IN,W1,ONLYMS,WIN);COVAL(IN,V1,ONLYMS,0.0)
IF (.NOT.HSTAG) THEN
+ COVAL(IN,H1,ONLYMS,HIN)
ENDIF
COVAL(IN,KE,ONLYMS,TKEIN)
IF(KWMOD) THEN
+COVAL(IN,OMEG,ONLYMS,EPSIN)
ELSE
+COVAL(IN,EP,ONLYMS,EPSIN)
ENDIF
** low free boundary
PATCH(LB,LOW,1,1,NY/2+1,NY,1,1,1,1)
COVAL(LB,P1,1.0E3,PEXIT)
IF (.NOT.HSTAG) THEN
+ COVAL(LB,H1,ONLYMS,HTOT)
ENDIF
TKEFRE=1.E-10
EPSFRE=0.09*TKEFRE*TKEFRE/(0.01*ENUL)
COVAL(LB,KE,ONLYMS,TKEFRE)
IF(KWMOD) THEN
+EPSFRE=TKEFRE/(0.01*ENUL)
+COVAL(LB,OMEG,ONLYMS,EPSFRE)
ELSE
+COVAL(LB,EP,ONLYMS,EPSFRE)
ENDIF
** north free boundary
PATCH(NB,NORTH,1,1,NY,NY,1,NZ,1,1)
COVAL(NB,P1,1.E3,PEXIT)
IF (.NOT.HSTAG) THEN
+ COVAL(NB,H1,ONLYMS,HTOT)
ENDIF
COVAL(NB,KE,ONLYMS,TKEFRE)
IF(KWMOD) THEN
+EPSFRE=TKEFRE/(0.01*ENUL)
+COVAL(LB,OMEG,ONLYMS,EPSFRE)
ELSE
+COVAL(NB,EP,ONLYMS,EPSFRE)
ENDIF
** High outflow boundary
PATCH(OUT,HIGH,1,1,1,NY,NZ,NZ,1,1)
COVAL(OUT,P1,1.0E3,PEXIT)
IF (.NOT.HSTAG) THEN
+ COVAL(OUT,H1,ONLYMS,SAME)
ENDIF
COVAL(OUT,KE,ONLYMS,SAME)
IF(KWMOD) THEN
+COVAL(OUT,OMEG,ONLYMS,SAME)
ELSE
+COVAL(OUT,EP,ONLYMS,SAME)
ENDIF
GROUP 15. Termination of sweeps
GROUP 16. Termination of iterations
LITER(P1)=20;DENPCO=T
GROUP 17. Under-relaxation devices
DTF=ZWLAST/(WIN*NZ);LSWEEP=1200
RELAX(P1,LINRLX,0.7)
RELAX(W1,FALSDT,DTF);RELAX(V1,FALSDT,DTF)
IF (.NOT.HSTAG) THEN
+ RELAX(H1,FALSDT,DTF)
ENDIF
RELAX(KE,FALSDT,DTF)
IF(KWMOD) THEN
+RELAX(OMEG,FALSDT,DTF)
ELSE
+RELAX(EP,FALSDT,DTF);KELIN=1
ENDIF
GROUP 18. Limits on variables or increments to them
VARMIN(V1)=-100.;VARMIN(W1)=-100.;VARMAX(V1)=100.;VARMAX(W1)=100.
VARMIN(RHO1)=1.0E-10*RHTOT;VARMIN(H1)=1.0E-10;VARMIN(P1)=1.0E-10
VARMAX(RHO1)=1.0E10*RHTOT;VARMAX(H1)=HTOT;VARMAX(P1)=1.0E10*PTOT
GROUP 21. Print-out of variables
OUTPUT(MACH,P,P,P,P,Y,P);OUTPUT(RHO1,P,P,P,P,Y,P)
IF (.NOT.HSTAG) THEN
+ OUTPUT(H1,Y,N,N,Y,Y,Y)
ENDIF
GROUP 22. Spot-value print-out
TSTSWP=-1
IF(DIFCUT.EQ.0.0) THEN
+ IZMON=78;IZPRF=43;IZPRL=55
ELSE
+ IZMON=78;IZPRF=53;IZPRL=60
ENDIF
NZPRIN=1;NYPRIN=1;IYPRF=1;IYPRL=5
GROUP 23. Field print-out and plot control
IF (.NOT.HSTAG) THEN
+ PATCH(PLOT3,PROFIL,1,1,1,1,1,NZ,1,1)
+ PLOT(PLOT3,H1,0.0,0.0)
ENDIF
PATCH(PLOT4,PROFIL,1,1,1,1,1,NZ,1,1);PLOT(PLOT4,W1,0.0,0.0)
PATCH(PLOT5,PROFIL,1,1,1,1,1,NZ,1,1);PLOT(PLOT5,MACH,0.0,0.0)
NPRINT=LSWEEP;NPLT=20;ITABL=3
GROUP 24. Preparations for continuation runs
SPEDAT(SET,GXMONI,PLOTALL,L,T)
SPEDAT(SET,OUTPUT,NOFIELD,L,T)
OUTPUT(ENUT,Y,N,Y,N,Y,Y)
DISTIL=T
CASE :CTURB: OF
WHEN KE,2
+EX(P1 )=1.490E-01;EX(V1 )=1.052E-01
+EX(W1 )=1.175E+00;EX(KE )=4.568E-02
+EX(EP )=3.175E-02;EX(H1 )=3.500E+00
+EX(TMP1)=7.579E-01;EX(VABS)=1.181E+00
+EX(LEN1)=5.947E-02;EX(ENUT)=6.952E-03
+EX(MACH)=1.226E+00;EX(RHO1)=1.960E-01
WHEN KWR,3
+EX(P1 )=1.486E-01;EX(V1 )=9.630E-02
+EX(W1 )=1.146E+00;EX(KE )=3.516E-02
+EX(EP )=1.940E-02;EX(H1 )=3.500E+00
+EX(DWDZ)=5.839E-01;EX(DWDY)=1.833E+00
+EX(DVDZ)=3.029E-01;EX(DVDY)=6.689E-01
+EX(DUDX)=4.111E-01;EX(GEN1)=1.559E+01
+EX(FBP )=8.682E-01;EX(XWP )=2.212E+04
+EX(OMEG)=4.393E+00;EX(TMP1)=7.679E-01
+EX(VABS)=1.152E+00;EX(LEN1)=8.968E-01
+EX(ENUT)=5.404E-03;EX(MACH)=1.191E+00
+EX(RHO1)=1.935E-01
WHEN KWS,3
+EX(P1 )=1.489E-01;EX(V1 )=1.032E-01
+EX(W1 )=1.165E+00;EX(KE )=3.353E-02
+EX(EP )=2.418E-02;EX(H1 )=3.500E+00
+EX(CDWS)=9.021E-01;EX(SIGW)=1.168E+00
+EX(SIGK)=1.000E+00;EX(GEN1)=1.352E+01
+EX(OMEG)=6.586E+00;EX(TMP1)=7.610E-01
+EX(VABS)=1.171E+00;EX(LEN1)=4.918E-02
+EX(ENUT)=4.971E-03;EX(MACH)=1.216E+00
+EX(RHO1)=1.952E-01
ENDCASE
STOP