Coal-Wood Combustion Model; Two-phase, with Slip.
Reactions: C (s) + 0.5 O2 > CO (exothermic )
CO + 0.5 O2 > CO2 (exothermic )
C(s) + CO2 > 2CO (endothermic)
C(s) + H2O > CO + H2 (endothermic)
H2 + 0.5 O2 > H2O (exothermic )
*asapbegin
*seegeo f
*asapend
trace=t
GROUP 1. Run title and other preliminaries
TEXT(coal- and wood-burning furnace
** from foot to meter
real(conv1); conv1=0.0254
** injection angles
real(pi,angle1,angle2);pi=3.14159
angle1=pi*33.5/180;angle2=pi*41.5/180
real(xcsw,ycsw,xcse,ycse,xcne,ycne,xcnw,ycnw)
xcsw=cos(angle1);ycsw=sin(angle1)
xcse=cos(angle2);ycse=sin(angle2)
xcne=cos(angle1);ycne=sin(angle1)
xcnw=cos(angle2);ycnw=sin(angle2)
Fuel composition and consequences for stoiciometry
** cincl & hincl are the mass fractions of carbon & hydrogen
in the coal
REAL(cincl,hincl,nincl)
cincl=0.86;hincl=0.05
In the following formulae:
0.232 is the mass of oxygen per unit mass of air
0.768 is the mass of nitrogen per unit mass of air
2.0, 12.0, 16.0, 18.0, 32.0 & 44.0 are molecular weights of
H2, C, O, H2O, O2 & CO2 respectively
REAL(FS,FS2)
** FS is the mass of fuel per unit mass of air/fuel mixture
to convert all carbon and oxygen to carbon monoxide.
FS=0.232/(0.232 + cincl*16.0/12.0)
** FS2 is the mass of fuel per unit mass of air/fuel mixture
to convert all carbon, hydrogen and oxygen to
carbon dioxide and water vapour.
FS2=0.232/(0.232 + CINCL*32.0/12.0 + HINCL*16./2.0)
Thermodynamic data
** hcco2 = heat of combustion for C + O2 --> CO2
hcco = " " " " C + 0.5 O2 --> CO
hhh2o = " " " " H2 + 0.5 O2 --> H2O
* the heat of reaction for C + O2 -> CO2, hcco2: 3.279E7
* the heat of reaction for C + 0.5*O2 -> CO , hcco : 9.208E6
* the heat of reaction for H2 + 0.5*O2 -> H2O, hhh2o: 1.209E6
* the specific heat at constant pressure, CP : 1.100E3
H = CP*T + HCHX*YCHX + hcoco2*YCO * HH2*YH2
REAL(HCCO2,HCCO,HHH2O,HCHX)
HCCO2=32.792E6; HCCO=9.208E6; HHH2O=120.9E6
HCHX=CINCL*HCCO2 + HINCL*HHH2O
HCHX
REAL(HGIN,GALF,HF,HA2,HO,HSIN)
** take cpsolid=cpgas=1.1e3
REAL(RHOGIN); CP1=1.1E3; CP2=CP1
Data concerning the inflows of fuel and air
REAL(FLOG,FLOS,VELS,VELG,CHATIM,RGIN,RSIN)
REAL(TGIN,TSIN)
** wood properties
real(hwood,hchar,hvol,hwdin,vlinwd,hinvl,cinvl,mlwtvl,moinwd)
vlinwd=0.5;hinvl=0.02;cinvl=1.0-hinvl;mlwtvl=100;moinwd=0.20
hvol =hhh2o*hinvl+hcco2*cinvl; hchar=hcco2
hwdin=cp1*tgin+hhh2o*vlinwd*hinvl+hcco2*(1-vlinwd+vlinwd*cinvl)
spedat(set,woodburn,hinvl,r,hinvl)
spedat(set,woodburn,vlinwd,r,vlinwd)
spedat(set,woodburn,mlwtvl,r,mlwtvl)
spedat(set,woodburn,moinwd,r,moinwd)
TGIN=500.;TSIN=500.;FLOS=1.0;VELS=1.
flos=20.0
vels= 2.0
HGIN=CP1*TGIN
HSIN=CP2*TSIN + HCHX
FLOG=2.25*FLOS
** burning rates. Note that woodburn is set below at chsoa
real(coalburn,woodburn,charburn)
coalburn=1.e4;woodburn=1.e2/5;charburn=1.e2/5
coalburn=1.;woodburn=1;charburn=1
GROUP 3. X-direction grid specification
GRDPWR(X,20,(32*12+11.00)*conv1,1.)
GROUP 4. Y-direction grid specification
GRDPWR(Y,20,(32*12+11.00)*conv1,1.)
GROUP 5. Z-direction grid specification
nregz=4
iregz=1;grdpwr(z,1,(14*12+ 7.00)*conv1,1.)
iregz=2;grdpwr(z,1,(10*12+4.875)*conv1,1.)
iregz=3;grdpwr(z,1,(18*12+ 9.00)*conv1,1.)
iregz=4;grdpwr(z,2,(23*12+6.125)*conv1,1.)
GROUP 6. Body-fitted coordinates or grid distortion
GROUP 7. Variables stored, solved & named
onephs=f
solve(p1,u1,u2,v1,v2,w1,w2,r1,r2,rs)
name(r1)=gas;name(r2)=fue;name(rs)=shad
solutn(p1,y,y,y,p,p,p)
** provide storage for inter-phase mass transfer etc.
store(mdot,cfip,den1)
** Solve additionally for the mixture fraction, i.e. the
quantity of phase-2 material which has entered phase 1.
solve(c1);name(c1)=wood;store(yco,yo2,yco2,yn2,yh2,yh2o)
solve(c3);name(c3)=fwd
solve(c5);name(c5)=mixf
solve(char);store(yvol)
** solve for enthalpy & store temperature
solve(h1,h2);store(tmp1,tmp2)
** k-e turbulence model
turmod(kemodl);store(enut);kelin=1
GROUP 8. Terms (in differential equations) & devices
eqdvdp=f;denpco=t;isolx=0;isoly=0;isolz=0
terms( h1 ,n,p,y,p,p,p);terms( h2 ,n,p,y,p,p,p)
terms(fwd ,n,p,y,p,p,y)
terms(wood,p,p,n,p,y,n);terms(char,p,p,p,p,y,n)
GROUP 9. Properties of the medium (or media)
rho1=7gases ; rho2=1.e3 ; press0=1.e5
temp0=0.0;rhogin=press0/(287.41*tgin)
rho1a=cincl ; rho1b=hincl
GROUP 10. Inter-phase-transfer processes and properties
** Set constant interphase friction factor and activate
the calculation of the interphase mass transfer by:
cfips=grnd1; cfipc=1.e8
cmdot=grnd3; cmdta=coalburn*100; cmdtc=fs
** Note that grnd3 is an mdot option, making the mass-
transfer rate proportional to (cmdtc - mixf), where
cmdtc stands for the saturation value of mixf, i.e. the
largest value which can be attained as a result of mass
transfer.
cint(mixf)=0.0 ; cint(fwd)=0.0
phint(h1)=7gases; phint(h2)=7gases
phint(mixf)=1.0 ; phint(fwd)=0.0
phint(wood)=0.0 ; phint(char)=0.0
GROUP 11. Initialization of variable or porosity fields
RSIN=FLOS/(RHO2*VELS)
RGIN=1.-RSIN;VELG=FLOG/(RHOGIN*RGIN)
fiinit(gas )=rgin
fiinit(fue )=rsin
fiinit(shad)=rsin
fiinit( u1 )=0.0 ;fiinit( u2 )=0.0
fiinit( v1 )=0.0 ;fiinit( v2 )=0.0
fiinit( w1 )=velg;fiinit( w2 )=vels
fiinit( h1 )=hgin;fiinit( h2 )=hsin
fiinit(tmp1)=tgin;fiinit(tmp2)=tsin
fiinit(mdot)=0.01*flos;fiinit(den1)=rhogin;fiinit(mixf)=0.1
real(kein,epin)
kein=0.0025*velg*velg
epin=0.1643*velg**1.5/(0.09*xulast/nx)
fiinit(ke)=kein
fiinit(ep)=epin
fiinit(enut)=0.001
GROUP 13. Boundary conditions and special sources
** injectors
patch(inletsw,west,1,1,1,1,#3,#3,1,1)
coval(inletsw, p1 ,fixflu, flog)
coval(inletsw, p2 ,fixflu, flos)
coval(inletsw, u1 ,onlyms, velg*xcsw)
coval(inletsw, u2 ,onlyms, vels*xcsw)
coval(inletsw, v1 ,onlyms, velg*ycsw)
coval(inletsw, v2 ,onlyms, vels*ycsw)
coval(inletsw,mixf,onlyms, 0.0 )
coval(inletsw, h1 ,onlyms, hgin)
coval(inletsw, h2 ,onlyms, hsin)
coval(inletsw, ke ,onlyms, kein)
coval(inletsw, ep ,onlyms, epin)
patch(inletse,south,nx,nx,1,1,#3,#3,1,1)
coval(inletse, p1 ,fixflu, flog)
coval(inletse, p2 ,fixflu, flos)
coval(inletse, v1 ,onlyms, velg*ycse)
coval(inletse, v2 ,onlyms, vels*ycse)
coval(inletse, u1 ,onlyms,-velg*xcse)
coval(inletse, u2 ,onlyms,-vels*xcse)
coval(inletse,mixf,onlyms, 0.0 )
coval(inletse, h1 ,onlyms, hgin)
coval(inletse, h2 ,onlyms, hsin)
coval(inletse, ke ,onlyms, kein)
coval(inletse, ep ,onlyms, epin)
patch(inletne,east,nx,nx,ny,ny,#3,#3,1,1)
coval(inletne, p1 ,fixflu, flog)
coval(inletne, p2 ,fixflu, flos)
coval(inletne, u1 ,onlyms,-velg*xcne)
coval(inletne, u2 ,onlyms,-vels*xcne)
coval(inletne, v1 ,onlyms,-velg*ycne)
coval(inletne, v2 ,onlyms,-vels*ycne)
coval(inletne,mixf,onlyms, 0.0 )
coval(inletne, h1 ,onlyms, hgin)
coval(inletne, h2 ,onlyms, hsin)
coval(inletne, ke ,onlyms, kein)
coval(inletne, ep ,onlyms, epin)
patch(inletnw,north,1,1,ny,ny,#3,#3,1,1)
coval(inletnw, p1 ,fixflu, flog)
coval(inletnw, p2 ,fixflu, flos)
coval(inletnw, v1 ,onlyms,-velg*ycnw)
coval(inletnw, v2 ,onlyms,-vels*ycnw)
coval(inletnw, u1 ,onlyms, velg*xcnw)
coval(inletnw, u2 ,onlyms, vels*xcnw)
coval(inletnw,mixf,onlyms, 0.0 )
coval(inletnw, h1 ,onlyms, hgin)
coval(inletnw, h2 ,onlyms, hsin)
coval(inletnw, ke ,onlyms, kein)
coval(inletnw, ep ,onlyms, epin)
** wood
patch(woodinsw,west,1,1,1,1,#3,#3,1,1)
coval(woodinsw, p1 ,fixflu,flog*0.01)
coval(woodinsw, fwd,onlyms,1.0)
coval(woodinsw,wood,onlyms,1.0)
coval(woodinsw,h1,onlyms,hwdin)
patch(woodinse,south,nx,nx,1,1,#3,#3,1,1)
coval(woodinse, p1 ,fixflu,flog*0.01)
coval(woodinse, fwd,onlyms,1.0)
coval(woodinse,wood,onlyms,1.0)
coval(woodinse,h1,onlyms,hwdin)
patch(woodinne,east,nx,nx,ny,ny,#3,#3,1,1)
coval(woodinne, p1 ,fixflu,flog*0.01)
coval(woodinne, fwd,onlyms,1.0)
coval(woodinne,wood,onlyms,1.0)
coval(woodinne,h1,onlyms,hwdin)
patch(woodinnw,north,1,1,ny,ny,#3,#3,1,1)
coval(woodinnw, p1 ,fixflu,flog*0.01)
coval(woodinnw, fwd,onlyms,1.0)
coval(woodinnw,wood,onlyms,1.0)
coval(woodinnw,h1,onlyms,hwdin)
**************************************************
**** Notes on combustion laws and constants ****
**************************************************
Each of (i) pyrolysis and (ii) char generation can be modelled either
as a power law or as an Arrhenius law. However, the same law must be
used for both.
Note also that currently only a single set of constants CHSOA ... CHSOE
are provided to cover both pyrolysis and char generation.
In the rate formulae below, "wood" means wood mass fraction, ditto for
"char", "co2", etc.
*********************
**** Pyrolysis ****
*********************
patch(chsow-,phasem,1,nx,1,ny,1,nz,1,lstep)
Pyrolysis (wood => volatiles) can be set either with a power law or an
Arrhenius relation (latter currently active).
Power law. Rate = -chsoa.chsoe.wood.T**chsod
=========
COVAL(CHSOW-,WOOD,GRND5,0.0)
Constants for power law.
CHSOA=1.
CHSOD=5.
CHSOE=10.0*(1.E-3)**CHSOD
CHSOB=0;CHSOC=0.0
Arrhenius. Rate = -chsoa.wood.e**(chsod/T)
=========
coval(chsow-,wood,grnd6,0.0)
Arrhenius rate constants.
real(acte,gascon)
chsoa=woodburn
Activation energy ACTE in kj/mol (hence GASCON div by 1000)
acte=2.50e+2
gascon=8.314
chsod=-acte/gascon
Multipliers for redundant terms in FN's in GXSOR (do not change these)
chsob=0.; chsoc=0.0; chsoe=0.
**************************
**** Char creation ****
**************************
patch(chsoc+,phasem,1,nx,1,ny,1,nz,1,1)
Char creation (wood => char) can be set either with a power law or an
Arrhenius relation (latter currently active). (This is similar to the
pyrolysis formation above, and the same choice must be used for both.)
Power law. Rate = -chsoa.chsoe.wood.T**chsod
=========
COVAL(CHSOC+,CHAR,FIXFLU,GRND1)
For constants see power law section of Pyrolysis (above).
Arrhenius. Rate = -chsoa.wood.*e**(chsod/T)
=========
coval(chsoc+,char,fixflu,grnd2)
spedat(set,chsoc+,react4,c,wood)
spedat(set,chsoc+,consta,r,chsoa*(1.0-vlinwd))
spedat(set,chsoc+,constb,r,0.0)
For constants see Arrhenius section of Pyrolysis (above).
*************************
**** Char burning ****
*************************
patch(chsoc-,phasem,1,nx,1,ny,1,nz,1,1)
Char burning is assumed to obey the following relation.
rate of burn = - char . ( C.co2 + D.h2o + E.o2)
The constants C, D, E are currently assumed to be equal and to have
the value CHARBURN.
Covals and spedats for CHSO1 patch.
coval(chsoc-,char,grnd1,0.0)
spedat(set,chsoc-,react1,c,yco2)
spedat(set,chsoc-,react2,c,yh2o)
spedat(set,chsoc-,react3,c,yo2)
spedat(set,chsoc-,constc,r,charburn)
spedat(set,chsoc-,constd,r,charburn)
spedat(set,chsoc-,conste,r,charburn)
**********************************************
**** End of combustion setting section ****
**********************************************
** walls
goto jump
patch(wallw,wwall,1,1,1,ny,1,nz,1,1)
coval(wallw,v1,grnd2,0.0)
coval(wallw,w1,grnd2,0.0)
coval(wallw,v2,grnd2,0.0)
coval(wallw,w2,grnd2,0.0)
coval(wallw,ke,grnd2,grnd2)
coval(wallw,ep,grnd2,grnd2)
patch(walle,ewall,nx,nx,1,ny,1,nz,1,1)
coval(walle,v1,grnd2,0.0)
coval(walle,w1,grnd2,0.0)
coval(walle,v2,grnd2,0.0)
coval(walle,w2,grnd2,0.0)
coval(walle,ke,grnd2,grnd2)
coval(walle,ep,grnd2,grnd2)
patch(walls,swall,1,nx,1,1,1,nz,1,1)
coval(walls,u1,grnd2,0.0)
coval(walls,w1,grnd2,0.0)
coval(walls,u2,grnd2,0.0)
coval(walls,w2,grnd2,0.0)
coval(walls,ke,grnd2,grnd2)
coval(walls,ep,grnd2,grnd2)
patch(walln,nwall,1,nx,ny,ny,1,nz,1,1)
coval(walln,v1,grnd2,0.0)
coval(walln,w1,grnd2,0.0)
coval(walln,v2,grnd2,0.0)
coval(walln,w2,grnd2,0.0)
coval(walln,ke,grnd2,grnd2)
coval(walln,ep,grnd2,grnd2)
patch(walll,lwall,1,nx,1,ny,1,1,1,1)
coval(walll,u1,grnd2,0.0)
coval(walll,v1,grnd2,0.0)
coval(walll,u2,grnd2,0.0)
coval(walll,v2,grnd2,0.0)
coval(walll,ke,grnd2,grnd2)
coval(walll,ep,grnd2,grnd2)
label jump
** Outlet at high end
real(outco1);outco1=1.e3;outco1=0.1
patch(outlet,high,6,16,6,16,nz,nz,1,1)
coval(outlet, p1,outco1*0.1 ,zero)
coval(outlet, p2,outco1*rho2,zero)
coval(outlet,mixf, onlyms ,same)
coval(outlet, h1 , onlyms ,same)
coval(outlet, h2 , onlyms ,same)
coval(outlet, ke , onlyms ,same)
coval(outlet, ep , onlyms ,same)
coval(outlet,wood, onlyms ,same)
coval(outlet,char, onlyms ,same)
coval(outlet,fwd , onlyms ,same)
GROUP 15. Termination of sweeps
LSWEEP=2000;SELREF=T;RESFAC=0.001
liter(p1)=100
endit(p1)=grnd1
liter(u1)= 2;liter(u2)= 2
liter(v1)= 2;liter(v2)= 2
liter(w1)= 2;liter(w2)= 2
GROUP 17. Underelaxation devices
relax(p1,linrlx,0.5)
relax(den1,linrlx,0.1)
relax(enut,linrlx,0.1)
chatim=0.001
relax(shad,linrlx,0.1);relax(fue ,linrlx,0.1)
relax(gas ,linrlx,0.1);relax(cfip,linrlx,0.3)
relax(den1,linrlx,0.1);relax(mdot,linrlx,0.3)
relax(enut,linrlx,0.3);relax(mixf,falsdt,.01)
relax(u1,falsdt,chatim)
relax(u2,falsdt,chatim)
relax(v1,falsdt,chatim)
relax(v2,falsdt,chatim)
relax(w1,falsdt,chatim)
relax(w2,falsdt,chatim)
relax(h1,falsdt,chatim*10)
relax(h2,falsdt,chatim*10)
relax(ke,linrlx,0.5)
relax(ep,linrlx,0.5)
relax(wood,falsdt,chatim*10)
relax(char,falsdt,chatim*10)
relax(fwd ,falsdt,chatim*10)
GROUP 18. Limits on variables or increments to them
varmax(mixf)=cmdtc;varmin(mixf)=0.0
varmax(den1)=10.;varmin(den1)=1.e-1
varmax(enut)=fiinit(enut)*10.0
varmin(enut)=fiinit(enut)*0.1
varmin(p1)=-press0+10
GROUP 21. Print-out of variables
output(mdot,y,y,y,y,y,y);output(fue ,y,n,n,n,n,n)
output(tmp1,y,n,n,n,n,n);output(mdot,y,n,n,n,n,n)
GROUP 22. Spot-value print-out
ixmon=nx/2; iymon=ny/2; tstswp=-1
izmon=nz/2; nzprin=2;nxprin=1;nyprin=1
GROUP 23. Field print-out and plot control
if(nz.eq.1) then
solutn(w1,n,n,n,n,n,n)
solutn(w2,n,n,n,n,n,n)
endif
restrt(all)
LIBREF=115
STOP