GROUP 1. Run title
TEXT(De Laval Nozzle  X Direction      
TITLE
  DISPLAY
    Transonic flow in a convergent-divergent nozzle is
  considered. The throat is situated at one-third of the
  length of the duct. The geometry of the nozzle is:
   convergent section, inlet area = 1.34; throat area= 1.0
   divergent section, outlet area = 1.65.
  The inlet Mach number is 0.5, and the design exit Mach
  number is 0.6 for which Pexit/Po =0.784.
  Consequently, a standing shock wave is created in the
  divergent section.
    The alteration of the cross-sectional area as a function
  of x is provided by means of the east-face porosity factor.
  ENDDIS
 
   Locally-defined parameters:
 
   GAMMA       Specific heats ratio
   GASCON      Perfect gas constant
   P0          Stagnation pressure
   T0          Stagnation temperature
   RHO0        Stagnation density
   MACHIN      Inlet MACH number
   RHOIN       Inlet density
   CIN         Inlet sound speed
   VELIN       Inlet gas velocity
   PIN         Inlet pressure
   PEXIT       Exit pressure
 
INTEGER(IXTHR);REAL(GAMMA,GASCON,PEXIT,GAM1,POWER)
REAL(TINDT0,RHO0,T0,P0,RHOIN,MACHIN,PIN,VELIN,CIN)
REAL(CMASS,VMASS)
 
   ** Gas constants
GAMMA=1.4;GASCON=287.0
GAM1=GAMMA-1.0;POWER=GAMMA/(GAMMA-1.0)
   ** Total pressure, temperature, density, inlet Mach Number
P0=1.0; T0=1.0;RHO0=P0/(GASCON*T0); MACHIN=0.5
   ** Inlet sound speed & velocity
CIN=(GAMMA*P0/RHO0)**.5; VELIN=MACHIN*CIN
   ** Inlet static pressure & density
TINDT0=1./(1.+GAM1/2.*MACHIN**2)
RHO1B=1./GAMMA; RHO1A=RHO0/(P0**RHO1B)
PIN=P0*TINDT0**POWER; RHOIN=RHO1A*PIN**RHO1B
   ** Exit pressure
PEXIT=P0*0.784
 
    GROUP 3. X-direction grid specification
GRDPWR(X,60,1.0,1.0)
 
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1);STORE(RHO1,EPOR)
 
    GROUP 8. Terms (in differential equations) & devices
   ** Remove viscous terms
TERMS(U1,P,P,N,P,P,P)
 
    GROUP 9. Properties of the medium (or media)
   ** Select isentropic gas law
RHO1=COMPRESS;DRH1DP=COMPRESS
 
    GROUP 11. Initialization of variable or porosity fields
INIADD=F
FIINIT(U1)=2.0*VELIN;FIINIT(RHO1)=0.5*RHOIN
FIINIT(P1)=0.5*(P0+PEXIT)
   ** Convergent section: inlet area = 1.34; throat area= 1.0
PATCH(CONV,LINVLX,1,NX/3,1,1,1,1,1,1)
INIT(CONV,EPOR,-1.0,1.34)
   ** Divergent section: outlet area = 1.65
PATCH(DIV,LINVLX,NX/3+1,NX,1,1,1,1,1,1)
INIT(DIV,EPOR,1.0,1.0)
 
    GROUP 13. Boundary conditions and special sources
   ** The prescription of the stagnation pressure at the
    nozzle inlet is obtained from the following relationship
    between the stagnation pressure and the inlet values
    of Mach number and density:
    p0=rhoin*[1+0.5*(gamma-1)*Min**2]**[gamma/(gamma-1)].
      When the inlet mass-flow per unit area is made the
    subject of this equation, the familiar co & val form
    results:
    rhoin*uin=co*(val-pP), where,
    co=2*gamma/[uin*(gamma-1)], val=p0*rhoin/rho0.
      For this case, the inlet velocity which appears in co
    ie. uin is assumed known. Strictly speaking this
    quantity is a function of the calculation, and should
    hence be re-calculated in GROUND as the calculation
    proceeds.
 
CMASS=2.*POWER/VELIN;VMASS=RHOIN*P0/RHO0
PATCH(INLET,WEST,1,1,1,1,1,1,1,1)
COVAL(INLET,P1,CMASS,VMASS)
COVAL(INLET,U1,ONLYMS,VELIN)
 
   ** The static pressure is prescribed at the outlet
PATCH(OUTLET,EAST,NX,NX,1,1,1,1,1,1)
COVAL(OUTLET,P1,10.,PEXIT)
COVAL(OUTLET,U1,ONLYMS,0.0)
 
    GROUP 15. Termination of sweeps
LSWEEP=100
 
    GROUP 17. Under-relaxation devices
RELAX(U1,FALSDT,0.05*XULAST/VELIN)
 
    GROUP 22. Spot-value print-out
IXMON=NX-5;TSTSWP=LSWEEP/10;ITABL=1
 
    GROUP 23. Field print-out and plot control
IXPRF=NX/2;NXPRIN=NX/10;ORSIZ=.4
PATCH(PLOT,PROFIL,1,NX,1,1,1,1,1,1);PLOT(PLOT,P1,0,0)
PLOT(PLOT,U1,0,0)