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.
The independent variables of the problem are the three components of cartesian coordinate system.
The main dependent (solved for) variables are:
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 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:
R_{e} = C_{e}HCOF(T_{e}-T_{1}) R_{wat}L^{-1}
where
T_{e} is a water saturation temperature;
R_{wat} stands for volume fraction of particle water in a reaction volume;
L is a latent heat of evaporation and
C_{e} 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:
R_{wat} = WAT2 MASS2 RHO1^{-1}/V_{cell}
wherein MASS2 is the mass of particulate phase, kg, and V_{cell} stands for cell volume, m^{3}.
The associated rate of water vapour generation is, obviously, equal to evaporation mass transfer:
R_{e, H2O} =R_{e}
The volatilisation of solid-fuel, raw coal, particles, is presumed to follow the one-step irreversible reaction type mechanism:
1 Kg Coal = Y_{V}Kg Volatiles + ( 1 - Y_{V} ) Kg Char
where Y_{V} 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, R_{coal}, in kg m^{-3}s^{-1}, is readily related to the rate of volatiles released through associated stoichiometric coefficient as follows:
R_{coal} = R_{V}/Y_{V}
where R_{V} is mass source of volatiles originating from the coal particles into the gas phase, kg m^{-3}s^{-1}.
The rate of char formation, R_{char}, in kg m^{-3}s^{-1}, is also calculated through the volatilisation rate as:
R_{char} = ( 1 - Y_{V} )R_{V}/Y_{V}
Volatilisation kinetics are coal-rank dependent. The volatilisation rate is modelled through an Arrhenius-type expression:
R_{V} = 2.10^{3}e^{-2829/T2}m_{V2}M_{VPV}
Here
m_{V2} = Y_{V}COL2
The mass of voilatiles per unit of reaction volume is given as:
M_{VPV} = RHO2 Y_{V}R2
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: CH_{4} + 1.5 ( O_{2} + 3.76N_{2}) = CO + 2H_{2}O + 5.64N_{2}
Step 2: CO + 0.5( O_{2} + 3.76N_{2} ) = CO_{2} + 1.88N_{2}
The reaction rates of combustion are obtained as the harmonic blend of a Arrhenius kinetics and eddy-dissipation rates:
R^{com}_{CH4} = - [ ( R^{k}_{CH4})^{-1} + (R^{e}_{CH4} )^{-1} ]^{-1} and
R^{com}_{CO} = - [ ( R^{k}_{CO} )^{-1} + ( R^{e}_{CO} )^{-1} ]^{-1}
where R^{k} and R^{e}, in kg/m^{3}s, are the kinetic and eddy-dissipation rates:
R^{e}_{CH4} = 4 RHO1 EP/KE min( YCH4, YO2/3)
R^{e}_{CO} = 4 RHO1 EP/KE min( YCO, YO2/0.57)
R^{k}_{CH4} =
1.15 10^{9}RHO1^{2}e^{-24444/T}
YCH4^{-0.3}YO2^{1.3}
R^{k}_{CO} =
5.42 10^{9}RHO1^{2}e^{-15152/T}
YO2^{0.25}YH2O^{0.5}YCO
The remaining rates are defined through associated stoichiometric coefficients:
Step 1:
R^{1}_{O2} =
3 R^{com}_{CH4}
R^{1}_{CO} =
-1.75 R^{com}_{CH4}
R^{1}_{H2O} =
-2.25 R^{com}_{CH4}
Step 2:
R^{2}_{O2} =
0.57 R^{com}_{CO}
R^{2}_{CO2} =
-1.57 R^{com}_{CO}
The net rates of species generation in volatile combustion are the balances of formation and combustion as appropriate:
R_{V,CH4} = R^{com}_{CH4}
R_{V,CO} = R^{com}_{CO} +
R^{1}_{CO}
R_{V,O2} = R^{1}_{O2} +
R^{2}_{O2}
R_{V,CO2} = R^{2}_{CO2}
R_{V,H2O} = R^{1}_{H2O}
The particle char burns creating carbon oxides according to the reaction:
C + 0.5(1 + w) O_{2} = ( 1 - w)CO + wCO_{2}
where w is the mole fraction of carbon oxides formed as CO_{2} which is related to the mass fractions as follows:
w = M_{CO}YCO2/( M_{CO}YCO2 + M_{CO2}YCO )
where M_{CO} = 28 and M_{CO2} =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, R_{C}, in kg m^{-3}s^{-1}, is obtained as an harmonic blend of a kinetic- and diffusion-controlled mass-transfer rates:
R_{C} = C_{c}A_{s} [ ( K^{k}_{C})^{-1} + (K^{d}_{C} )^{-1} ]^{-1}P_{O2} R_{CHA}/R2
where
A_{s} = 6R2/SIZE, is volumetric particle surface area, 1/m;
P_{O2} is partial pressure of O_{2}, N/m^{2};
C_{c} is an empirical constant accounting for the difference between particle surface area and effective char one exposed to the gas phase and
R_{CHA} is a volume fraction of a char in the reaction volume.
The latter is calculated as:
R_{CHA} = CHA2 MASS2 RHO2^{-1}/ V_{cell}
The partial pressure P_{O2} is calculated using Dalton's law:
P_{O2} = x_{O2}P1
Wherein, the mole fraction of O_{2}, x_{O2}, is related to mass fraction, YO2, through molecular masses of O_{2}, M_{O2}, and the mixture, M_{mix}, as
x_{O2} = M_{mix}/M_{O2}YO2
Kinetic mass transfer coefficient, K^{k}_{C}, in s/m, is calculated as:
K^{k}_{C} = 0.1309 e^{-26850/T2}
Diffusion mass transfer coefficient, K^{d}_{C}, in s/m, is calculated as follows:
K^{d}_{C} = Sh D_{O2}M_{C}/ ( R_{gas}T_{2}SIZE )
where
R_{C,O2} = - 0.5(1 + w) M_{O2}/M_{C}
The rate of CO produced by char combustion, in kg m^{-3}s^{-1}, is related to R_{C} through its own stoichiometric coefficient:
R_{C,CO} = (1 - w) M_{CO}/M_{C}
The rate of CO_{2} produced by char combustion, in kg m^{-3}s^{-1}, is related to R_{C} through the stoichiometry as follows:
R_{C,CO2} = wM_{CO2}/M_{C}
Carbon monoxide emerging from char further burns to carbon dioxide in the gas phase as a participant of the second step of volatile combustion.
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, NO_{2}, and nitrous oxide, N_{2}O.
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.
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/m^{3}/s, is given by:
R_{T,NO} = 2 RHO1 K_{1} YN2 [O] M_{NO}/ M_{N2}
where M_{NO} = 30, is NO molecular mass,
M_{N2} =28, is N_{2} molecular mass, and,
K_{1} = 1.8 10^{8}e^{-38370/T1},
m^{3}mol^{-1}s^{-1},
stands for the reaction-rate constant of the forward
reaction:
N_{2} + O --> NO +H
The equilibrium O-atom concentration, [O], can be obtained from the expression:
[O] = 3.97 10^{5}T1^{-0.5}(RHO1 YO2/M_{O2}) ^{0.5}e^{-31090/T1}, mol/m^{3}
wherein M_{O2} =32, is the molecular mass of oxygen.
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:
It is assumed that all char-bound nitrogen converts to HCN. Thus,
R_{C,HCN} = R_{C}m_{NC}M_{HCN}/M_{N}
where
R_{C} is char burnout rate, kg m^{-3}s^{-1} ;
m_{NC} is a mass fraction of nitrogen in char
( specified input parameter ) ;
M_{HCN} = 27, and M_{N} = 14, are the molecular
masses of HCN and N, respectively.
The source of HCN from volatiles is related to the rate of volatile release:
R_{V,HCN} = R_{V}m_{NV}M_{HCN}/M_{N}
where
R_{V} is the mass source of volatiles,
kg m^{-3}s^{-1} ;
m_{NV} is a mass fraction of nitrogen in the volatiles
( specified input parameter ).
HCN is oxidised to NO via reaction:
4HCN + 5O_{2} = 4NO +2H_{2}O + 4CO
The reaction rate, in kg m^{-3}s^{-1}, is given as:
R_{NO+} = RHO1 10^{11}e^{-33678/T1}YO2^{a} YHCN
where the oxygen reaction order, a, is assumed to be unity.
The homogeneous reaction of NO reduction is modelled as:
4HCN + 6NO = 5N_{2} + 2H_{2}O + 4CO
The reaction rate, in kg m^{-3}s^{-1}, is calculated as
R_{NO-1} = RHO1 3.10^{12}e^{-30069/T1}YHCN YNO
The heterogeneous reaction rate, in kg m^{-3}s^{-1}, of NO reduction on the char surface is modelled as:
R_{NO-2} = 4.10^{-4}e^{-18042/T2} A_{s}P_{NO}
where P_{NO} is NO partial pressure, N/m^{2}.
The partial pressure P_{NO} is calculated using Dalton's law:
P_{NO} = x_{NO}P1
Wherein, the mole fraction of NO, x_{NO}, is related to mass fraction, YNO, through molecular masses of NO, M_{NO}, and the mixture, M_{mix}, as
x_{NO} = M_{mix}/M_{NO}YNO
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, R_{I} in W m^{-2}, accounting for radiating particles and gases together can be derived as:
div ( G_{rad}gradR_{I} ) + a_{g} ( 4sT_{1}^{4} - R_{I}) + a_{p} ( 4sT_{2}^{4} - R_{I}) = 0
where s = 5.68 10^{-8}, is Steffan-Boltzman constant, W m^{-2} K^{-4}, T_{1} and T_{2} are the gas and particle temperatures, K.
The exchange coefficient, G_{rad}, is expressed by:
G_{rad} = ( 3(a_{g}+s_{g}) + 3(a_{p}+s_{p}) - C_{g}s_{g} )^{-1}
where
a_{p} = e_{p}A_{p}
where e_{p} is the emmisivity of particle and A_{p} is the volumetric particle projected area, m^{-1}.
The latter is calculated from particulate volume fraction, r_{2} and current particle diameter, d_{p} as follows:
A_{p} = 1.5r_{2}d_{p}^{-1}
The equivalent particle scattering coefficient is defined as:
s_{p} = (1 - s_{p}) (1 - e_{p}) A_{p}
where s_{p} stands for particle scattering factor.
The s_{p} and e_{p} are related by the "incident-radiation-sharing" equation:
s_{p} + e_{p} = 1
The symmetry factor, C_{g}, is used to model anisotroping scattering by means of a linear-anisotropic scattering phase function. C_{g} 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 C_{g}=1 corresponding to complete forward scattering.
A negative value means that more radiant energy is scattered backward than forward with C_{g}=-1 standing for complete backward scattering.
A zero value of C_{g} 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:
S_{H1,rad} = a_{g} (R_{I} - 4sT_{1}^{4}) + a_{p} (R_{I} - 4sT_{2}^{4})
The volumetric heat source, due to particle radiation, included in the particulate phase energy equations is as follows:
S_{H2,rad} = e_{p}A_{s} (0.25R_{I} - sT_{2}^{4})
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:
S_{R, wall} = 0.5 e_{w} (4sT_{w}^{4} - R_{I})(2 - e_{w})^{-1}
where e_{w} 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.
Only the compilation of the model sources which required non-standard settings is provided below.
Interphase mass transfer, kg/s
CMDOT = V_{cell}( R_{V} + R_{C} + R_{e})
2nd phase composition, kg m^{-3}s^{-1}
1st phase composition, kg m^{-3}s^{-1}
Incident radiation, W m^{-3}
The source term for the transport equation of R_{I} is:
S_{R} = a_{g} ( 4sT_{1}^{4} - R_{I}) + a_{p} ( 4sT_{2}^{4} - R_{I})
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:
S_{H1} = S_{H1,com} + S_{H1,rad } + S_{H1,int}
The contribution from the heats of reactions is:
S_{H1,com} = R^{com}_{CH4}H°_{CH4} + R^{com}_{CO}H°_{CO} + R_{C}H°_{C}
The interphase heat transfer contribution is:
S_{H1,int} = HCOF( C_{S}(T_{2} - T_{1}) R_{sol} + C_{e}(T_{e} - T_{1}) R_{wat} )
where
C_{S} is an empirical constant accounting for the difference between particle surface area and effective solid one exposed to the gas phase and
R_{sol} is the volume fraction of solid components of the particulate phase.
The latter is given by:
R_{sol} = R2 - R_{wat}
2nd phase energy, W m^{-3}
The source term for 2nd phase specific enthalpy, H2, is caculated as:
S_{H2} = S_{H2,rad } + S_{H2,int}
where
S_{H2,int} = C_{S}HCOF( T_{1} - T_{2} ) R_{sol} - R_{e}L
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
The specific enthalpies are related to the phase temperatures and specific heats:
H1 = C_{P,1}T_{1} and
H2 = C_{P,2}T_{2}
The solid-phase density is taken as constant, as the ideal gas law is used for the gas phase as follows:
RHO1 = M_{mix}R^{-1}_{gas}T_{1}^{-1}P1
where the molecular mass of the mixture is composition-dependent:
M^{-1}_{mix} = 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, C_{P,1}, is assumed to be equal for all gas components and is a function of gas temperature as follows:
C_{P,1} = 1059 + 0.25( T_{1} - 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:
H°_{C} = wH°_{C2CO2} + (1 - w)H°_{C2CO}
where H°_{C2CO2} and H°_{C2CO} are the heats of carbon oxidation to CO_{2} and CO, J/kgC, respectively.
Interphase particle heat transfer is represented by a sphere-law correlation:
NUSS = 2. +0.65 REYN^{0.5}
The heat transfer coefficient, HCOF, per unit of particulate volume, is then evaluated as:
HCOF = A_{s}k_{gas}SIZE^{-1}NUSS/R2
where, k_{gas} 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.
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 :