TEXT(2D Sonic Underexpanded Round Jet
TITLE
  DISPLAY
  The problem considered is the near-field of a sonic under-
  expanded turbulent round jet issuing into stagnant surroundings.
  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 the k-e model, and the IPARAB=5 option of the
  parabolic solver is used to perform the marching integration.
  The calculations are started at the jet discharge plane, and
  carried out until 6 diameters downstream. The calculations are
  made with 35 radial grid cells and a forward step size of 1%
  (DZW1) of the local grid width. The main input parameters are
  the inlet and free-stream Mach numbers and the jet-to-ambient
  static pressure ratio. This case has been tested for inlet
  Mach numbers up to 2 and pressure ratios up to 2.5.
  ENDDIS
   So as to allow a direct computation of dimensionless flow
   variables, the flow equations are normalised such that the flow
   variables can be interpreted as: P/Po; RHO/RHOo; T/To; U/Uref;
   H/Href; KE/Uref**2; EP*L/Uref**3; and ENUT/(Uref*L).
 
   Here: Po, RHOo and To are the inlet total pressure, density
   and temperature; Uref=Ao/SQRT(GAMMA); L is the nozzle diameter;
   and Href=(gam-1.)*Ho/gam, where Ho is the inlet total enthalpy
   and gam is the specific heat ratio. Ao is the acoustic velocity
   at To (see Palacio et al, Int.J.Heat Mass Transfer, Vol.33,
   No.6, p1193, [1990] ).
 
 
  PHOTON USE
   P
   PARPHI
 
 
   VEC X 1 SH
   pau;cl
   con mach x 1 fi;.5
   pau;cl
   con p1 x 1 fi;.5
   pau cl;
   con tmp1 x 1 fi;.5
  ENDUSE
 
REAL(AIN,CP,GAM,GM1,PTOT,HTOT,RHTOT,PAMB,HAMB,MIN,PIN,HIN,RHOIN)
REAL(WIN,EPSIN,TKEIN,DTF,DN,RAD,PRAT,RCON,RGAM,TTOT,TIN,UAC,RHAMB)
REAL(FLOWIN,AOIN,UREF,DYGDZ,TAMB,MAMB,WAMB)
CHAR(CTURB);BOOLEAN(GEXPAN)
   **  GEXPAN=T activates a linear y-grid expansion with z.
       A linear expansion is inappropriate, but it suffices to
       test for convergence.
GEXPAN=F
   **  Gas properties
GAM=1.4;GM1=GAM-1.;RCON=1.0;CP=RCON*GAM/GM1;RGAM=1./GAM
 
   **  INLET CONDITIONS & NOZZLE DIAMETER
PTOT=1.0;RHTOT=1.0;TTOT=1.0;MIN=1.0;DN=1.0;HTOT=CP*TTOT
 
RAD=0.5*DN
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
AOIN=(GAM*PTOT/RHTOT)*0.5;UREF=AOIN/GAM**0.5
    ** Set static pressure ratio, i.e. PIN/PAMB
PRAT=1.85
    ** Free-stream Mach number
MAMB=1.E-3
    ** Ambient conditions
PAMB=PIN/PRAT;HAMB=HTOT
TAMB=HAMB/CP;RHAMB=PAMB/(RCON*TAMB)
WAMB=MAMB*(GAM*PAMB/RHAMB)**0.5
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
PARAB=T
    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;NREGY=2
IF(GEXPAN) THEN
+ DYGDZ=0.0875;AZYV=1.0;ZWADD=DN/DYGDZ
  ** The following line activates parabolic enhancement
     which accounts for grid-angle inclination in the
     calculation of the radial convection fluxes.
+ V1AD=GRND1
ENDIF
IREGY=1;GRDPWR(Y,20,RAD,1.0);IREGY=2;GRDPWR(Y,15,2.*RAD,2.0)
    GROUP 5. Z-direction grid specification
NZ=400;AZDZ=PROPY;DZW1=0.01
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1)
  ** Point-by-point velocity solution for IPARAB=5
SOLUTN(W1,P,P,Y,Y,P,P);SOLUTN(V1,P,P,Y,Y,P,P)
  ** Provide storage for the density.
STORE(RHO1,MACH,ENUT,LEN1,H1,TMP1)
  ** Activate option to handle supersonic flow with a
     subsonic free stream by the following line of commands.
IPARAB=5;STORE(MACZ);RMACHZ=1.0
  ** For MIN>1 set RMACHZ=100.0 for convergence
IF(MIN.GT.1) THEN
+ RMACHZ=100.
ELSE
+ RMACHZ=1.0
ENDIF
 
MESG( Enter the required turbulence model:
MESG(  KE   -  Standard k-e model (Default)
MESG(  LAM  -  Laminar model
MESG(
READVDU(CTURB,CHAR,KE)
CASE :CTURB: OF
WHEN KE,2
+ MESG(Standard k-e model
+ TURMOD(KEMODL);KELIN=3
WHEN LAM,3
+ MESG(Laminar model
+ ENUT=0.
+ RMACHZ=1.E3
ENDCASE
 
    GROUP 8. Terms (in differential equations) & devices
DIFCUT=0;DENPCO=T
  ** Deactivate radial diffusion of V1 for IPARAB=5
TERMS(V1,P,P,N,P,P,P)
    GROUP 9. Properties of the medium (or media)
RHO1=IDEALGAS;RHO1B=1./RCON;PRESS0=0.0
DRH1DP=IDEALGAS;RHO1C=RGAM
TMP1=VARSTAGH;CP1=CP
ENUL=1.8E-5/(UREF*DN)
    GROUP 10. Inter-phase-transfer processes and properties
    GROUP 11. Initialization of variable or porosity fields
TKEIN=(0.05*WIN)**2;EPSIN=0.1643*TKEIN**1.5/(0.1*RAD/DN)
FIINIT(EP)=EPSIN; FIINIT(KE)=TKEIN
FIINIT(W1)=WIN;FIINIT(RHO1)=RHOIN;FIINIT(P1)=PIN
FIINIT(H1)=HTOT;FIINIT(MACZ)=MIN
FIINIT(TMP1)=TIN
INIADD=F
PATCH(INITFS,INIVAL,1,1,#2,#2,1,1,1,1)
COVAL(INITFS,W1,0.0,0.0);COVAL(INITFS,P1,0.0,PAMB)
COVAL(INITFS,TMP1,0.0,TAMB);COVAL(INITFS,RHO1,0.0,RHAMB)
COVAL(INITFS,MACZ,0.0,MAMB)
    GROUP 13. Boundary conditions and special sources
INLET(IN,LOW,1,NX,#1,#1,1,1,1,1)
VALUE(IN,P1,RHOIN*WIN);VALUE(IN,W1,WIN)
VALUE(IN,KE,TKEIN);VALUE(IN,EP,EPSIN)
 
INLET(INFS,LOW,1,NX,#2,#2,1,1,1,1)
VALUE(INFS,P1,RHAMB*WAMB);VALUE(INFS,W1,WAMB)
 
PATCH(NB,NORTH,1,1,NY,NY,1,NZ,1,1)
COVAL(NB,P1,1.E3,PAMB);COVAL(NB,W1,ONLYMS,WAMB)
 
    GROUP 16. Termination of iterations
LITER(P1)=100
FLOWIN=RHOIN*WIN*AIN
SELREF=F;RESREF(P1)=1.E-12*FLOWIN
RESREF(W1)=1.E-12*FLOWIN*WIN
RESREF(KE)=1.E-12*FLOWIN*TKEIN;RESREF(EP)=1.E-12*FLOWIN*EPSIN
RESREF(V1)=RESREF(W1)
    GROUP 17. Under-relaxation devices
DTF=100.*DZW1*YVLAST/(WIN+UAC)
RELAX(V1,LINRLX,0.5);RELAX(W1,LINRLX,0.5)
RELAX(KE,LINRLX,0.3);RELAX(EP,LINRLX,0.3)
RELAX(P1,LINRLX,0.7);RELAX(RHO1,LINRLX,1.0)
RELAX(MACZ,LINRLX,0.8)
    GROUP 18. Limits on variables or increments to them
VARMIN(W1)=1.E-10;VARMIN(P1)=1.E-4*PTOT
VARMIN(RHO1)=0.001*RHTOT;VARMIN(H1)=1.0E-10
VARMAX(RHO1)=5.0*RHTOT;VARMAX(H1)=HTOT;VARMAX(P1)=1.0E10*PTOT
    GROUP 19. Data communicated by satellite to GROUND
    GROUP 20. Preliminary print-out
OUTPUT(MACH,P,P,P,P,Y,P)
    GROUP 21. Print-out of variables
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=1;IZMON=1
    GROUP 23. Field print-out and plot control
ITABL=2;NPLT=2;NYPRIN=2;NZPRIN=NZ/5
IF(NZ.GT.1) THEN
+ IDISPA=2;IDISPB=1;IDISPC=NZ
ENDIF
LITHYD=20;TSTSWP=-1;IPLTL=LITHYD