```
text(1D coal + wood combustion
#cls
Title
DISPLAY
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 )

GROUP 1. Run title and other preliminaries

Fuel composition and consequences for stoichiometry

** cincl  & hincl are the mass fractions of carbon & hydrogen
in the coal
ENDDIS
#pause
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
HCHX
HSIN
HGIN
FLOG=2.25*FLOS

** burning rates. Note that woodburn is set below at chsoa
real(coalburn,woodburn,charburn)
coalburn=1.;woodburn=1.e2;charburn=1.e3

GROUP 3. X-direction grid specification
GRDPWR(X,20,10.,1.)

GROUP 7. Variables stored, solved & named
onephs=f
solve(p1,u1,u2,r1,r2,rs)

** 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)

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

RSIN=FLOS/(RHO2*VELS)
RGIN=1.-RSIN;VELG=FLOG/(RHOGIN*RGIN)

enut=0.01*xulast*velg

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.e5
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

fiinit( u1 )=0.0 ;fiinit( u2 )=0.0
fiinit( h1 )=hgin;fiinit( h2 )=hsin
fiinit(tmp1)=tgin;fiinit(tmp2)=tsin

fiinit(mdot)=0.01*flos;fiinit(den1)=rhogin;fiinit(mixf)=0.1

GROUP 13. Boundary conditions and special sources

** injectors
patch(inletsw,west,1,1,1,1,1,1,1,1)
coval(inletsw, p1 ,fixflu, flog)
coval(inletsw, p2 ,fixflu, flos)
coval(inletsw, u1 ,onlyms, velg)
coval(inletsw, u2 ,onlyms, vels)
coval(inletsw,mixf,onlyms, 0.0 )
coval(inletsw, h1 ,onlyms, hgin)
coval(inletsw, h2 ,onlyms, hsin)

** wood
patch(woodinsw,west,1,1,1,1,1,1,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)

**************************************************
****  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   ****
**********************************************

** Outlet at high end
real(outco1);outco1=1.e3;outco1=0.1
patch(outlet,east,nx,nx,1,ny,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,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)=30
liter(u1)= 2;liter(u2)= 2

GROUP 17. Underelaxation devices

chatim=0.0002

relax(gas ,linrlx,0.1);relax(cfip,linrlx,0.3)
relax(den1,linrlx,0.1);relax(mdot,linrlx,0.1)
relax(mixf,falsdt,.01)

relax(u1,falsdt,chatim)
relax(u2,falsdt,chatim)
relax(h1,falsdt,chatim*10)
relax(h2,falsdt,chatim*10)
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
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
1-d plots
*
PATCH(AXIS,PROFIL,1,NX,1,1,1,1,1,1)
COVAL(AXIS,MIXF,0,0.1)
COVAL(AXIS,FWD,0,0.1)
COVAL(AXIS,WOOD,0,0.1)
COVAL(AXIS,CHAR,0,0.1)

PATCH(AXIS2,PROFIL,1,NX,1,1,1,1,1,1)
COVAL(AXIS2,YO2,0,0)
COVAL(AXIS2,YCO2,0,0)
COVAL(AXIS2,YCO,0,0)
COVAL(AXIS2,YH2O,0,0)
COVAL(AXIS2,YH2,0,0)

PATCH(AXIS3,PROFIL,1,NX,1,1,1,1,1,1)
COVAL(AXIS3,FUE,0,0)
COVAL(AXIS3,MDOT,0,0)

PATCH(AXIS4,PROFIL,1,NX,1,1,1,1,1,1)
COVAL(AXIS4,H1,0,0)
COVAL(AXIS4,H2,0,0)
COVAL(AXIS4,TMP1,0,0)
COVAL(AXIS4,TMP2,0,0)

PATCH(AXIS5,PROFIL,1,NX,1,1,1,1,1,1)
COVAL(AXIS5,3,0,0)
COVAL(AXIS5,4,0,0)
COVAL(AXIS5,1,0,0)

ORSIZ=0.2;ABSIZ=0.6
TSTSWP=-1
VARMIN(FUE)=1.E-10
DENPCO=T
STOP
```