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
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
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)
if(:ans:.eq.y) then
sponig=t
goto skip
endif
mesg(Do you want to simulate the fully-stirred reactor? (y/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
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