TEXT(2D Flow In A Supersonic Diffuser
TITLE
  DISPLAY
  The problem considered is 2D supersonic flow in a 'supersonic
  diffuser', i.e. a converging duct. The problem is solved using
  the IPARAB=4 option of the parabolic solver. The flow enters
  axially at MACH 2, and then passes over a 8o corner where a weak
  oblique shock wave is generated of angle 37o, which is in turn
  reflected from the bottom of the duct to arrive at the top corner
  of the exit plane. The PHOENICS results are compared with the
  results of shock theory below:
 
                M1    M2    M3   |  P1/Po1   P2/Po1  P3/Po1
    Theory     2.0   1.71  1.42  |  0.128    0.197   0.293
    PHOENICS   2.0   1.75  1.45  |  0.128    0.192   0.293
 
  Here, M, P and Po denote Mach number, static pressure and total
  pressure, respectively; and 1, 2 and 3 denote inlet, post
  shock, and post reflected-shock, respectively.
  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;
   and U/Uref.
 
   Here: Po, RHOo and To are the inlet total pressure, density
   and temperature; Uref=Ao/SQRT(gam); 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
 
 
 
    con mach x 1 fi;.01
    pau;cl
    vec x 1 sh
    pau;cl
    con p1 x 1 fi;.01
    pau;cl
    con rho1 x 1 fi;.01
   ENDUSE
REAL(GASCON,GAMMA,PTOTAL,TTOTAL,RHOTOT,MACHI,PEXRAT,AGAM1,RGAM)
REAL(DTF,PIN,TIN,POWER,WIN,RHOIN,PEXIT,CHORD)
REAL(AIN,FLOWIN,ANG1,PI,TANA,YIN,YOUT,ZLEN)
PI=3.1415927
 
GASCON=1.0;GAMMA=1.4;PTOTAL=1.0;TTOTAL=1.0;RHOTOT=1.0
 
MACHI=2.0;YIN=1.0;ZLEN=2.3
   ** Corner angle
ANG1=8.0;ANG1=ANG1*PI/180.;TANA=-TAN(ANG1)
 
   ** Calculation of inlet velocity
AGAM1=GAMMA-1.;RGAM=1./GAMMA;POWER=GAMMA/AGAM1
PIN=PTOTAL/(1.+AGAM1*MACHI*MACHI/2.)**POWER
RHOIN=RHOTOT/(PTOTAL/PIN)**RGAM
WIN=MACHI*(GAMMA*PIN/RHOIN)**0.5
   ** Calculation of Inlet Temperature
TIN=PIN/(GASCON*RHOIN)
    GROUP 1. Run title and other preliminaries
    GROUP 2. Transience; time-step specification
   ** activate wholly-supersonic 'parabolic' solver
PARAB=T;IPARAB=4
    GROUP 4. Y-direction grid specification
GRDPWR(Y,80,YIN,1.0)
AZYV=1.0;ZWADD=YIN/TANA
    GROUP 5. Z-direction grid specification
GRDPWR(Z,320,ZLEN,1.0)
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);STORE(RHO1,MACH)
    GROUP 8. Terms (in differential equations) & devices
TERMS(V1,P,P,N,P,P,P);TERMS(W1,P,P,N,P,P,P)
V1AD=GRND1 ; DENPCO=T
    GROUP 9. Properties of the medium (or media)
   ** Use Isentropic Density Law
RHO1=COMPRESS;RHO1A=RHOTOT/PTOTAL**RGAM;RHO1B=RGAM
RHO1C=0.;PRESS0=0.;DRH1DP=COMPRESS
    GROUP 11. Initialization of variable or porosity fields
FIINIT(P1)=PIN;FIINIT(W1)=WIN;FIINIT(RHO1)=RHOIN
    GROUP 13. Boundary conditions and special sources
INLET(IN,LOW,1,NX,1,NY,1,1,1,1)
VALUE(IN,P1,RHOIN*WIN);VALUE(IN,W1,WIN)
    GROUP 16. Termination of iterations
AIN=YVLAST;FLOWIN=WIN*RHOIN*AIN
RESREF(P1)=1.E-12*FLOWIN
RESREF(V1)=RESREF(P1)*WIN;RESREF(W1)=RESREF(V1)
LITER(P1)=20;LITHYD=10
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,0.7)
DTF=ZLEN/WIN
RELAX(V1,FALSDT,DTF);RELAX(W1,FALSDT,DTF)
    GROUP 18. Limits on variables or increments to them
    GROUP 19. EARTH Calls To GROUND Station
    GROUP 22. Spot-value print-out
IYMON=NY/2;NPLT=2;TSTSWP=-1
    GROUP 23. Field print-out and plot control
ITABL=2;NYPRIN=2;NZPRIN=NZ/4
IF(NZ.GT.1) THEN
+ IDISPA=2;IDISPB=1;IDISPC=NZ
ENDIF