PHOTON USE
  p
 
 
 
 
  gr x 1
  msg Grid
  vec x 1 sh
  MSG Velocity vectors
  msg Press return to redraw
  pause
  gr cl; gr ou x 1; red
  msg Press return to plot pressure contours
  pause
  cont p1 x 1 fil;.01
  msg Type e to End
  ENDUSE
 
    GROUP 1. Run title
TEXT(External flow over a rocket:        B529
TITLE
  DISPLAY
   This case represents supersonic flow along an axisymmetric
  body of rotation whose axis is parallel to the free stream.
  The grid is generated using transfinite interpolation.
  Density is calculated from temperature assuming constant
  stagnation enthalpy and is thus valid for adiabatic
  non-isentropic flow.
  ENDDIS
REAL(Z1,DZ1,DZ2,DZ3,RR1,RR2,RR3,RR4,PWRY,PZ1,PZ2,PZ3,PZ4,ZZ,HJ,DJ)
REAL(PZ5,STH,CTH,Z5)
INTEGER(K1,K2,K3,K4,KJ)
STH=0.2;CTH=(1.0-STH*STH)**0.5
REAL(GA,RGAS,AM,WMOL,DENS,RTEMP,CP,WVEL,COEF)
 
    GROUP 6. Body-fitted coordinates or grid distortion
BFC=T;NONORT=T
GSET(D,1,14,39,1.0,2.0,4.0)
PWRY=1.7
PZ1=1.8;PZ2=1.5;PZ3=1.0;PZ4=1.0;PZ5=1.5
 
     Z1     Z-coordinate of nose from start of domain
     DZ1    Z-distance from nose to large end of first cone
     DZ2    Length of second cone
     DZ3    Length of third cone
     RR1    Nose radius
     RR2    Base radius of first cone
     RR3    base radius of second cone
     RR4    base radius of third cone and cylinder
Z1=0.5;DZ1=.2874;DZ2=0.9;DZ3=0.4;RR1=0.03;RR2=0.1;RR3=0.2;RR4=0.35
Z5=0.2271
 
   ** Coordinates of join between nose and first cone; from front
DJ=0.03
K1=8;K2=15;K3=23;K4=29;KJ=10
GSET(P,P0,0.0,0.0,0.0)
GSET(P,P5,0.0,0.0,Z5)
   ** Set nose position
GSET(P,P1,0.0,0.0,Z1)
   ** set join between nose and first cone
GSET(P,PJ,0.0,RR1,Z1+DJ)
ZZ=Z1+DZ1;GSET(P,P2,0.0,RR2,ZZ)
ZZ=ZZ+DZ2;GSET(P,P3,0.0,RR3,ZZ)
ZZ=ZZ+DZ3;GSET(P,P4,0.0,RR4,ZZ)
GSET(P,PN,0.0,RR4,ZWLAST)
   ** Set all points along south boundary
GSET(L,L05,P5,P0,2,:PZ1:)
GSET(L,L51,P1,P5,K1-3,:PZ1:)
GSET(P,P7,0.0,0.03*0.707,0.53-0.03*0.707)
GSET(L,L1J,P1,PJ,KJ-K1,1.0,ARC,P7)
GSET(L,LJ2,PJ,P2,K2-KJ,:PZ2:)
GSET(L,L23,P2,P3,K3-K2,:PZ3:)
GSET(L,L34,P3,P4,K4-K3,:PZ4:)
GSET(L,L4N,P4,PN,NZ+1-K4,:PZ5:)
   ** Set points along low and high boundaries
GSET(P,PA,0.0,YVLAST,0.0);GSET(P,PB,0.0,YVLAST,ZWLAST)
GSET(P,P6,0.0,YVLAST,Z5)
GSET(L,L0A,P0,PA,NY,:PWRY:);GSET(L,LNB,PN,PB,NY,:PWRY:)
GSET(L,LA6,P6,PA,2,:PZ1:);GSET(L,L6B,P6,PB,NZ-2,1.0)
GSET(L,L56,P5,P6,NY,:PWRY:)
 
GSET(F,F1,P0,-,PA,-,P6,-,P5,-)
GSET(M,F1,+J+K,1,1,1,TRANS)
GSET(F,F2,P5,-,P6,-,PB,-,PN,P4.P3.P2.PJ.P1)
GSET(M,F2,+J+K,1,1,3,LAP10.FFFTFT)
   ** Adjust near axis points at entry
GSET(T,J3,F,J1,1,1,1,2,1.0)
   ** Transform to 2D polar domain
GSET(C,I2,F,I1,1,NY,1,NZ,RZ,-0.2,0.0,0.0)
 
    GROUP 7. Variables stored, solved & named
SOLVE(P1,V1,W1);SOLUTN(P1,Y,Y,Y,N,N,N);STORE(DEN1)
NCRT=2;STORE(VCRT,WCRT)
    GROUP 9. Properties of the medium (or media)
   ** Ratio of specific heats
GA=1.4
   ** Gas constant
RGAS=8305.6
   ** Molecular weight
WMOL=29.
   ** Specific heat
CP=1000.
   ** Free-stream Mach number
AM=1.2
   ** F.s. static temperature
RTEMP=288.
   ** Datum pressure (free-stream)
PRESS0=1.E5
   ** Density from perfect gas law in GREX3
RHO1=IDEALGAS
   ** Local static temperature from GREX3
TMP1=CONSTAGH
   ** DRH1DP for compressible adiabatic flow.
DRH1DP=IDEALGAS; RHO1C=1./GA
   ** F.s. density
DENS=PRESS0*WMOL/(RGAS*RTEMP)
   ** Stagnation enthalpy
TMP1A=RTEMP*CP*(1.+(GA-1.)*AM*AM*.5)
CP1=CP
   ** Constant in ideal gas law
RHO1B=WMOL/RGAS; RHO1C=1./GA
   ** Free stream velocity
WVEL=AM*(GA*RGAS*RTEMP/WMOL)**.5
   ** Kinematic viscosity
ENUL=1.34E-5;ENUT=2.E-4
    GROUP 11. Initialization of variable or porosity fields
FIINIT(W1)=WVEL;FIINIT(DEN1)=DENS
    GROUP 13. Boundary conditions and special sources
   ** Small-wave theory b.c for supersonic flow
COEF=(AM*AM-1.)**.5/WVEL
PATCH(NB,NORTH,#1,#NREGX,#NREGY,#NREGY,#1,#NREGZ,1,1)
COVAL(NB,P1,COEF,0.0)
COVAL(NB,W1,ONLYMS,WVEL);COVAL(NB,V1,ONLYMS,SAME)
 
WALL (SB,SOUTH,#1,#NREGX,#1,#1,#3,#NREGZ,1,1)
 
INLET(BFCIN,LOW,#1,#NREGX,#1,#NREGY,#1,#1,1,1)
VALUE(BFCIN,P1,GRND1);VALUE(BFCIN,W1,GRND1)
VALUE(BFCIN,WCRT,WVEL)
  *  Transfer density for GXBFC subroutine
BFCA=DENS
 
PATCH(OUTLET,HIGH,#1,#NREGX,#1,#NREGY,#NREGZ,#NREGZ,1,1)
      `Weak` exit boundary condition for supersonic flows.
  This allows an oblique shock to cross the exit plane by permitting
  a pressure gradient normal to the flow. This is done by making the
  outflow from each cell at the exit boundary equal to that from its
  low-z neighbour.
COVAL(OUTLET,P1,FIXFLU,GRND1)
COVAL(OUTLET,V1,ONLYMS,0.)
COVAL(OUTLET,W1,ONLYMS,WVEL)
 
    GROUP 15. Termination of sweeps
LSWEEP=100
   ** At least 150 sweeps will be required for a converged solution.
 
    GROUP 16. Termination of iterations
LITER(P1)=10;LITER(U1)=1;LITER(V1)=1;LITER(W1)=1
    GROUP 17. Under-relaxation devices
RELAX(P1,LINRLX,.35)
RELAX(V1,FALSDT,.002)
RELAX(W1,FALSDT,.0005)
RELAX(DEN1,LINRLX,1.)
    GROUP 19. Data communicated by satellite to GROUND
USEGRD=F
    GROUP 21. Print-out of variables
OUTPUT(P1,Y,Y,Y,Y,Y,Y)
OUTPUT(DEN1,Y,N,N,N,N,N)
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=3;IZMON=28;TSTSWP=-1
SELREF=T; RESFAC=0.01
    GROUP 23. Field print-out and plot control
NPRINT=LSWEEP;NYPRIN=1;NZPRIN=2
PATCH(YZ,CONTUR,1,1,1,NY,1,NZ,1,1)
PLOT(YZ,P1,0.0,20.0)