DISPLAY
 
  FLAME PHENOMENA (Propagation, ignition in a gas)
 
  2-dimensional (x-y), cartesian, steady, elliptic simulation
 
  Demonstration of some steady one-phase flame phenomena
  namely   :-  1 D spontaneous ignition in a long duct
               1 D laminar flame propagation
 
    ----------------------------------------------------
    r=1.e-4        <--| flame-front                    r=1.
    -----------------------------------------------------
                          x----->
 
  In order to enable parametric studies to be made quickly,
  and to facilitate display, each problem is made 2 D,
  the additional dimensions being used for the varied
  parameters.
 
 
  enddis
 
  PHOTON USE
  p
 
 
 
  msg Reactedness contours
  con rctd z 1 fi;0.005
  msg( Press RETURN to continue
  pause
  *msg Velocity vectors. Press RETURN for grid
  *set ref vec 10;vec z 1
  *pause
  msg Grid
  gr z 1
  msg( Press RETURN for velocity contours
  pause
  con off;vec off;red
  msg Contours of gas velocity
  con u1 z 1 fi;0.01
  msg Press e to END
  enduse
 
       GROUP 1. Run title and other preliminaries
TEXT(Steady 1D Flame Phenomena         
TITLE
mesg(PC486/50 time last reported as appx. 2. min
    GXCHSO is employed for this calculation (see group 13
    of GREX3).
    -----------------------------------------------------
    r=1.e-4        <--| flame-front                    r=1.
    -----------------------------------------------------
                          x----->
real(u1in,rho1in,rctdin)
BOOLEAN(propag,sponig,stirred)
propag=f;sponig=f;stirred=f
mesg(Do you want to simulate spontaneous ignition? (y/n)
readvdu(ans,char,n)
if(:ans:.eq.y) then
 sponig=t
 goto skip
endif
  mesg(Do you want to simulate the fully-stirred reactor? (y/n)
  readvdu(ans,char,n)
  if(:ans:.eq.y) then
   stirred=t
   goto skip
   endif
mesg(Steady laminar flame propagation will be simulated
propag=t
label skip
propag
sponig
   stirred
u1in=1.0;rho1in=1.0;rctdin=0.01
if(sponig) then
 rctdin=0.2
endif
    GROUP 2. Transience; time-step specification
    GROUP 3. X-direction grid specification
GRDPWR(X,100,1.0,1.0)
  if(stirred) then
   nx=2
  endif
    GROUP 4. Y-direction grid specification
grdpwr(y,10,0.5,1.0)
    GROUP 5. Z-direction grid specification
    GROUP 6. Body-fitted coordinates or grid distortion
    GROUP 7. Variables stored, solved & named
SOLVE(P1,U1,H1)
if(ny.gt.1) then
 solve(v1)
 store(vpor,epor,npor)
endif
NAME(H1)=RCTD
store(enul,tmp1,rho1)
    GROUP 8. Terms (in differential equations) & devices
  ** Cut out built-in enthalpy source (viscous dissipation)
TERMS(RCTD,N,Y,Y,Y,Y,Y)
    GROUP 9. Properties of the medium (or media)
TMP1=LINH;tmp1a=0.0;cp1=1.0

if(propag) then
 ENUL=LINTEM;enula=0.2;enulb=0.6
endif
rho1=recscal;rho1a=1.0/rho1in;rho1b=5.0
    GROUP 10. Inter-phase-transfer processes and properties
    GROUP 11. Initialization of variable or porosity fields
fiinit(u1)=u1in;fiinit(v1)=0.0;fiinit(rctd)=rctdin
if(ny.gt.1) then
 fiinit(npor)=0.0
 iniadd=f
 patch(poros,linvly,1,nx,1,ny,1,1,1,lstep)
 init(poros,vpor,0.7,0.12);init(poros,epor,0.7,0.12)
endif
  if(stirred) then
   fiinit(rctd)=0.9
   fiinit(u1)=0.0;fiinit(v1)=0.0
  endif
    GROUP 12. Patchwise adjustment of terms in PDEs
    GROUP 13. Boundary conditions and special sources
   ** RCTD is held to unity at IX=NX
integer(nxin)
nxin=nx
  if(stirred) then
   nxin=1
  endif
integer(nyout)
nyout=ny
PATCH(IXEQNX,CELL,NXin,NX,1,nyout,1,1,1,LSTEP)
if(propag) then
+ COVAL(IXEQNX,RCTD,FIXVAL,1.0)
endif
 COVAL(IXEQNX,P1,FIXP,0.0)
   ** Inflow value of RCTD is specified at IX=1
integer(nyin)
nyin=ny
nxin=1
  if(stirred) then
   nxin=nx
  endif
 PATCH(IXEQ1,west,1,nxin,1,nyin,1,1,1,LSTEP)
 coval(ixeq1,p1,fixflu,u1in*rho1in)
 COVAL(IXEQ1,RCTD,onlyms,rctdin)
  if(.not.stirred) then
 coval(ixeq1,u1,onlyms,u1in)
  endif
   ** 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,1,1,LSTEP)
COVAL(CHSOTERM,RCTD,POLYNOM,1.0)
CHSOA=5.0e3;CHSOB=6.0
     The source is therefore CHSOA*(1.0-RCTD)*RCTD**CHSOB
    GROUP 15. Termination of sweeps
LSWEEP=100
    GROUP 16. Termination of iterations
    GROUP 17. Under-relaxation devices
relax(v1,falsdt,1.e-20)
relax(rho1,linrlx,0.5)
relax(p1,linrlx,0.5);relax(u1,falsdt,0.02*xulast/u1in)
  if(stirred) then
   relax(u1,falsdt,1.e-20)
  endif
    GROUP 18. Limits on variables or increments to them
    GROUP 22. Spot-value print-out
ixmon=nx/2;iymon=ny/2;itabl=1
    GROUP 23. Field print-out and plot control
NXPRIN=NX/5;NTPRIN=LSTEP/4;NPLT=1
   ** Plot profiles over the length of the domain
PATCH(XWISE,PROFIL,nx/2,NX,1,1,1,1,1,LSTEP)
PLOT(XWISE,RCTD,0.0,0.0)
if(ny.gt.1) then
 patch(map,contur,1,nx,1,ny,1,1,1,lstep)
 plot(map,rctd,0.0,10.0)
endif
xulast
enul
chsoa
real(reactno,peclet)
reactno=chsoa*xulast/(rho1in/u1in)
if(propag) then
 peclet=u1in*xulast*prndtl(rctd)/enula
  * Print to screen some interesting dimensionless quantities
 mesg(Peclet no = :peclet:
 mesg("reaction no"  chsoa*xulast/(rho1in*u1in) = :reactno:
endif
tstswp=-1;selref=t; resfac=1.e-2
libref=851