Fuel-dust flames


BY : Dr S V Zhubrin, CHAM Ltd

DATE : March 2002

FOR : Multi-physics case for V3.4


The case presents the implementation in PHOENICS of an IPSA based model for non-equilibrium two-phase reactive flow of combustible particles finely dispersed in the carrying air stream. The model uses the Eulerian description of the two interpenetrating continua with the transfer of heat, mass and momentum between them. The processes accounted for are particle drying, the volatilisation of dispersed phase material, gaseous exothermic oxidation of volatile products, heterogeneuos combustion of particle char and NOx formation.

The particle drying is represented by an assumption that the heat transferred to the particle is used both for heating-up its solid components and as a latent heat of phase change for water evaporation.

The devolatilisation of particulate phase is considered as kinetically driven.

Turbulent combustion of volatiles is modelled via two-step reaction of hydrocarbon oxidation, in which carbon monoxide is an intermediate product.

The char combustion is represented by blended mechanism of the oxygen diffusion to the particle and chemical kinetic. The relative amount of char combustion products are assumed to be dependent on the particle temperature.

The NOx formation is represented by several submodels, such as oxidation of nitrogen present in the combustion air and that contained in the fuel. The reduction of nitric oxides both in gas phase and on char surface is also taken into account.

The turbulence is accounted for by conventional K-e model.

The radiation modelling and its contribution to particle heat-up is based on P-1 approximation extended and suitably generalized to handle the effects of absorbing, emmiting, and scattering particles on anisotropic scattering in absorbing, emmiting, and scattering media.

The model is applied to the pulverized coal combustion in a wall-fired model furnace.


Conservation equations

The independent variables of the problem are the three components of cartesian coordinate system.

The main dependent (solved for) variables are:

The main auxiliary variables are:

Main features of the model

For the second, particulate, phase, raw coal, COL2, is consumed during devolatilisation, and char, CHA2, solid, combustible residue of volatilisation, appears as COL2 disappears. The char then disappears as it burnt in the char combustion. The rise of the particle temperature also causes the coal to dry out and transfer its water contents, WAT2, to the gas, as water vapour. The ash contents of the phase, ASH2, increase to value 1 in completely burned dry particles.

For the first, gas, phase volatile fuel, YCH4, appears as a consequence of volatilisation, and then disappears as it burned in gas-phase combustion. The oxygen contents, YO2, decrease as it is consumed in both volatile and char combustion.

Particle drying

Particle drying is accounted for by assuming that the heat transferred from the gas phase is used both for heating-up the solid components of the particle and as a latent heat of phase change for the evaporation of the water.

The associated evaporation mass transfer, in kg s-1 m-3, is then given by:

Re = CeHCOF(Te-T1) RwatL-1


Te is a water saturation temperature;
Rwat stands for volume fraction of particle water in a reaction volume;
L is a latent heat of evaporation and
Ce is an empirical constant accounting for the difference between particle surface area and effective water one exposed to the gas phase.

The volume fraction of particle water is calculated as:

Rwat = WAT2 MASS2 RHO1-1/Vcell

wherein MASS2 is the mass of particulate phase, kg, and Vcell stands for cell volume, m3.

The associated rate of water vapour generation is, obviously, equal to evaporation mass transfer:

Re, H2O =Re

Volatilisation mechanism

The volatilisation of solid-fuel, raw coal, particles, is presumed to follow the one-step irreversible reaction type mechanism:

1 Kg Coal = YVKg Volatiles + ( 1 - YV ) Kg Char

where YV stands for the mass fraction of volatiles in coal which is coal-rank dependent and is, normally, the user specified input.

As a consequence, the rate of coal consumption, Rcoal, in kg m-3s-1, is readily related to the rate of volatiles released through associated stoichiometric coefficient as follows:

Rcoal = RV/YV

where RV is mass source of volatiles originating from the coal particles into the gas phase, kg m-3s-1.

The rate of char formation, Rchar, in kg m-3s-1, is also calculated through the volatilisation rate as:

Rchar = ( 1 - YV )RV/YV

Devolatilisation kinetics

Volatilisation kinetics are coal-rank dependent. The volatilisation rate is modelled through an Arrhenius-type expression:

RV = 2.103e-2829/T2mV2MVPV


The mass fraction of volatiles remaining in a second, particulate, phase is related to the mass fraction of raw coal:

mV2 = YVCOL2

The mass of voilatiles per unit of reaction volume is given as:


Volatiles combustion

The volatiles are taken as 100% methane. Combustion of the volatiles is treated as a two-step irreversible chemical reaction of methane oxidation as follows:

Step 1: CH4 + 1.5 ( O2 + 3.76N2) = CO + 2H2O + 5.64N2

Step 2: CO + 0.5( O2 + 3.76N2 ) = CO2 + 1.88N2

The reaction rates of combustion are obtained as the harmonic blend of a Arrhenius kinetics and eddy-dissipation rates:

RcomCH4 = - [ ( RkCH4)-1 + (ReCH4 )-1 ]-1 and

RcomCO = - [ ( RkCO )-1 + ( ReCO )-1 ]-1

where Rk and Re, in kg/m3s, are the kinetic and eddy-dissipation rates:

ReCH4 = 4 RHO1 EP/KE min( YCH4, YO2/3)
ReCO = 4 RHO1 EP/KE min( YCO, YO2/0.57)
RkCH4 = 1.15 109RHO12e-24444/T YCH4-0.3YO21.3
RkCO = 5.42 109RHO12e-15152/T YO20.25YH2O0.5YCO

The remaining rates are defined through associated stoichiometric coefficients:

Step 1:
R1O2 = 3 RcomCH4
R1CO = -1.75 RcomCH4
R1H2O = -2.25 RcomCH4

Step 2:
R2O2 = 0.57 RcomCO
R2CO2 = -1.57 RcomCO

The net rates of species generation in volatile combustion are the balances of formation and combustion as appropriate:

RV,CH4 = RcomCH4
RV,CO = RcomCO + R1CO
RV,O2 = R1O2 + R2O2
RV,CO2 = R2CO2
RV,H2O = R1H2O

Combustion of char

The particle char burns creating carbon oxides according to the reaction:

C + 0.5(1 + w) O2 = ( 1 - w)CO + wCO2

where w is the mole fraction of carbon oxides formed as CO2 which is related to the mass fractions as follows:


where MCO = 28 and MCO2 =44, are the molecular masses of carbon monoxide and carbon dioxide, respectively.

The mass split between carbon oxides is assumed to depend on the particle temperature as:

YCO/YCO2 = 2500e-6249/T2

The char burnout rate, RC, in kg m-3s-1, is obtained as an harmonic blend of a kinetic- and diffusion-controlled mass-transfer rates:

RC = CcAs [ ( KkC)-1 + (KdC )-1 ]-1PO2 RCHA/R2


As = 6R2/SIZE, is volumetric particle surface area, 1/m;
PO2 is partial pressure of O2, N/m2;
Cc is an empirical constant accounting for the difference between particle surface area and effective char one exposed to the gas phase and
RCHA is a volume fraction of a char in the reaction volume.

The latter is calculated as:

RCHA = CHA2 MASS2 RHO2-1/ Vcell

The partial pressure PO2 is calculated using Dalton's law:

PO2 = xO2P1

Wherein, the mole fraction of O2, xO2, is related to mass fraction, YO2, through molecular masses of O2, MO2, and the mixture, Mmix, as

xO2 = Mmix/MO2YO2

Kinetic mass transfer coefficient, KkC, in s/m, is calculated as:

KkC = 0.1309 e-26850/T2

Diffusion mass transfer coefficient, KdC, in s/m, is calculated as follows:

KdC = Sh DO2MC/ ( RgasT2SIZE )


The rate of O2 consumed by char combustion, in kg m-3s-1, is related to RC through stoichiometric coefficient:

RC,O2 = - 0.5(1 + w) MO2/MC

The rate of CO produced by char combustion, in kg m-3s-1, is related to RC through its own stoichiometric coefficient:

RC,CO = (1 - w) MCO/MC

The rate of CO2 produced by char combustion, in kg m-3s-1, is related to RC through the stoichiometry as follows:


Carbon monoxide emerging from char further burns to carbon dioxide in the gas phase as a participant of the second step of volatile combustion.

The formation of NOx

NOx is collectively referred to the number of nitrogen oxidation products present in the gases emerging from combustion. They are considered as hazardous pollutants.

NOx emmision consists of mostly nitric oxide, NO. Less significant are nitrogen oxide, NO2, and nitrous oxide, N2O.

Because the species involved are normally present in the low concentarion trace quantities, the NOx modeling is assumed not to contribute neither to the gas mixture properties nor to the mass conservation of mixture components.

The formation of NOx is attributed to thermal NOx formation and fuel NOX formation with subsequent reduction.

Thermal NOx

The rate of the oxidation of atmospheric nitrogen present in the combustion air is modelled here via the steady-state simplification of Zeldovich mechanism with partial-equilibrium assumptions.

The thermal NOx formation rate, in kg/m3/s, is given by:

RT,NO = 2 RHO1 K1 YN2 [O] MNO/ MN2

where MNO = 30, is NO molecular mass,
MN2 =28, is N2 molecular mass, and,
K1 = 1.8 108e-38370/T1, m3mol-1s-1,
stands for the reaction-rate constant of the forward reaction:

N2 + O --> NO +H

The equilibrium O-atom concentration, [O], can be obtained from the expression:

[O] = 3.97 105T1-0.5(RHO1 YO2/MO2) 0.5e-31090/T1, mol/m3

wherein MO2 =32, is the molecular mass of oxygen.

Fuel NOx

The extent of conversion of fuel nitrogen to NOx is dependent on the local combustion characteristics and the initial concentration of nitrogen-bound compounds. Fuel-bound nitrogen-containing compounds are released into the gas phase when the fuel particles are heated and devolatilized. From the thermal decomposition of these compounds the intermediates are formed. For simplicity, the nitrogen-containing intermediates are here grouped to be hydrogen cyanide, HCN, only.

The fuel-NOx mechanism employed also assumes that fuel nitrogen is distributed between volatiles and char.

The source terms for NOx and HCN are determined as follows:

  1. HCN from char

    It is assumed that all char-bound nitrogen converts to HCN. Thus,



    RC is char burnout rate, kg m-3s-1 ;
    mNC is a mass fraction of nitrogen in char ( specified input parameter ) ;
    MHCN = 27, and MN = 14, are the molecular masses of HCN and N, respectively.

  2. HCN form volatiles

    The source of HCN from volatiles is related to the rate of volatile release:



    RV is the mass source of volatiles, kg m-3s-1 ;
    mNV is a mass fraction of nitrogen in the volatiles ( specified input parameter ).

  3. NOx formation

    HCN is oxidised to NO via reaction:

    4HCN + 5O2 = 4NO +2H2O + 4CO

    The reaction rate, in kg m-3s-1, is given as:

    RNO+ = RHO1 1011e-33678/T1YO2a YHCN

    where the oxygen reaction order, a, is assumed to be unity.

  4. NOx reduction in a gas phase

    The homogeneous reaction of NO reduction is modelled as:

    4HCN + 6NO = 5N2 + 2H2O + 4CO

    The reaction rate, in kg m-3s-1, is calculated as

    RNO-1 = RHO1 3.1012e-30069/T1YHCN YNO

  5. NOx reduction on char surface

    The heterogeneous reaction rate, in kg m-3s-1, of NO reduction on the char surface is modelled as:

    RNO-2 = 4.10-4e-18042/T2 AsPNO

    where PNO is NO partial pressure, N/m2.

    The partial pressure PNO is calculated using Dalton's law:

    PNO = xNOP1

    Wherein, the mole fraction of NO, xNO, is related to mass fraction, YNO, through molecular masses of NO, MNO, and the mixture, Mmix, as

    xNO = Mmix/MNOYNO

Thermal radiation

Thermal radiation is modelled by the expanding the radiation intensity in terms of first order spherical harmonics.

Model equations

Assuming that only four terms representing the moments of the intensity are used, the conservation equation of incident radiation, RI in W m-2, accounting for radiating particles and gases together can be derived as:

div ( GradgradRI ) + ag ( 4sT14 - RI) + ap ( 4sT24 - RI) = 0

where s = 5.68 10-8, is Steffan-Boltzman constant, W m-2 K-4, T1 and T2 are the gas and particle temperatures, K.

The exchange coefficient, Grad, is expressed by:

Grad = ( 3(ag+sg) + 3(ap+sp) - Cgsg )-1


The equivalent particle absorption coefficient is defined as:

ap = epAp

where ep is the emmisivity of particle and Ap is the volumetric particle projected area, m-1.

The latter is calculated from particulate volume fraction, r2 and current particle diameter, dp as follows:

Ap = 1.5r2dp-1

The equivalent particle scattering coefficient is defined as:

sp = (1 - sp) (1 - ep) Ap

where sp stands for particle scattering factor.

The sp and ep are related by the "incident-radiation-sharing" equation:

sp + ep = 1

The symmetry factor, Cg, is used to model anisotroping scattering by means of a linear-anisotropic scattering phase function. Cg ranges from -1 to +1 and represents the amount of radiation scattered in forward direction.

A positive value indicates that more radiant energy is scattered forward than backward with Cg=1 corresponding to complete forward scattering.

A negative value means that more radiant energy is scattered backward than forward with Cg=-1 standing for complete backward scattering.

A zero value of Cg defines the scattering that is equally likely in all directions, i.e. isotropic scattering

Phase-energy source terms

The volumetric source term, in W m-3, for the gas mixture enthalpy, due to radiation, is given by:

SH1,rad = ag (RI - 4sT14) + ap (RI - 4sT24)

The volumetric heat source, due to particle radiation, included in the particulate phase energy equations is as follows:

SH2,rad = epAs (0.25RI - sT24)

Boundary conditions

For symmetry planes and perfectly-reflecting boundaries, the radiation boundary conditions are assumed to be zero-flux type.

For the incident radiation equation, the following boundary sources per unit area are used at the walls:

SR, wall = 0.5 ew (4sTw4 - RI)(2 - ew)-1

where ew is the wall emmisivity.

The sources for incident radiation at the inlets and outlets are computed in the manner similar to the walls.

Often, it is safe to assume that the emmisivity of all flow inlets and outlets is unity (black body absorption). If the temperature outside the inlet or outlet considerably differs from that in the enclosure, the different temperatures should be used for radiation and convection fluxes at inlet and outlets.


Model sources

Only the compilation of the model sources which required non-standard settings is provided below.

Interphase mass transfer, kg/s

CMDOT = Vcell( RV + RC + Re)

2nd phase composition, kg m-3s-1

1st phase composition, kg m-3s-1

Incident radiation, W m-3

The source term for the transport equation of RI is:

SR = ag ( 4sT14 - RI) + ap ( 4sT24 - RI)

The symmetry plane, the perfectly-reflecting walls and inlets are assumed to be zero-source regions.

1st phase energy, W m-3

The reaction heats, i.e. the heats of combustions for volatile methane, H°CH4, carbon monoxide, H°CO, and char, H°C, are released with the rates of their consumptions.

The heat is also transfered due to thermal radiation and from the 2nd phase via interphase heat transfer.

In terms of the transport equation for H1, the resulting source term is:

SH1 = SH1,com + SH1,rad + SH1,int

The contribution from the heats of reactions is:

SH1,com = RcomCH4CH4 + RcomCOCO + RCC

The interphase heat transfer contribution is:

SH1,int = HCOF( CS(T2 - T1) Rsol + Ce(Te - T1) Rwat )


CS is an empirical constant accounting for the difference between particle surface area and effective solid one exposed to the gas phase and
Rsol is the volume fraction of solid components of the particulate phase.

The latter is given by:

Rsol = R2 - Rwat

2nd phase energy, W m-3

The source term for 2nd phase specific enthalpy, H2, is caculated as:

SH2 = SH2,rad + SH2,int


SH2,int = CSHCOF( T1 - T2 ) Rsol - ReL

Inert components

The mass fraction of inert ash is related to the other particulate phase components by "phase-sharing" equation:

ASH2 = 1 - COL2 - CHA2 - WAT2

The mass fraction of nitrogen in a gas phase is similarly calculated as:

YN2 = 1 - YO2 - YCO - YCO2 - YCH4 - YH2O

Temperature calculations

The specific enthalpies are related to the phase temperatures and specific heats:

H1 = CP,1T1 and
H2 = CP,2T2

Densities, specific- and reaction-heats

The solid-phase density is taken as constant, as the ideal gas law is used for the gas phase as follows:

RHO1 = MmixR-1gasT1-1P1

where the molecular mass of the mixture is composition-dependent:

M-1mix = YCH4/16+YO2/32+YN2/28+YCO/18+YCO2/44+YH2O/18

The solid-phase specific heats are taken as constant, and the gas-phase specific heat, CP,1, is assumed to be equal for all gas components and is a function of gas temperature as follows:

CP,1 = 1059 + 0.25( T1 - 300 )

The combustion heats of volatiles and carbon monoxide are taken as constants. The heat of combustion for char is assumed to depend on the product split ratio:

C = wC2CO2 + (1 - w)H°C2CO

where H°C2CO2 and H°C2CO are the heats of carbon oxidation to CO2 and CO, J/kgC, respectively.

Interphase heat transfer

Interphase particle heat transfer is represented by a sphere-law correlation:

NUSS = 2. +0.65 REYN0.5

The heat transfer coefficient, HCOF, per unit of particulate volume, is then evaluated as:


where, kgas is a gas heat conductivity taken as constant.



The furnace chamber chosen for the present computations, idealised though, retains many features of real-life industrial equipment: twenty four coal-fired burners are located at the lower part on the furnace front wall; air-pulverised coal particles are injected axially into the combustion chamber; they are ignited; steady state combustion takes place; combustion products flow towards the outlet at the upper part of the chamber.

The furnace is symmetrical, and only half of it is simulated and represented.

Boundary conditions

At all inlets, values are given of all dependent variables together with the prescribed flow rates.

At the outlet, the fixed exit pressure is set equal to zero and the computed pressures are relative to this pressure.

At the walls, standard "wall-functions" are used for the gas velocities, and the condition of zero flux is assumed for all other dependent variables.


All model settings have been made by PIL commands and PLANT settings of PHOENICS 3.4

The relevant Q1 file can be downloaded by clicking here.


The plots show the flow distribution, mixture composition as represented by the model, gas temperature and velocity within the furnace.

Pictures are as follows :