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)
name(r1)=gas;name(r2)=fue;name(rs)=shad

    ** 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(gas )=rgin;fiinit(fue )=rsin;fiinit(shad)=rsin

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(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.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,SHAD,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