talk=t;run(1,1)
  PHOTON USE
  AUTOPLOT
  upause 1
  file; phi 5

  cl; da 1; yco2; da 1; yo2; da 1; yco; scale; cola 1; col3 2
  colf 3
  msg mass fractions of CO2 - yellow, O2 - blue, CO - red   ....more
  pause; cl; da 1; yh2o; da 1; mixf; scale; col3 1; colf 2
  msg mass fractions of H2O - blue, coal-derived gas - red  ....more
  pause; cl; da 1; tmp1; da 1; tmp2; scale; colf 1; cola 2
  msg temp of gas - red, of coal - yellow.                  ....more
  pause; cl; da 1; gvel z 1 19; da 1; cvel z 1 19; col3 1; col8 2
  msg gas velocity - blue, coal velocity - green.
  pause; end
  enduse

  DISPLAY
   Coal-fines combustion model; two-phase with slip, 1D,
                                one space, interphase heat
                                and mass transfer.

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

  Fuel composition and consequences for stoichiometry
  ** CINCL & HINCL are the mass fractions of carbon & hydrogen
     in the coal
REAL(CINCL, HINCL, NINCL)
CINCL=1.0; HINCL=0.1;NINCL=1.-CINCL-HINCL
     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)
  ** 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)

REAL(HCCO2,HCCO,HHH2O,HCHX,HCOCO2)
 HCCO2=32.792E6; HCCO=9.208E6; HHH2O=120.9E6
 HCOCO2=(12.0/28.0)*(HCCO2-HCCO)
 HCHX=CINCL*HCCO2 + HINCL*HHH2O

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,VELO,VELG,CHATIM,LENGTH,RGIN,RSIN)
REAL(TGIN,TSIN,BURNRATE)
TGIN=500.; TSIN=350.; FLOS=1.0; VELO=1.;LENGTH=10.0
BURNRATE=10000.
HGIN=CP1*TGIN
HSIN=CP2*TSIN + HCHX

FLOG=2.25*FLOS
TEXT(1D Fine coal-particle combustion
TITLE
NZ=20; GRDPWR(Z,NZ,LENGTH,1.)
ONEPHS=F
    GROUP 7. Variables stored, solved & named
SOLUTN(1,Y,Y,Y,P,P,P); OUTPUT(1,Y,Y,Y,Y,Y,Y)
BOOLEAN(NXNYNZ1);NXNYNZ1=T
IF(.NOT.ONEPHS) THEN
 SOLUTN(9,Y,Y,n,P,P,P); OUTPUT(9,Y,Y,Y,Y,Y,Y)
 SOLUTN(10,Y,Y,N,P,P,P); OUTPUT(10,Y,Y,Y,Y,Y,Y)
ENDIF
IF(NX.GT.1) THEN
NXNYNZ1=F
 SOLUTN(3,Y,Y,N,P,P,P); OUTPUT(3,Y,Y,Y,Y,Y,Y)
 IF(.NOT.ONEPHS) THEN
+  SOLUTN(4,Y,Y,N,P,P,P); OUTPUT(4,Y,Y,Y,Y,Y,Y)
 ENDIF
ENDIF
IF(NY.GT.1) THEN
NXNYNZ1=F
 SOLUTN(5,Y,Y,N,P,P,P); OUTPUT(5,Y,Y,Y,Y,Y,Y)
 IF(.NOT.ONEPHS) THEN
+  SOLUTN(6,Y,Y,N,P,P,P); OUTPUT(6,Y,Y,Y,Y,Y,Y)
 ENDIF
ENDIF
IF(NZ.GT.1) THEN
NXNYNZ1=F
+ SOLUTN(7,Y,Y,N,P,P,P); OUTPUT(7,Y,Y,Y,Y,Y,Y)
 IF(.NOT.ONEPHS) THEN
+ SOLUTN(8,Y,Y,N,P,P,P); OUTPUT(8,Y,Y,Y,Y,Y,Y)
 ENDIF
ENDIF
IF(NXNYNZ1) THEN
mesg(nx=ny=nz=1, so no velocities stored.
ENDIF

NAME(7)=GVEL; NAME(8)=CVEL
    GROUP 7. Variables stored, solved & named
SOLUTN(11,Y,Y,N,P,P,P)
NAME(9)=GAS; NAME(R2)=FUEL; NAME(11)=SHAD
  ** provide storage for inter-phase mass transfer
STORE(MDOT,CFIP)
  ** Solve additionally for the mixture fraction, i.e. the quantity
     of phase-2 material which has entered phase 1.
SOLUTN(C1,Y,Y,Y,P,P,P); NAME(C1)=MIXF
STORE(YCO,YO2,YCO2,YN2,YH2,YH2O,RHO1,RHO2)
SOLUTN(H1,Y,Y,N,P,P,P); SOLUTN(H2,Y,Y,N,P,P,P)
STORE(TMP1,TMP2)

    GROUP 8. Terms (in differential equations) & devices
TERMS(H1,N,P,P,P,Y,Y); TERMS(H2,N,P,P,P,N,Y)

    GROUP 9. Properties of the medium (or media)
RHO1= 1.2
RHO2=1.E3 ; PRESS0=1.E5

    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

    GROUP 11. Initialization of variable or porosity fields
RSIN=FLOS/(RHO2*VELO);  RHOGIN=PRESS0/(287.41*TGIN)
RGIN=1.-RSIN;           VELG=FLOG/(RHOGIN*RGIN)
FIINIT(GAS)=RGIN;       FIINIT(FUEL)=RSIN; FIINIT(SHAD)=RSIN
FIINIT(GVEL)=VELG;      FIINIT(CVEL)=VELO; FIINIT(MDOT)=0.01*FLOS
FIINIT(RHO1)=RHOGIN;    FIINIT(C4)=0.0;    FIINIT(C6)=0.0
FIINIT(H1)=HGIN;        FIINIT(H2)=HSIN;   FIINIT(MIXF)=0.
FIINIT(TMP1)=TGIN;FIINIT(TMP2)=TSIN
    GROUP 13. Boundary conditions and special sources
REAL(OUTCO1); OUTCO1=1.E3
REAL(MULT);MULT=2.0

PATCH(INLETS,LOW,1,1,1,1,1,1,1,1); PATCH(INLETG,LOW,1,1,1,1,1,1,1,1)
PATCH(OUTLET,CELL,1,1,1,1,NZ,NZ,1,1)

COVAL(INLETS,CVEL,ONLYMS,VELO); COVAL(INLETS,H2,ONLYMS,HSIN)
COVAL(INLETS,P2,FIXFLU,FLOS)  ; COVAL(INLETG,P1,FIXFLU,FLOG)
  gas phase
COVAL(INLETG,GVEL,ONLYMS,VELG); COVAL(INLETG,MIXF,ONLYMS,0.)
COVAL(INLETG,H1,ONLYMS,HGIN)
  ** multiply OUTCO1 by estimated value of exit density
COVAL(OUTLET,P1,OUTCO1*0.1,0.0); COVAL(OUTLET,P2,OUTCO1*RHO2,0.0)
    GROUP 15. Termination of sweeps
LSWEEP=2000; SELREF=T; RESFAC=0.001

    GROUP 17. Under-relaxation devices
CHATIM=MULT*LENGTH/(VELG*NX*NY*NZ)

RELAX(P1,LINRLX,0.3);   RELAX(SHAD,LINRLX,0.3)
RELAX(FUEL,LINRLX,0.3); RELAX(GAS,LINRLX,0.3)

RELAX(GVEL,FALSDT,0.001); RELAX(CVEL,FALSDT,0.001)

RELAX(MIXF,FALSDT,0.01)
RELAX(H1,FALSDT,0.1); RELAX(H2,FALSDT,0.1)

RELAX(RHO1,LINRLX,0.5)

    GROUP 18. Limits on variables or increments to them
VARMIN(MIXF)=0.; VARMAX(MIXF)=FS
VARMIN(YO2) =0.; VARMAX(YO2) =1.0
VARMIN(YCO) =0.; VARMAX(YCO) =1.0
VARMIN(YCO2)=0.; VARMAX(YCO2)=1.0
VARMIN(YH2O)=0.; VARMAX(YH2O)=1.0
VARMIN(YH2) =0.; VARMAX(YH2) =1.0
VARMIN(YN2) =0.; VARMAX(YN2) =1.0
VARMIN(FUEL)=1.E-9

    GROUP 21. Print-out of variables
OUTPUT(SHAD,Y,Y,Y,Y,Y,Y); OUTPUT(GAS ,Y,Y,Y,Y,Y,Y)
OUTPUT(MIXF,Y,Y,Y,Y,Y,Y); OUTPUT(FUEL,Y,Y,Y,Y,Y,Y)
OUTPUT(H1  ,Y,Y,Y,Y,Y,Y); OUTPUT(H2  ,Y,Y,Y,Y,Y,Y)
    GROUP 22. Spot-value print-out
IYMON=NY-2
TSTSWP=-1
    GROUP 23. Field print-out and plot control
ITABL=1; ORSIZ=0.2; Nplt=1

PATCH(PROFIL1,PROFIL,1,NX,1,NY,1,NZ,1,1)
COVAL(PROFIL1,YCO,0,0)
COVAL(PROFIL1,YCO2,0,0)
COVAL(PROFIL1,MIXF,0,0)
COVAL(PROFIL1,FUEL,0,0)

PATCH(PROFIL2,PROFIL,1,NX,1,1,NY,NZ,1,1)
COVAL(PROFIL2,P1,0,0)
COVAL(PROFIL2,GVEL,0,0)
COVAL(PROFIL2,5,0,0)
COVAL(PROFIL2,RHO1,0,0)

PATCH(PROFIL3,PROFIL,1,NX,1,NY,1,NZ,1,1)
COVAL(PROFIL3,TMP1,0,0)
COVAL(PROFIL3,TMP2,0,0)

   ========================================================
             6 gases model
   ========================================================

namsat=mosg

REAL(AIRO2,AIRN2,GASCON)
REAL(MN2,MC,MO2,MH2,MCO,MCO2,MH2O)
STORE(RMIX,HSUB,YN2,YH2,YO2,YCO,YCO2,YH2O,YSUM)
STORE(RHO1,RHO2,FLIM,FRAC,GO,GC,GH,GOFU,GOPA)

    Gas constant:
GASCON=8.3143e3

AIRO2=0.232;AIRN2=0.768
    Molecular masses:
MN2=28.; MC=12.; MO2=32.; MH2=2.; MCO=28.; MCO2=44.; MH2O=18.

RHO1=GRND
        DEN1=PRESS0/(RMIX*TMP1+tiny)
        DEN1=AMIN1(VARMAX(142),AMAX1(0.0,DEN1,VARMIN(142)))
CMDOT=GRND
        INTMDT=RCAR*VOL
CINT(H1)=GRND
        COI1(H1)=0.0
CINT(H2)=GRND
        COI2(H2)=RCAR*VOL
PHINT(H1)=GRND
        FII1(H1)=:CP2:*TMP2+:HCHX:
PHINT(H2)=GRND
        FII2(H2)=:CP2:*TMP1+:HCHX:

    Model settings
    ==============

    Sub-model 1:  Coal carbon oxidation
    -----------------------------------

   (1)  Gas mixture composition parameters

    FLIM=0.232/(0.232+:CINCL:*32./12.+$
                          :HINCL:*32./(2.*2.))
    GO=:AIRO2:*(1-MIXF)
    GC=:CINCL:*MIXF
    GH=:HINCL:*MIXF
    GOPA=GC*32./(2.*12.)/(1-GO+GC*32./(2*12.)+TINY)
    GOFU=(GH*32./(2.*2.)+GC*32./12.)/$
            (1.-GO+GH*32./(2.*2.)+GC*32./12.+TINY)
    FRAC=(GO-GOPA)/(GOFU-GOPA+TINY)

  (2) Mass fraction of nytrogen

    YN2=:NINCL:*MIXF+:AIRN2:*(1.-MIXF)

  (3) ** Region 1**  containing O2, CO2 & H2O

    YH2O=:HINCL:*MIXF*18./2.
   IF(MIXF.LE.FLIM)
    YCO2=:CINCL:*MIXF*44./12.
   IF(MIXF.LE.FLIM)
    YO2 =:AIRO2:*(1.-MIXF)-:CINCL:*MIXF*32./12.-$
                 :HINCL:*MIXF*32./(2.*2.)
   IF(MIXF.LE.FLIM)
    YCO=0.0
   IF(MIXF.LE.FLIM)
    YH2=0.0
   IF(MIXF.LE.FLIM)

  (4) ** Region 2 **  containing CO2, H2O, H2 & CO

    YH2O=:HINCL:*MIXF*18./2.*FRAC*(1-GOFU)/(1-GO+TINY)
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)
    YCO2=:CINCL:*MIXF*44./12.*FRAC*(1-GOFU)/(1-GO+TINY)
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)
    YO2=0.0
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)
    YCO=:CINCL:*MIXF*28./12.*(1-FRAC)*$
                 (1-GOPA)/(1-GO+TINY)
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)
    YH2=:HINCL:*MIXF*(1-FRAC)*(1-GOPA)/(1-GO+TINY)
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)

  (5) ** Region 3 ** containing H2 & CO.

    YH2O=0.0
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)
    YCO2=0.0
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)
    YO2=0.0
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)
    YCO=:AIRO2:*(1-MIXF)*2*28./32.
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)
    YH2=:HINCL:*MIXF
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)
       *************************************

    YH2=AMAX1(0.,YH2)
    YH2O=AMAX1(0.,YH2O)
    YCO=AMAX1(0.,YCO)
    YCO2=AMAX1(0.,YCO2)
    YN2=AMAX1(0.,YN2)
    YO2=AMAX1(0.,YO2)


    HSUB=0.0
   IF(MIXF.LE.FLIM)
    RMIX=:GASCON:*(YO2/32.+YH2O/18.+YCO2/44.+$
                 YN2/28.)
   IF(MIXF.LE.FLIM)

    HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)
    RMIX=:GASCON:*(YH2O/18.+YCO/28.+YCO2/44.+$
                 YH2/2.+YN2/28.)
   IF(MIXF.GT.FLIM.AND.FRAC.GE.0.)

    HSUB=YCO*:HCOCO2:+YH2*:HHH2O:
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)
    RMIX=:GASCON:*(YCO/28.+YH2/2.+YN2/28.)
   IF(MIXF.GT.FLIM.AND.FRAC.LT.0.)

    TMP1=(H1-HSUB)/1100.
    TMP2=(H2-:HCHX:)/1100.
    TMP1=AMIN1(VARMAX(140),AMAX1(100.,TMP1,VARMIN(140)))
    TMP2=AMIN1(VARMAX(139),AMAX1(100.,TMP2,VARMIN(139)))
    YSUM=YN2+YO2+YCO+YCO2+YH2O+YH2

    Interphase transport.
    --------------------
STORE(RCAR)

    Carbon mass transfer and related sources.
    -----------------------------------------

PATCH(car2gas,VOLUME,1,NX,1,NY,1,NZ,1,1)

  (1) Transfer of mass leading to increase of gas flow rate:

     VAL = RCAR
    COVAL(car2gas,P1,FIXFLU,GRND)

  (2) Transfer of carbon leading to increase of mixture
      fraction at the same rate:

     CO = RCAR
COVAL(car2gas,MIXF,GRND,1.)

  (3) Transfer of enthalpy and heat leading to increase of
      gas enthalpy at the same rate:

      VAL=:CP2:*TMP2 + :HCHX:
     COVAL(car2gas,H1,ONLYMS,GRND)

        *  Carbon from  coal patch
PATCH(carFcoal,VOLUME,1,NX,1,NY,1,NZ,1,1)

      1. Sink of second phase mass
    VAL= -RCAR
    COVAL(carFcoal,P2,FIXFLU,GRND)

      2. Heat From coal patch
PATCH(heatFcol,VOLUME,1,NX,1,NY,1,NZ,1,1)
   CO = RCAR
   VAL= TMP1*:CP2: + :HCHX:
  COVAL(heatFcol,H2,GRND,GRND)

   RCAR=:BURNRATE:*FUEL*(:FS:-MIXF)

nyprin=1;nzprin=1
dmpstk=t
 LIBREF=101
STOP