****** TO LOAD CASE 111:TYPE L(N111) *****
    GROUP 1. Run title
TEXT(1D SHOCK-FREE TRANSONIC NOZZLE FLOW
TITLE
  DISPLAY
   This case concerns steady 1d inviscid transonic flow through a
   laval nozzle ( M.R.Malin, MSc, London Univ., [1977] ). The nozzle
   geometry is designed to produce a linear distribution of the Mach
   number M under shock-free conditions. The flow is asymmetric
   about the throat with subsonic axial inflow and supersonic
   outflow. The nozzle length is X=3 and the throat area is unity.
   The inlet conditions are prescribed total pressure Po and total
   temperature To at M=0.5. The design outlet Mach number is 2.0 for
   which the exit pressure is 0.1278 times Po. Library case N112
   considers the same case, but with a back pressure of 0.744 times
   Po, implying a normal shock at X=2.4. Computations are made with
   the bounded higher-order UMIST scheme, and compressibility
   corrections are used so as to obtain acceptable shock predictions
   (see Malin & Sanchez, PHOENICS J., Vol.1, No.2, p214, [1988]).
   For shock-free flow the UMIST scheme offers no advantage, but it
   produces a much sharper shock than the UDS scheme.
  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*SQRT(GAMMA)/Ao. Ao is the acoustic velocity at To (see Palacio
   et al, Int.J.Heat Mass Transfer, Vol.33, No.6, p1193, [1990] ).
 
   For case N112, PHOENICS without the compressibility corrections
   produces a shock at the exit plane because the in-line convection
   flux in the momentum equation uses an averaged density rather
   than the upwind density. If the upwind density were used, as in
   GXHOCS, the shock would be predicted somewhat down of its
   physical location. The compressibility corrections arrange that
   this flux is computed by averaging the continuity fluxes, thereby
   ensuring that the shock is predicted in the correct location.
  PHOTON USE
  AUTOPLOT
  file
  phi 5
 
  cl;d 1 p1;d 1 pa;col9 1;blb2 2;level x 2.4
  msg computed and analytical pressure distribution
  msg press  to continue
  pause
  cl;d 1 mach;d 1 ma;col9 1;blb2 2;level x 2.4
  msg computed and analytical Mach-number distribution
  msg press  to continue
  pause
  ENDUSE
REAL(AEX,AIN,AREA,ATH,CMASS,CP)
REAL(RCON,GAM,GP,GM1,GP1,GRHO,G1,G2,G3,GDUDX,PTOT,MFEX,TTOT)
REAL(RHTOT,MI,RGM,MO,MX,PIN,TIN,POW1,WIN,RHIN,PEX,WEX,RHEX,DTF)
REAL(XCORD,VMASS,MFIN,XTH,XH,DMDX,MX2)
INTEGER(NUP,NDOWN);CHAR(DIRH,VELH,SCHM,PORA)
  ** Define various constants
GAM=1.4;GP1=GAM+1.;GM1=GAM-1.;RGM=1./GAM;POW1=GAM/GM1
G1=2./GP1;G2=GM1/GP1;G3=GP1/(2.0*GM1);RCON=1.0;CP=RCON*GAM/GM1
  ** Define inlet total conditions & Mach number
PTOT=1.0;TTOT=1.0;RHTOT=1.0;MI=0.5
  ** Define throat area, nozzle length,design oulet Mach number
     and throat area location
XCORD=3.0;ATH=1.0;MO=2.0;XTH=1.0
  ** Compute inlet and exit areas for design condition
AIN=((G1+G2*MI*MI)**G3)*ATH/MI
AEX=((G1+G2*MO*MO)**G3)*ATH/MO
  ** Calculation of inlet static values
PIN=PTOT/(1.+GM1*MI*MI/2.)**POW1;RHIN=RHTOT/(PTOT/PIN)**RGM
WIN=MI*(GAM*PIN/RHIN)**0.5;MFIN=RHIN*WIN*AIN
TIN=PIN/(RCON*RHIN)
  ** Calculation of outlet design static values
PEX=PTOT/(1.+GM1*MO*MO/2.)**POW1;RHEX=RHTOT/(PTOT/PEX)**RGM
WEX=MO*(GAM*PEX/RHEX)**0.5;MFEX=MFIN/AEX
    GROUP 5. Z-direction grid specification
MESG( Enter required coordinate direction X,Y or Z
MESG( Default:  X
READVDU(DIRH,CHAR,X)
CASE :DIRH: OF
WHEN X,1
+ VELH=U1;PORA=EPOR
WHEN Y,1
+ VELH=V1;PORA=NPOR
WHEN Z,1
+ VELH=W1;PORA=HPOR
ENDCASE
NREG:DIRH:=2;NUP=20;NDOWN=40
IREG:DIRH:=1;GRDPWR(:DIRH:,NUP,XTH,1.0)
IREG:DIRH:=2;GRDPWR(:DIRH:,NDOWN,XCORD-XTH,1.0)
DTF=0.5*XCORD/(WIN*N:DIRH:)
 
MESG(NOZZLE OPERATING AT DESIGN CONDITIONS
MESG(Expected Mexit=2.0 and Pexit/POin=0.128
    GROUP 7. Variables stored, solved & named
SOLVE(P1,:VELH:);SOLUTN(P1,Y,Y,Y,N,N,N)
STORE(RHO1,TMP1,MACH,PA,MA);STORE(:PORA:)
    GROUP 8. Terms (in differential equations) & devices
  ** Cut-out diffusion terms
TERMS(:VELH:,P,P,N,P,P,P)
MESG( Enter required convection scheme
MESG(  Default: UMIST Differencing Scheme
MESG( The alternative is:
MESG(  UDS - Upwind differencing scheme
READVDU(SCHM,CHAR,UMIST)
CASE :SCHM: OF
WHEN UDS,3
+ DIFCUT=0.
+ MESG(Upwind-differencing scheme
WHEN UMIST,3
+ MESG(UMIST differencing scheme
+ SCHEME(UMIST,ALL)
ENDCASE
    GROUP 9. Properties of the medium (or media)
  ** Perfect Gas Law
RHO1=IDEALGAS;RHO1B=1./RCON;TMP1=CONSTAGH;TMP1A=CP*TTOT;CP1=CP
PRESS0=0.0;DRH1DP=IDEALGAS
    GROUP 11. Initialization of variable or porosity fields
DMDX=(MO-MI)/XCORD;GDUDX=(WEX-WIN)/ZWLAST
INIADD=F
  ** Initial guess plus analytical solution for M and P1; strictly
     M and P1 stored at node, but for simplicity this is ignored
     and so M and P1 are computed at cell face.
DO KK=1,N:DIRH:
IF(:DIRH:.EQ.Z) THEN
+ PATCH(INA:KK:,INIVAL,1,NX,1,NY,:KK:,:KK:,1,LSTEP)
+ XH=ZFRAC(KK)*ZWLAST
ENDIF
IF(:DIRH:.EQ.Y) THEN
+ PATCH(INA:KK:,INIVAL,1,NX,:KK:,:KK:,1,NZ,1,LSTEP)
+ XH=YFRAC(KK)*YVLAST
ENDIF
IF(:DIRH:.EQ.X) THEN
+ PATCH(INA:KK:,INIVAL,:KK:,:KK:,1,NY,1,NZ,1,LSTEP)
+ XH=XFRAC(KK)*XULAST
ENDIF
+ MX=DMDX*XH+MI;MX2=MX*MX;AREA=((G1+G2*MX2)**G3)*ATH/MX
+ GP=PTOT/(1.+0.5*GM1*MX2)**POW1;GRHO=RHTOT*(GP/PTOT)**RGM
+ INIT(INA:KK:,:PORA:,ZERO,AREA);INIT(INA:KK:,P1,ZERO,GP)
+ INIT(INA:KK:,RHO1,ZERO,GRHO)
  ** Store analytical values
+ INIT(INA:KK:,PA,ZERO,GP);INIT(INA:KK:,MA,ZERO,MX)
ENDDO
    GROUP 13. Boundary conditions and special sources
  ** Total pressure inlet condition; neglects effect on momentum,
     i.e. the updating of WIN
VMASS=RHIN*PTOT/RHTOT;CMASS=2.*POW1*AIN/WIN
CASE :DIRH: OF
WHEN X,1
+ PATCH(INLET, CELL, 1, 1,1,NY,1,NZ,1,1)
+ PATCH(OUTLET,CELL,NX,NX,1,NY,1,NZ,1,1)
WHEN Y,1
+ PATCH(INLET, CELL,1,NX, 1, 1,1,NZ,1,1)
+ PATCH(OUTLET,CELL,1,NX,NY,NY,1,NZ,1,1)
WHEN Z,1
+ PATCH(INLET, CELL,1,NX,1,NY, 1, 1,1,1)
+ PATCH(OUTLET,CELL,1,NX,1, 1,NZ,NZ,1,1)
ENDCASE
COVAL(INLET,P1,CMASS,VMASS);COVAL(INLET,:VELH:,ONLYMS,WIN)
  ** Exit boundary condition
COVAL(OUTLET,P1,1.E4*MFEX*AEX/PEX,PEX)
COVAL(OUTLET,U1,ONLYMS,0.0);COVAL(OUTLET,:VELH:,ONLYMS,0.0)
    GROUP 15. Termination of sweeps
LSWEEP=500
    GROUP 16. Termination of iterations
LITER(P1)=15
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,1.0);DENPCO=T
RELAX(:VELH:,FALSDT,DTF)
    GROUP 18. Limits on variables or increments to them
VARMIN(RHO1)=0.1*RHIN;VARMAX(RHO1)=RHTOT;VARMIN(:VELH:)=-0.1
VARMAX(:VELH:)=50.;VARMIN(P1)=0.01*PIN;VARMAX(P1)=PTOT
    GROUP 21. Print-out of variables
OUTPUT(RHO1,Y,P,P,P,Y,Y)
    GROUP 22. Spot-value print-out
I:DIRH:MON=NUP+20
    GROUP 23. Field print-out and plot control
NPRINT=LSWEEP;TSTSWP=-1;ITABL=3;NPLT=10
    GROUP 23. Field print-out and plot control