file name fireball (dbs  31.08.92)   27.01.94
 
  PHOTON USE
  ext;p1;;
 
  view z;  norm
  msg  Contours of reactedness
  gr ou z 1;  con rctd z 1 sh;int 25;  upause 3
  ext;p2;;
 
  view z;  norm
  msg  Contours of reactedness
  gr ou z 1;  con rctd z 1 sh;int 25;  upause 3
  ext;p3;;
 
  view z;  norm
  msg  Contours of reactedness
  gr ou z 1;  con rctd z 1 sh;int 25;  upause 3
  ext;p4;;
 
  view z;  norm
  msg  Contours of reactedness
  gr ou z 1;  con rctd z 1 sh;int 25;  upause 3
  ext;p5;;
 
  view z;  norm
  msg  Contours of reactedness
  gr ou z 1;  con rctd z 1 sh;int 25;  upause 3
  ext;phi;;
 
  view z;  norm
  msg  Contours of reactedness
  ext;phi;;;
 
  view z;  norm;  gr ou z 1
  msg  Contours of reactedness with velocity vectors
  con rctd z 1 sh;int 25
  msg  Velocity vectors
  vec z 1;
  msg     -
  msg Press e to END
  ENDUSE
 
    GROUP 1. Run title and other preliminaries
TEXT(Spherical Flame Propagation         
TITLE
  DISPLAY
  A sphere of hot gas spreads flame into a cold combustible mixture
  by diffusion, heat conduction and chemical reaction.
  Gravity also causes relative motion.
  The reaction is single-step, with rate dependent only on
  reactedness.
  The spherical geometry is contrived by use of special porosity
  settings.
  Dimensionless diffusion, reaction-rate and buoyancy parameters
  may be varied.
  PHOTON USE commands are supplied for display purposes.
  ENDDIS
REAL(U1IN,RHO1IN,RCTDIN,RHOREF,BUOSPD,BUOY,VALU,RATIO)
 
    GROUP 2. Transience; time-step specification
STEADY=F
GRDPWR(T,50,0.01,1.0)
    GROUP 3. X-direction grid specification
CARTES=F
GRDPWR(X,12,3.1416,1.0)
    GROUP 4. Y-direction grid specification
GRDPWR(Y,12,0.01,1.0)
    GROUP 5. Z-direction grid specification
ZWLAST=YVLAST
 
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,V1,H1)
NAME(H1)=RCTD
SOLUTN(P1,Y,Y,Y,N,N,N)
STORE(RHO1,EPOR,VPOR,NPOR)
 
    GROUP 8. Terms (in differential equations) & devices
  ** Cut out built-in enthalpy source (viscous dissipation)
TERMS(RCTD,N,Y,Y,Y,Y,Y)
  ** denpco=t introduces density variations into pressure-correction
     coefficients
DENPCO=T
 
    GROUP 9. Properties of the medium (or media)
TMP1=LINH;TMP1A=0.0;CP1=1.0
RHOREF=1.0
RHO1=RECSCAL;RHO1A=1.0/RHOREF;RHO1B=5.0
ENUL=0.001
REAL(DIFSPD,DELT,DELY,FACTOR)
DELT=LSTEP
DELT=TLAST/DELT
DELY=NY
DELY=YVLAST/DELY
DIFSPD=ENUL*DELT/DELY**2
mesgm(radial cell dimension, dely = :dely:
mesgm(time interval        , delt = :delt:
mesgm(dimensionless diffusion speed = enul*delt/dely**2
mesg(Its value is    :difspd:  If not OK, type required value
READVDU(DIFSPD,REAL,DIFSPD)
ENUL=DIFSPD*DELY**2/DELT
mesgm(dimensionless diffusion speed = :difspd:
 
    GROUP 11. Initialization of variable or porosity fields
FIINIT(V1)=0.0;FIINIT(RCTD)=0.0
REAL(RCTDNY)
RCTDNY=0.0
PATCH(SPHERE1,INIVAL,1,NX,1,NY/4,1,1,1,1)
INIT(SPHERE1,RCTD,0.0,1.0)
mesgm(propagation will be outward. If not OK press n
ans=y
READVDU(ANS,CHAR,Y)
IF(:ANS:.EQ.N) THEN
 PATCH(SPHERE1,INIVAL,1,NX,3*NY/4+1,NY,1,1,1,1)
 RCTDNY=1.0
 RHOREF=RHOREF/(1.0+RHO1B)
 mesgm(propagation is from outer hot gas to inner unburned
ENDIF
 
LSWEEP=10
 
    GROUP 13. Boundary conditions and special sources
PATCH(OUTER,NORTH,1,NX,NY,NY,1,1,1,LSTEP)
COVAL(OUTER,P1,FIXP,0.0);COVAL(OUTER,RCTD,ONLYMS,RCTDNY)
 
   ** A non-linear source of RCTD is present; for this
      purpose, the subroutine GXCHSO is called from GREX3,
      group13 by setting COefficient to POLYNOM (GRND7) in the COVAL
      command; the subroutine is entered only when the patch
      name begins with the character CHSO.
      CO=GRND7 selects the following option:
           COefficient=chsoa*RCTD**chsob
      VAL=1.0 enforces that the source falls to zero when
      RCTD equals unity.
PATCH(CHSOTERM,PHASEM,1,NX,1,NY,1,NZ,1,LSTEP)
COVAL(CHSOTERM,RCTD,POLYNOM,1.0)
CHSOA=1.0E5;CHSOB=6.0
     The source is therefore chsoa*(1.0-RCTD)*RCTD**chsob
REAL(REACRATE,REAC)
REACRATE=(DELY**2/ENUL)*CHSOA/(CHSOB+1.0)/(CHSOB+2.0)
mesga(dimensionless reaction speed is
mesg(   (dely**2/enul)*chsoa/(chsob+1.0)/(chsob+2.0)
mesgb(Its value is  :reacrate:  If not OK, insert value
READVDU(REACRATE,REAL,REACRATE)
CHSOA=REACRATE*(CHSOB+1.0)*(CHSOB+2.0)*ENUL/DELY**2
mesga(dimensionless reaction speed is :reacrate:
 
    ...BUOYANCY
BUOYB=-200.0; BUOYD=RHOREF; buoy=buoyb
mesgm(dimensionless buoyancy number = -buoyb*enul/dely**3
BUOSPD = -BUOYB*ENUL/DELY**3
mesg(Its value is = :buospd: If not OK type required value
READVDU(BUOSPD,REAL,BUOSPD)
BUOYB=-BUOSPD*DELY**3/ENUL
IF(BUOYB.EQ.0) THEN
 NX=2
 SOLUTN(U1,Y,N,N,N,N,N)
 BUOYANCY=SKIP
ELSE
 IF(BUOYB.LT.BUOY) THEN
 mesgm(Reduce time step? (y/n)
 READVDU(ANS,CHAR,N)
 IF(:ANS:.EQ.Y) THEN
+ RATIO=BUOY/BUOYB
+ mesgm(time steps reduced in ratio :ratio:
+ DELAY(500)
+ TLAST=TLAST*RATIO
+ DELT=DELT*RATIO
 ENDIF
 ENDIF
ENDIF
 
PATCH(BUOYANCY,PHASEM,1,NX,1,NY,1,1,1,LSTEP)
COVAL(BUOYANCY,U1,FIXFLU,DENSDIFF)
COVAL(BUOYANCY,V1,FIXFLU,DENSDIFF)
mesgm(dimensionless buoyancy number = :buospd:
 
    GROUP 15. Termination of sweeps
SELREF=T;RESFAC=0.1
    GROUP 16. Termination of iterations
LITER(U1)=20;LITER(V1)=20
    GROUP 17. Under-relaxation devices
RELAX(RCTD,LINRLX,0.25)
    GROUP 18. Limits on variables or increments to them
VARMIN(RCTD)=0.0;VARMAX(RCTD)=1.0
    GROUP 19. Data communicated by satellite to GROUND
  ** next setting introduces porosities having the effect of
     spherical coordinates
IPORIB=1
    GROUP 22. Spot-value print-out
IXMON=1;IYMON=1;ITABL=1
    GROUP 23. Field print-out and plot control
NXPRIN=NX/5;NTPRIN=LSTEP/4;NPLT=1
OUTPUT(P1,Y,Y,Y,Y,Y,Y);OUTPUT(RCTD,Y,Y,Y,Y,Y,Y)
OUTPUT(V1,Y,Y,Y,Y,Y,Y);OUTPUT(U1,Y,Y,Y,Y,Y,Y)
TSTSWP=-1
IDISPB  =1 ;IDISPC =50;IDISPD=LSTEP/5