****** TO LOAD CASE 113:TYPE L(N113) *****
    GROUP 1. Run title and other preliminaries
TEXT(2D TRANSONIC UNDEREXPANDED JET: N113
TITLE
mesg(PC486/50 time last reported as appx. 2 hrs
  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 is represented by means of the k-e
  turbulence model. The calculation may be made with the van-Leer
  higher-order scheme or the hybrid differencing scheme. A
  relatively coarse mesh of 30 radial by 90 axial cells is used in
  the computations.
  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
     PHI 5
 
     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)
CHAR(SCHM);BOOLEAN(HSTAG)
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,VABS);TURMOD(KEMODL)
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)
    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
MESG( Enter required convection scheme
MESG(  Default: VANL1 for all variables
MESG( The alternative is:
MESG(  HYB - Hybrid differencing for all variables
READVDU(SCHM,CHAR,HOC)
CASE :SCHM: OF
WHEN HYB,3
+ MESG(Hybrid-differencing scheme
+ DTF=ZWLAST/(WIN*NZ);DIFCUT=0.5;LSWEEP=850
WHEN HOC,3
+ MESG(VANL1 for all variables
+ SCHEME(VANL1,ALL)
+ DTF=0.25*ZWLAST/(NZ*WIN);LSWEEP=1600
ENDCASE
    GROUP 9. Properties of the medium (or media)
RHO1=IDEALGAS;RHO1B=1./RCON;PRESS0=0.0;DRH1DP=IDEALGAS;RHO1C=RGAM
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)
FIINIT(EP)=EPSIN; 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
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);COVAL(IN,EP,ONLYMS,EPSIN)
 
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
COVAL(LB,KE,ONLYMS,1.E-10)
COVAL(LB,EP,ONLYMS,.09*1.E-20/(0.01*ENUL))
 
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,1.E-10)
COVAL(NB,EP,ONLYMS,0.09*1.E-20/(0.01*ENUL))
 
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);COVAL(OUT,EP,ONLYMS,SAME)
    GROUP 15. Termination of sweeps
    GROUP 16. Termination of iterations
LITER(P1)=10;DENPCO=T
    GROUP 17. Under-relaxation devices
KELIN=1;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);RELAX(EP,FALSDT,DTF)
    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