GROUP 1. Run title and other preliminaries
TEXT(2D Supersonic Rocket Exhaust Plume
TITLE
  DISPLAY
  The problem considered is an axisymmetric, supersonic under-
  expanded rocket exhaust plume discharging into a subsonic moving
  air stream. The case is the same as library 911, except that
  here parabolic marching solution is obtained in a single sweep
  using the IPARAB=5 option, rather than an elliptic multi-sweep
  solution. 
  
  The exhaust gas has a specific-heat ratio of 1.21, a Mach number 
  of 2.39, a total pressure of 5.786MPa, and a total temperature of 
  3462 K.
  
  The free-stream has a Mach number of 0.029, a static temperature
  of 288K, and a static pressure of 1.01325 bar. The nozzle-to-
  ambient static pressure ratio is 3.8.
  ENDDIS
   PHOTON USE
   P
   PARPHI
 
 
   VEC X 1 SH
   pau;cl
   con mach x 1 fi;.1
   pau;cl
   con p1 x 1 fi;.1
   pau cl;
   con tmp1 x 1 fi;.1
  ENDUSE
    
   AUTOPLOT USE
     FILE
     PARPHI 5
 
     D 1 MACH Y 1;PLOT;REDR;LEVEL Y 1
     msg Mach number distribution along flow axis
     msg Press  to continue
     pause
     cl
     D 1 P1 Y 1;PLOT
     shift y
     1.01325e5 1
     scale;redr
     msg static pressure distribution along flow axis
     msg Press  to continue
     pause
     cl
     D 1 RHO1 Y 1;PLOT;REDR
     msg static density distribution along flow axis
     msg Press  to continue
     pause
     cl
 
   ENDUSE
 
REAL(AIN,FLOWIN,HTOTJ,PA,HA,HTOTA,MAJ,PJ,HJ,RHOJ,RHOA,TKEA,EPSA)
REAL(WJ,TJ,EPSJ,TKEJ,DTF,DN,RAD,PRAT,GASCJ,RGAM,TTOT,FLOWJ)
REAL(MAA,WA,TA,GASCU,GASCA,MWTJ,MWTA,GAMJ,GAMA,CPJ,CPA,AVJ)
   **  Gas properties
GASCU=8314.43
   **  Jet inlet conditions
RAD=0.036325;DN=2.*RAD
WJ=2246.34; TJ=2163.0; PJ=3.851E5; MWTJ=24.66; GAMJ=1.21  ! jet-gas properties
GASCJ=GASCU/MWTJ; CPJ=GASCJ*GAMJ/(GAMJ-1.)
RHOJ=PJ/(GASCJ*TJ)
HJ=CPJ*TJ;HTOTJ=HJ+0.5*WJ*WJ;AVJ=(GAMJ*GASCJ*TJ)**0.5
MAJ=WJ/AVJ
    ** Ambient-air conditions
PA=1.01325E5;WA=10.0;TA=288.0
MWTA=28.76; GAMA=1.4                                      ! air properties
GASCA=GASCU/MWTA;CPA=GASCA*GAMA/(GAMA-1.)
RHOA=PA/(GASCA*TA)
HA=CPA*TA;HTOTA=HA+0.5*WA*WA
MAA=WA/(GAMA*GASCA*TA)**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
NREGY=2
IREGY=1;GRDPWR(Y,20,RAD,1.0);IREGY=2;GRDPWR(Y,20,2.*DN,1.4)
    GROUP 5. Z-direction grid specification
NZ=500;AZDZ=PROPY;DZW1=0.01
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1,H1,C1)
STORE(RHO1,MACH,ENUT,LEN1,TMP1);TURMOD(KEMODL)
SOLUTN(V1,P,P,Y,Y,P,N);SOLUTN(W1,P,P,Y,Y,P,N)
SOLUTN(KE,P,P,P,P,P,N);SOLUTN(EP,P,P,P,P,P,N)
IPARAB=5;STORE(MACZ)
   ** Subsonic-suppression device
RMACHZ=100.0
    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,P,P,P,P,P)
DIFCUT=0.0
    GROUP 9. Properties of the medium (or media)
RHO1=IDEALGAS;RHO1B=1./GASCA;PRESS0=PA    ! IDEALGAS = GRND5
RHO1A=GASCA; RHO1B=GASCJ                  ! gas constants of air and jet
DRH1DP=IDEALGAS; RHO1C=1./GAMA
TMP1=VARSTAGH;                            ! VARSTAGH = GRND6; so temperature
                                          ! is computed from the computed
                                          ! stagnation enthalpy
CP1=GRND7; CP1A=CPA; CP1B=CPJ             ! CP in linear in C1
ENUL=1.E-5                                  
    GROUP 11. Initialization of variable or porosity fields
TKEJ=(0.1*WJ)**2; EPSJ=0.1643*TKEJ**1.5/(0.1*RAD)
FIINIT(EP)=EPSJ; FIINIT(KE)=TKEJ;FIINIT(W1)=WJ;FIINIT(V1)=0.0
FIINIT(RHO1)=RHOJ;FIINIT(P1)=PJ-PA;FIINIT(H1)=HTOTJ
FIINIT(TMP1)=TJ;FIINIT(C1)=1.0;FIINIT(MACZ)=MAJ
INIADD=F
PATCH(INITA,INIVAL,1,1,#2,#2,1,NZ,1,1)
COVAL(INITA,P1,zero,0.0);COVAL(INITA,W1,ZERO,WA)
COVAL(INITA,C1,ZERO,0.);COVAL(INITA,H1,ZERO,HTOTA)
COVAL(INITA,RHO1,ZERO,RHOA);COVAL(INITA,TMP1,ZERO,TA)
TKEA=1.E-10;EPSA=.09*TKEA*TKEA/(0.01*ENUL)
COVAL(INITA,KE,ZERO,TKEA);COVAL(INITA,EP,ZERO,EPSA)
COVAL(INITA,MACZ,ZERO,MAA)
    GROUP 13. Boundary conditions and special sources
PATCH(IN,LOW,1,1,#1,#1,1,1,1,1)
COVAL(IN,P1,FIXFLU,RHOJ*WJ)
COVAL(IN,W1,ONLYMS,WJ);COVAL(IN,V1,ONLYMS,0.0)
COVAL(IN,H1,ONLYMS,HTOTJ);COVAL(IN,C1,ONLYMS,1.0)
COVAL(IN,KE,ONLYMS,TKEJ);COVAL(IN,EP,ONLYMS,EPSJ)
 
PATCH(LB,LOW,1,1,#2,#2,1,1,1,1)
COVAL(LB,P1,FIXFLU,RHOA*WA);COVAL(LB,W1,ONLYMS,WA)
COVAL(LB,H1,ONLYMS,HTOTA)
COVAL(LB,KE,ONLYMS,TKEA);COVAL(LB,EP,ONLYMS,EPSA)
 
PATCH(NB,NORTH,1,1,NY,NY,1,NZ,1,1)
COVAL(NB,P1,1.E3,0.)
COVAL(NB,W1,ONLYMS,WA);COVAL(NB,C1,ONLYMS,0.0)
COVAL(NB,H1,ONLYMS,HTOTA)
COVAL(NB,KE,ONLYMS,TKEA);COVAL(NB,EP,ONLYMS,EPSA)
    GROUP 15. Termination of sweeps
LITHYD=10
    GROUP 16. Termination of iterations
DENPCO=T
FLOWJ=RHOJ*WJ*AIN;FLOWIN=FLOWJ
SELREF=F;RESREF(P1)=1.E-12*FLOWIN
RESREF(W1)=1.E-12*FLOWIN*WJ;RESREF(C1)=RESREF(P1)
RESREF(KE)=1.E-12*FLOWIN*TKEJ;RESREF(EP)=1.E-12*FLOWIN*EPSJ
RESREF(V1)=RESREF(W1);RESREF(U1)=RESREF(W1)
RESREF(H1)=1.E-12*FLOWIN*HTOTJ
    GROUP 17. Under-relaxation devices
DTF=2.*ZWLAST/(WJ+AVJ)/NZ
KELIN=3;RELAX(P1,LINRLX,0.5)
RELAX(W1,FALSDT,DTF);RELAX(V1,FALSDT,DTF)
RELAX(H1,FALSDT,DTF);RELAX(C1,LINRLX,1.0)
RELAX(KE,FALSDT,DTF);RELAX(EP,FALSDT,DTF)
RELAX(RHO1,LINRLX,1.0);RELAX(MACZ,LINRLX,0.8)
    GROUP 18. Limits on variables or increments to them
VARMIN(RHO1)=1.E-5;VARMIN(TMP1)=200.
VARMIN(C1)=1.E-10;VARMAX(C1)=1.0
    GROUP 21. Print-out of variables
OUTPUT(MACH,P,P,P,P,Y,P);OUTPUT(RHO1,P,P,P,P,Y,P)
    GROUP 22. Spot-value print-out
TSTSWP=-1;IYMON=NY/2
    GROUP 23. Field print-out and plot control
ITABL=3;NPLT=1;NYPRIN=1;NZPRIN=NZ/5
IF(NZ.GT.1) THEN
+ IDISPA=nz/20;IDISPB=1;IDISPC=NZ
ENDIF
    GROUP 24. Preparations for continuation runs