A Multi-Fluid Model is presented for turbulent diffusion flames.
The single discretised population-distinguishing attribute is the 'mixture fraction, i.e. the mass of material from the fuel-bearing stream which is present in unit mass of the local mixture.
Each thus-distinguished fluid is allowed also to possess several continuously-varying attributes, for each of which a balance equation is solved.
The model takes account of finite-rate chemical reactions for both adiabatic and non-adiabatic systems including thermal radiation.
It is applicable to "diffusion-flame" combustion proceeding through multi-step mechanisms of homo-and heterogeneous reactions followed by chemistry-rate-limited NOx formation.
The computational nature of the model is described in detail and is then applied to:
The model employed is of the general kind described in the PHOENICS documentation as a multi-fluid model.
It thus exhibits the following general features:
Summary of the particular assumptions
The PDA is uniformly discretized; therefore each fluid of the population has its own average value of mixture fraction uniquely pre-defined as follows:
f_{mix,k} = ( k-1 )/( N_{fluid} - 1 )
where k is a current fluid number and N_{fluid} stands for the total number of fluids in the population.
The implication of the above in the present context is that the lowest-numbered fluid, k=1, is always fuel-free, f_{mix,1} = 0, while the largest-number one represents the pure fuel, f_{mix,N} =1.
The choice of CVAs, Continuously-Varying Attributes, depends upon the situation in question. It forms the modeller's options and will be specified in the sections corresponding to each of the cases considered.
A few guidelines as mentioned in [ 1 ] are as follows:
Unlike PDA, the values of CVA are not pre-defined. They are calculated from their own conservation equations as described in the later section.
The mass fraction of each fluid, or its "presence probability", m_{k}, in a multi-fluid population, is assumed to be a conserved quantity. Its value at each point in the flow domain is computed by PHOENICS through solution of the following conservation equations of conventional type:
d(rm_{k})/dt+ div(rVm_{k}- G_{t} grad m_{k} ) = R_{m,k}
R_{m,k}, the net rate of k-fluid generation, is the balance of micromixing rate, R_{mix,k}, and interphase transfer, S_{p,k}:
R_{m,k} = R_{mix,k} + S_{p,k}
The source term, S_{p,k}, is due solely to transfer of mass into the gas phase from reacting particles (e.g. coal). In all other cases there are no such a source.
The term, R_{mix,k}, results from micromixing of the fluids as they move past, or collide with, each other in their turbulent motion. It is expressed, for uniformly-divided population, as:
R_{mix,k} = r S_{i} S_{j} F_{k,i,j} m_{i}m_{j} T_{i,j}
wherein:
T_{i,j} = C_{mix}e/K
with K standing for the kinetic energy of turbulence, e for its dissipation rate and C_{mix} for an empirical constant.
The fractional loss of mass is computed by following rules:
F_{k,i,j} = -0.5 for k=i or k=j and j greater than i+1, = 0.0 for k less than i or k greater than j or j=i+1, = 1/(j-i-1) for all other values of i, j and k.
These hypotheses are what is called, in MFM parlance, the Promiscuous-Mendelian coupling/splitting scheme[ 2 ].
The sources that result from the above scheme, applied to the interactions between,say, 5 fluids with T=r=1, would be in fact as follows:
R_{mix,1} = -0.5(m_{3}+
m_{4}+
m_{5})m_{1}
R_{mix,2} = m_{1}m_{3}+
m_{1}m_{4}/2+
m_{1}m_{5}/3-
0.5(m_{4}+
m_{5})m_{2}
R_{mix,3} = m_{4}(m_{1}/2+m_{2})+
m_{5}(m_{1}/3+m_{2}/2)-
0.5(m_{1}+
m_{5})m_{3}
R_{mix,4} =
m_{3}m_{5}+
m_{2}m_{5}/2+
m_{1}m_{5}/3-
0.5(m_{1}+
m_{2})m_{4}
R_{mix,5} = -0.5(m_{1}+
m_{2}+
m_{3})m_{5}
Let C_{k} be the value of a continuously-varying attribute of fluid k. The conservation equation for C_{k} takes the following general form:
d(rC_{k})/dt+ div(rVC_{k}- G_{t} grad C_{k} ) = Rc_{k} + Sc_{k,p}
where Rc_{k} is within-fluid mass rate of creation and depletion of C_{k} and Sc_{k,p} is the rate of creation by addition from the dispersed phase, if any.
The net rate of within-fluid generation is given by the balance of the sources resulting from the ij encounters, Rc_{mix,k}, and the source of CVA due to its in-fluid generation and/or dissipation, Rc_{gen,k}:
Rc_{k} = Rc_{mix,k} + Rc_{gen,k}
The contributions resulting from micro-mixing by fluid encounters are written as:
Rc_{mix,k}= S_{i}M_{k,i}(C_{i}-C_{k})
where M_{k,i} is the micro-mixing mass transfer which enters the fluid k from fluid i.
It is calculated as the i-related portion of the total mass transfer to each fluid in ij encounters:
S_{i}M_{k,i}= rS_{i} S_{j} Fc_{k,i,j} m_{i}m_{j}T_{i,j}
where
Fc_{k,i,j} = 0.0 for k less or equal than i or k greater or
equal than j or j=i+1,
= 1/(j-i-1) for all other values of i, j and k.
For the five-fluids population with T=r=1, the resulting sources are, in fact, as follows:
Rc_{mix,1} = 0 ;
Rc_{mix,2} =
(m_{1}m_{3}/2+m_{1}m_{4}/4+m_{1}m_{5}/6)
(C_{1}-C_{2})
+m_{1}m_{3}/2
(C_{3}-C_{2})
+m_{1}m_{4}/4
(C_{4}-C_{2})
+m_{1}m_{5}/6
(C_{5}-C_{2}) ;
Rc_{mix,3} =
(m_{1}m_{4}/4+m_{1}m_{5}/6)
(C_{1}-C_{3})
+(m_{2}m_{4}/2+m_{2}m_{5}/4)
(C_{2}-C_{3})
+(m_{1}m_{4}/4+m_{2}m_{4}/2)
(C_{4}-C_{3})
+(m_{1}m_{5}/6+m_{2}m_{5}/4)
(C_{5}-C_{3}) ;
Rc_{mix,4} =
m_{1}m_{5}/6
(C_{1}-C_{4})
+m_{2}m_{5}/4
(C_{2}-C_{4})
+m_{3}m_{5}/2
(C_{3}-C_{4})
+(m_{3}m_{5}/2+m_{2}m_{5}/4+m_{1}m_{5}/6)
(C_{5}-C_{4}) ;
Rc_{mix,5} = 0
The generation/dissipation rates, Rc_{gen,k}, that appear as source terms of CVA are usually problem-specific.
For in-fluid chemical reaction, they can be computed from Arrhenius rate expressions, using the eddy-dissipation concept or blending of two, as appropriate.
For example, employing the eddy-dissipation model gives the following reaction-rate relation for in-fluid mass fraction of unreacted fuel as CVA:
Rc_{gen,k} = R_{fu,k} = - Are/K min( C_{fu,k}, C_{ox,k}/s ) m_{k}
where
If the heat losses, say, to the cold walls can not be neglected compared with the heat release through reaction, then the specific gas enthalpy can be treated as CVA, H1_{k}, of the fluid.
The source term of the within-fluid heat losses to wall can be expressed as:
Rc_{gen,k} = R_{H1,k} = S_{H1,k}m_{k}
where net rates of heat transfer to the wall, S_{H1,k}, can be readily accounted for via near-wall heat transfer coefficients, either explicitly or in terms of wall functions employed.
More examples of Rc_{gen,k} formulations will be shown in later sections.
For the systems considered here, the boundaries will be of three general types: walls, fuel/oxidant supply inlets, including distributed ones (e.g., interphase transfer from dispersed phase), and outlets. The following boundary conditions are generally applied:
For the m_{k}'s conservation equations,the unity values of either m_{1} or m_{N} with zero values for the rest of the population are assumed to control the latter at the oxidant and fuel supply regions, correspondingly.
Often the fluid attributes, either PDA or CVA, require further processing in order to derive other interesting quantities, for example:
fª_{mix}= S_{k}m_{k}f_{mix,k}
where m_{k} are used as the weighting coefficients to determine the fluid-averaged mean values.
The population-average value of a CVA, such as C_{k}, is given by a similar relation:
Cª_{k}= S_{k}m_{k}C_{k}
The root-mean-square of the fluctuations for any attribute, f can be written as
RMS_{f} = (S_{k} (fª_{mix}-f _{mix,k})²m_{k})½
An example of auxiliary attribute is a fluid temperature, T_{k}. It can be calculated via two within-fluid CVAs, namely, total enthalpy, H1_{k}, and mass fraction of fuel, C_{fu,k} as follows:
T_{k} = H1_{k} - C_{fu,k}H°
wherein H° represents the specific heat of fuel combustion.
In some cases, e.g. non-adiabatic "mixed-is-burned" combustion, T_{k} manifests itself as a joint function of both PDA, f_{mix,k}, and CVA, H1_{k}, whereas, for the adiabatic infinitely-fast-reaction, combustion system it becomes a function of mixture fraction only.
The Multi-Fluid Model is here applied to steady-state flow, heat and mass transfer in a burner of typical design. Special attention is given to the formation of carbon monoxide and nitric oxides.
The objectives of the study are:
The case considered has the following features:
The task is to calculate the temperatures and composition of combustion gases, along with all related flow properties.
The independent variables of the problem are the three cartesian coordinate X, Y and Z.
A five-fluid population is introduced to represent the actual gas mixture, in which the 5th fluid is associated with inlet fuel-bearing material while the 1st one is considered to be the inlet oxidant substance.
The remaining fluids represent the intermediate states, uniformly distributed over the population space.
The micro-mixing factor, C_{mix}, is supposed to be constant and equal to 10.
The within-fluid properties form the set of CVAs for each fluid, which are solved along with APA, All-Population Attributes, and FMF, Fluid Mass Fractions.
The dependent (solved for) variables are:
The transport equations for APA are solved for the mixture as a whole. The conservation equations for FMF and CVA are solved for all fluids simultaneously.
The nitric oxides are supposed to be presented only in trace quantities, with negligible influence on the state of fluids and within-fluid mass conservation.
For APAs, the transport equations which govern their conservation are of conventional type. They are not described here.
The following formulations are used for:
The source terms featuring in the latter will now be presented in turn.Within-fluid enthalpy-source terms are defined as:
S_{H1,k} = (R^{com}_{CH4,k}H°_{CH4} + R^{com}_{CO,k}H°_{CO})F_{k}
where
Within-fluid methane-source terms are calculated from the rate of combustion as:
S_{CH4,k} = R^{com}_{CH4}F_{k}
Within-fluid oxygen-source terms are defined via net rates as:
S_{O2,k} = R^{net}_{O2,k} F_{k}
Within-fluid carbon-dioxide-source terms are defined via net rates as:
S_{CO2,k} = R^{net}_{CO2,k}F_{k}
Within-fluid carbon-monoxide-source terms are defined via net rates as:
S_{CO,k} = R^{net}_{CO,k}F_{k}
Within-fluid water-vapour-source terms are defined via net rates as follows:
S_{H2O,k} = R^{net}_{H2O,k}F_{k}
Within-fluid nitric-oxide-source terms are defined via formation rates as:
S_{NO,k} = R_{NO,k}F_{k}
Combustion taking place within k-fluid 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 limiting blend of Arrhenius-kinetics and eddy-dissipation rates:
R^{com}_{CH4,k} = - min ( R^{a}_{CH4,k}, R^{e}_{CH4,k} ) and
R^{com}_{CO,k} = - min ( R^{a}_{CO,k}, R^{e}_{CO,k} ),
where R^{a}_{k} and R^{e}_{k}, in kg/m^{3}/s, are the kinetic and eddy-dissipation rates:
R^{e}_{CH4,k} = 4 RHO1 EP/KE
min( FU_{k}, OX_{k}/3)
R^{e}_{CO,k} = 4 RHO1 EP/KE
min( CO_{k}, OX_{k}/0.57)
R^{a}_{CH,k4} =
1.15 10^{9}RHO1^{2}e^{-24444/Tk}
FU_{k}^{-0.3}OX_{k}^{1.3}
R^{a}_{CO,k} =
5.42 10^{9}RHO1^{2}e^{-15152/Tk}
OX_{k}^{0.25}HO_{k}^{0.5}CO_{k}
The remaining rates are defined through associated stoichiometric coefficients:
Step 1:
R^{1}_{O2,k} =
3 R^{com}_{CH4,k}
R^{1}_{CO,k} =
-1.75 R^{com}_{CH4,k}
R^{1}_{H2O,k} =
-2.25 R^{com}_{CH4,k}
Step 2:
R^{2}_{O2,k} =
0.57 R^{com}_{CO}
R^{2}_{CO2,k} =
-1.57 R^{com}_{CO}
The net rates of species generation are the balances of formation and combustion as appropriate:
R^{net}_{CH4,k} =
R^{com}_{CH4,k}
R^{net}_{CO,k} = R^{com}_{CO,k} +
R^{1}_{CO,k}
R^{net}_{O2,k} =
R^{1}_{O2,k} +
R^{2}_{O2,k}
R^{net}_{CO2,k} =
R^{2}_{CO2,k}
R^{net}_{H2O,k} =
R^{1}_{H2O,k}
NOX formation is calculated by incorporating the simple realistic model for the rate of the oxidation of atmospheric nitrogen present in the combustion air. It is known as the steady-state simplification of Zeldovich mechanism with partial-equilibrium assumptions.
The NOX formation rate within k-fluid, in kg/m^{3}/s, is given by:
R_{NO,k} = 2 RHO1 K_{1,k} NN2_{k} [O]_{k} M_{NO}/M_{N2}
where M_{NO} = 30 is NO molecular mass,
M_{N2} =28 is N_{2} molecular mass and
K_{1,k} = 1.8 10^{8}e^{-38370/Tk},
stands for the reaction-rate constant for the forward
reaction:
N_{2} + O ==> NO +H
The equilibrium O-atom concentration, [O]_{k}, can be obtained from the expression:
[O]_{k} = 3.97 10^{5}T_{k}^{-0.5} (RHO1 OX_{k}/M_{O2}) ^{0.5}e^{-31090/Tk}
wherein M_{O2} =32 is the molecular mass of oxygen.
The gas mixture density is computed from the ideal-gas law.
The fluid-specific enthalpies are related to fluid temperatures, T_{k}, and fluid specific heat:
H1_{k} = C_{P,k}T_{k}
The k-fluid specific heat, C_{P,k}, is assumed to be equal for all gas species and is a function of fluid temperature as follows:
C_{P,k} = 1059 + 0.25( T_{k} - 300 )
The population-average values of general fluid attribute, Fª, is computed as:
Fª= SF_{k} F_{k}
The RMS of the attribute fluctuations is calculated as
RMS_{F} = (S (Fª-F _{k})²F_{k})½
The following boundary conditions are applied:
The plots show the distribution of temperatures, velocities and the gas/fluid compositions within the burner.
The results of single-fluid simulation of the same case are used for comparisons. The details of the single-fluid case can be found here.
Pictures are as follows :
To procure steady monotonic convergence, "false-time-step" relaxation was applied to all dependent variables: the false time step was the order of the time of residence of fluid in typical cell. In addition, linear under-relaxation was applied to the pressure corrections and density.
Around 100 iterative sweeps of the calculation domain were required for reasonable convergence, by which is meant that the maximum residual was less than 0.1%.
An extension of MFM for simulating turbulent combustion was presented. The model utilises the solution for CVA to describe the behaviour of realistic chemistry in turbulent reactive flow together with the prediction of presence probability of material from the fuel-bearing stream cased on a one-dimensional discretisation of fluid-population space.
The results appear qualitatively realistic; and mixture temperatures, velocities and gas compositions are within the expected range.
The energy and mass conservations are strictly maintained both within fluid and for the whole-population behaviour.
The highly non-linear nature of reaction rates appears to be more realistically represented by population averaging rather than operating with the single-fluid mean values of conventional treatment.
It contains PLANT statements, and therefore requires the recompilable version of PHOENICS.
In this section, an extension of MFM to the simulation of thermal radiation in combusting systems is presented and exemplified.
The objective of the study is threefold:
The multi-fluid variant is presented here.
The detailed description of its single-fluid counterpart, and subsequent PHOENICS implementation; can be found these two hyperlinked references: gaseous and particulate combustion.
As usual in MFM, the gas mixture is represented by a multi-fluid population.
For the kth-fluid of the population, the conservation equation for the incident radiation, R_{I,k} in Wm^{-2} is written as:
div ( G_{rad,k}gradR_{I,k} ) + a_{k}( 4sT_{k}^{4} - R_{I,k})m_{k} = 0
where s = 5.68 10^{-8}, is the Stefan-Boltzman constant, Wm^{-2}K^{-4}, and m_{k} is the fluid mass fraction, representing its "presence probability".
The within-fluid exchange coefficient, G_{rad,k}, is expressed by:
G_{rad,k} = ( 3(a_{k}+s_{k})- C_{g,k}s_{k} )^{-1}
where
S_{H1,rad,k} = a_{k} (R_{I,k} - 4sT_{k}^{4})
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,k,wall} = 0.5 e_{w} (4sT_{w}^{4} - R_{I,k})(2 - e_{w})^{-1} m_{k}
where e_{w} is the wall emissivity.
The sources for incident radiation at the inlets and outlets are computed in the manner similar to the walls.
The population-average value of incident radiation, Rª_{I}, is computed as:
Rª_{I} = Sm_{k}R_{I,k}
The use of the model will now be exemplified.
Carbon monoxide burns in a stream of oxygen-rich air within a rectangular furnace chamber.
The fuel and oxidant are introduced through separate streams:
The gases are ignited on the entry and steady-state combustion takes place.
The combustion products, with variable absorption coefficient and constant scattering coefficient, flow in a swirling motion towards four exhaust openings located at the corners of the top wall.
The furnace is symmetrical, wherefore only one quarter of it is simulated and as shown here.
The independent variables of the problem are the two components of cartesian coordinate system, namely X and Y.
The five-fluid population is introduced to represent the actual gas mixture, in which the 5th fluid is associated with inlet fuel-bearing material while the 1st one is considered to be the inlet oxidant substance. The remaining fluids represent the intermediate states uniformly distributed over the population space. The micro-mixing factor, C_{mix}, is supposed to be constant and equal to 5.
The within-fluid properties form the set of CVAs (Continuously-Varying Attributes), for each fluid, which are solved along with APAs (All-Population Attributes) and FMFs, (Fluid Mass Fractions).
The dependent (solved-for) variables are:
The transport equations for APA are solved for the mixture as a whole. The conservation equations for FMF and CVA are solved for all fluids simultaneously.
For APA, the transport equations which govern their conservation are of conventional type. They are not desribed here.
The following formulations are used for:
The source terms featuring in the latter will now be presented in turn.
Within-fluid total-enthalpy-source terms are defined as:
S_{H1,k} = a (CR_{k} - 4 face="symbol">sT_{k}^{4})
Within-fluid fuel-source terms are calculated from the rate of combustion as:
S_{FU,k} = R_{FU,k}F_{k}
Within-fluid oxidant-source terms are defined as:
S_{OX,k} = R_{OX,k}F_{k}
Within-fluid product-source terms are defined as:
S_{PR,k} = R_{PR,k}F_{k}
Note that the solution for PR_{k} is not strictly necessary, for it may be calculated from within-fluid composition as:
PR_{k} = 1 - FU_{k} - OX_{k}
Here, however, the solution for PR_{k} is retained as a check for within-fluid mass conservation.
The within-fluid radiation exchange coefficient, G_{rad}, is defined as:
G_{rad} = ( 3(a+s) )^{-1}
Combustion taking place within k-fluid is treated as a one-step irreversible chemical reaction of carbon monoxide oxidation as follows:
CO + 0.5( O_{2} + 0.86N_{2} ) = CO_{2} + 0.43N_{2}
The stoichiometric ratio, Stoic, is 1, and the reaction rates of combustion, in kgm^{-3}s^{-1}, are obtained from the eddy-dissipation rates:
R_{FU,k} = 4 RHO1 min(FU_{k},OX_{k}) EP/KE
The oxidant is consumed at the same rate:
R_{OX,k} = StoicR_{FU,k} = R_{FU,k}
The rate of product generation is twice as much:
R_{PR,k} = (1+Stoic)R_{FU,k} = 2 R_{FU,k}
The within-fluid absorbtion coefficient is expressed as a linear function of fuel and product mass fractions:
a_{k} = 0.2FU_{k} + 0.1PR_{k}
The scattering is supposed to be isotropic with:
s_{k} = 0.1 and
C_{g,k} = 0
The gas mixture density is computed from the fluid-average temperature and ideal-gas law.
The fluid total enthalpies are related to fluid temperatures, T_{k}, fluid specific heat and fuel mass fractions:
H1_{k} = C_{P}T_{k} + H^{o}FU_{k}
where H^{o} is the heat of combustion, J/kg.
The specific heat, C_{P}, is assumed to be equal for all gas species and constant.
The population-average values of general fluid CVA, Fª, is computed as:
Fª= SF_{k} F_{k}
wherin Fs stand for CR_{k}, H1_{k}, FU_{k}, OX_{k} and PR_{k}, as appropriate.
The following boundary conditions are applied:
The plots show the flow distribution, gas temperature and mixture composition within the furnace, as represented by the MFM and single-fluid treatment.
Pictures are as follows :
An extension of MFM for simulating turbulent combustion with significant radiation has been presented.
The CVA concept has been used to describe the behaviour of chemistry radiation.
A one-dimensional five-fluid model has been used.
The results appear qualitatively realistic; and mixture temperatures, velocities and gas compositions are within the expected range.
The energy and mass conservations are strictly maintained both within fluid and for the whole-population behaviour.
The highly non-linear nature of reaction rates and radiation related sources appears to be more realistically represented by population-averaging rather than operating with the single-fluid mean values of conventional treatment.
It contains PLANT statements, and therefore requires the recompilable version of PHOENICS.
In section 4 of this report, MFM is applied to flow, heat and mass transfer in the combustion chamber of a wall-fired furnace. It is relevant to coal-fired boilers and pays special attention to NOX-pollutant formation.
The section presents the implementation in PHOENICS of an IPSA-based Multi-Fluid Model for non-equilibrium two-phase reactive flow of combustible particles finely dispersed in the carrying air stream.
The processes accounted for are:
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 a blended mechanism involving oxygen diffusion to the particle and chemical kinetics.
The relative amount of char combustion products are assumed to be dependent on the particle temperature.
The NOx formation is represented by several simplified sub-models, 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 hydrodynamic turbulence is accounted for by conventional K-e model.
The radiation is modelled via a composite-radiosity model, appropriately modified to account for radiating particles and gases together.
The independent variables of the problem are the three cartesian coordinates X,Y and Z.
The main dependent (solved-for) variables are:
A five-fluid population is introduced to represent the actual gas mixture, in which:
The within-fluid properties form the set of CVAs, for each fluid, which are solved along with APAs (All-Population Attributes), and FMFs (Fluid Mass Fractions).
The dependent (solved for) variables are:
The following formulations are used for:
The source terms featuring in the latter will be presented in the later sections.The transport equations for APA are solved for the each of the thermodynamic phases as a whole.
The conservation equations for FMF and CVA are solved for all fluids in gas phase simultaneously.
The nitric oxides and its intermediates are supposed to be presented only in trace quantities, not contributing to the state of gas mixture fluids and within-fluid mass conservation.
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 is burnt in the char combustion. The ash contents of the phase, ASH2, increase to value 1 in a completely burned particles.
For the first, gaseous phase volatile fuel, CH4, appears as a consequence of devolatilisation, and then disappears as it is burned in gas-phase combustion. The oxygen contents, O2, decrease as it is consumed in both volatile and char combustion.
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 taking place within k-fluid of a gas phase 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 limiting blend of Arrhenius kinetics and eddy-dissipation rates:
R^{com}_{CH4,k} = - min ( R^{a}_{CH4,k}, R^{e}_{CH4,k} ) and
R^{com}_{CO,k} = - min ( R^{a}_{CO,k}, R^{e}_{CO,k} ),
where R^{a}_{k} and R^{e}_{k}, in kg/m^{3}/s, are the kinetic and eddy-dissipation rates:
R^{e}_{CH4,k} = 4 RHO1 EP/KE
min( FU_{k}, OX_{k}/3)
R^{e}_{CO,k} = 4 RHO1 EP/KE
min( CO_{k}, OX_{k}/0.57)
R^{a}_{CH,k4} =
1.15 10^{9}RHO1^{2}e^{-24444/Tk}
FU_{k}^{-0.3}OX_{k}^{1.3}
R^{a}_{CO,k} =
5.42 10^{9}RHO1^{2}e^{-15152/Tk}
OX_{k}^{0.25}HO_{k}^{0.5}CO_{k}
The remaining rates are defined through associated stoichiometric coefficients:
Step 1:
R^{1}_{O2,k} =
3 R^{com}_{CH4,k}
R^{1}_{CO,k} =
-1.75 R^{com}_{CH4,k}
R^{1}_{H2O,k} =
-2.25 R^{com}_{CH4,k}
Step 2:
R^{2}_{O2,k} =
0.57 R^{com}_{CO}
R^{2}_{CO2,k} =
-1.57 R^{com}_{CO}
The net rates of species generation are the balances of formation and combustion as appropriate:
R^{net}_{CH4,k} =
R^{com}_{CH4,k}
R^{net}_{CO,k} = R^{com}_{CO,k} +
R^{1}_{CO,k}
R^{net}_{O2,k} =
R^{1}_{O2,k} +
R^{2}_{O2,k}
R^{net}_{CO2,k} =
R^{2}_{CO2,k}
R^{net}_{H2O,k} =
R^{1}_{H2O,k}
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 a harmonic blend of a kinetic- and diffusion-controlled mass-transfer rates:
R_{C} = A_{p} [ ( K^{k}_{C})^{-1} + (K^{d}_{C} )^{-1} ]^{-1}P_{O2}
where
A_{p} = 6R2/SIZE, is volumetric particle surface area, 1/m, and P_{O2} is partial pressure of O_{2}, N/m^{2}.
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.
Thermal radiation is modelled by a multi-fluid variant of the P-1 radiation model. The detailed description of its single-fluid counterpart and subsequent PHOENICS implementation can be found elsewhere both for gaseous and particulate combustion.
For k-fluid of the population the incident radiation, CR_{k} or R_{I,k} in Wm^{-2}, the conservation equation of the P-1 model is written as:
div ( G_{rad,k}gradR_{I,k} ) + a_{k}( 4sT_{k}^{4} - R_{I,k})m_{k} = 0
where s = 5.68 10^{-8}, is Stefan-Boltzman constant, Wm^{-2}K^{-4}, and m_{k} is the fluid mass fraction, representing its "presence probability".
The within-fluid exchange coefficient, G_{rad,k}, is expressed by:
G_{rad,k} = ( 3(a_{k}+s_{k})- C_{g,k}s_{k} )^{-1}
where
NOx is the term used for anumber of nitrogen oxidation products present in the gases emerging from combustion. They are considered as hazardous pollutants.
NOx emission 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 low concentration trace quantities, the NOx modeling is assumed 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
NOX formation is calculated by incorporating the simple realistic model for the rate of the oxidation of atmospheric nitrogen present in the combustion air. It is known as the steady-state simplification of Zeldovich mechanism with partial-equilibrium assumptions.
The NOX formation rate within k-fluid, in kg/m^{3}/s, is given by:
R_{NO,k} = 2 RHO1 K_{1,k} NN_{k} [O]_{k} M_{NO}/M_{N2}
where M_{NO} = 30 is NO molecular mass,
M_{N2} =28 is N_{2} molecular mass and
K_{1,k} = 1.8 10^{8}e^{-38370/Tk},
stands for the reaction-rate constant for the forward
reaction:
N_{2} + O ==> NO +H
The equilibrium O-atom concentration, [O]_{k}, can be obtained from the expression:
[O]_{k} = 3.97 10^{5}T_{k}^{-0.5} (RHO1 OX_{k}/M_{O2}) ^{0.5}e^{-31090/Tk}
wherein M_{O2} =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:
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_{p}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
Only the model sources which required non-standard settings are described below.
Interphase mass transfer, kg/s
CMDOT = V_{cell}( R_{V} + R_{C} )
Wherein V_{cell} is a cell volume, m^{3}
2nd phase composition, kg m^{-3}s^{-1}
1st phase composition, kg m^{-3}s^{-1}
Composite radiosity, W m^{-3}
The source term for the transport equation of R_{CR} is:
S_{R} = 4( a_{p}E_{p} + a_{g}E_{g}) - 4 ( a_{p} + a_{g} )R_{CR}
At the inlets the following boundary-source term per unit area is also applied:
S_{R,in} = ( a_{p} + a_{g} + s_{p} ) sT_{in}^{4} -R_{CR}
The symmetry plane and perfectly-reflecting walls 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}
Enthalpy-source term due to radiation is:
S_{H1,rad} = 4 a_{g}( R_{CR} - E_{g})
The interphase heat transfer contribution is:
S_{H1,int} = HCOF( T_{2} - T_{1} )
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,rad} = 4 a_{p}( R_{CR} - E_{p}) and
S_{H2,int} = HCOF( T_{1} - T_{2} )
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}
Volumetric heat transfer coefficient, HCOF, in W m^{-3} K^{-1}, is then evaluated as:
HCOF = A_{p}k_{gas}SIZE^{-1}NUSS
where, k_{gas} is gas heat conductivity taken as constant.
The furnace chamber chosen for the present computations, idealised though it is, retains many features of real-life industrial equipment:
Because the furnace is symmetrical, only half of it is simulated.
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.
The plots show the flow distribution, mixture composition as represented by the model, gas temperature and velocity within the combustion chamber.
Pictures are as follows :
An extension of MFM for simulating turbulent combustion of pulverised coal has been presented.
The model utilises the solution for CVAs to describe the behaviour of realistic chemistry in turbulent reactive flow, together with the prediction of presence probability of material from the fuel-bearing stream on the grounds of one-dimensional discretisation of fluid population space.
The results appear qualitatively realistic; mixture temperatures, velocities and gas compositions are within the expected range.
The energy and mass conservations are strictly maintained both within fluid and for the whole-population behaviour.
The highly non-linear nature of reaction rates appears to be more realistically represented by population-averaging rather than operating with the single-fluid mean values of conventional treatment.
It contains PLANT statements, and therefore requires the recompilable version of PHOENICS.
BY : Dr S V Zhubrin, CHAM Ltd
DATE : December 2001
FOR : Demonstration case for V3.4
The case presents the implementation in PHOENICS of P-1 radiation model. It is the simplest case of the more general P-N model based on the first order spherical harmonic expansion of the radiation intensity. The model is capable of taking into account the effect of anisotropic scattering in absorbing, emmiting and scattering media.
The model is easy to implement and solve with little computing efforts. It can also be applied to complex BFC geometries and is easily extendable to handle the particulate effects.
The model offers the advantages of simplicity, high computational efficiency and relatively good accuracy, if the optical thickness is not too small.
The model is applied to the reactive flow in a three-dimensional gas-fired stove.
The independent variables of the problem are the three components of cartesian coordinate system, namely X, Y and Z.
The main dependent (solved for) variables are:
The K-epsilon model, KEMODL, closed by wall functions is used to calculate the distribution of turbulence energy and its dissipation rate from which the turbulence viscosity is deduced.
Combustion is treated as a single-step irreversible diffusion-controlled chemical reaction with a infinitely fast rate between fuel and oxidant. The gas composition and its enthalpy are related to the mixture fraction according to the Simple Chemical Reaction Scheme, SCRS, concept.
The gas density is computed from the local pressures, gas temperatures and local mixture molecular masses.
The specific enthalpies are related to gas temperatures, fuel mass fraction and the heat of combustion.
Thermal radiation is modelled by the expanding the radiation intensity in terms of first order spherical harmonics. Assuming that only four terms representing the moments of the intensity are used, this leads to a incident radiation, CRAD or R_{I}, W m^{-2},conservation equation:
div ( G_{rad}gradR_{I} ) + a ( 4sT^{4} - R_{I}) = 0
where s = 5.68 10^{-8}, is Stefan-Boltzman constant, W m^{-2} K^{-4}.
The exchange coefficient, G_{rad}, is expressed by:
G_{rad} = ( 3(a+s)- C_{g}s )^{-1}
where
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
The volumetric source term for the mixture enthalpy, due to radiation, is then given by:
S_{H1,rad} = a (R_{I} - 4sT^{4})
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 emissivity.
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 emissivity 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.
The furnace chamber comprised a rectangular enclosure. The fuel, CH_{4}, and oxidant, air, are introduced through separate streams. The fuel is introduced through the duct on the middle of the bottom wall, and the heated air is introduced through the inlet slots surrounding the fuel duct.
The gases are ignited on the entry and steady state combustion takes place. The combustion products with constant absorption coefficient of 0.4 and constant scattering coefficient of 0.1 flow, in a swirling motion, towards four exhaust openings located at the corners of the top wall.
The furnace is symmetrical, and only quarter of it is simulated and represented by computational domain as shown here.
The walls, inlets and outlets are considered as perfectly reflecting boundaries.
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.
The plots show the flow distribution, gas temperature and mixture composition within the furnace, as represented by the model in the absence of radiation and with radiation included.
Pictures are as follows :
All model settings have been made by PIL commands and PLANT settings of PHOENICS 3.4
BY : Dr S V Zhubrin, CHAM Ltd
DATE : January 2002
FOR : Demonstration case for V3.4
The implementation in PHOENICS of a radiation model which is based on P-1 approximation and accounts for particulate effects is presented. The P-1 is the simplest case of the more general P-N modelling approach based on the first order spherical harmonic expansion of the radiation intensity.
The current model is an extension of P-1 approximation suitably generalized to handle the effects of absorbing, emmiting, and scattering particles on anisotropic scattering in absorbing, emmiting, and scattering media.
The model is easy to implement and solve with little computing efforts. It is readily applied to complex BFC geometries.
The model offers the advantages of simplicity, high computational efficiency and relatively good accuracy, if the optical thickness is not too small.
The implementation of the radiation model for the pulverized coal-combustion in a wall-fired furnace is demonstrated.
Thermal radiation is modelled by the expanding the radiation intensity in terms of first order spherical harmonics.
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_{g}^{4} - R_{I}) + a_{p} ( 4sT_{p}^{4} - R_{I}) = 0
where s = 5.68 10^{-8}, is Stefan-Boltzman constant, W m^{-2} K^{-4}, T_{p} and T_{g} 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 emissivity 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
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_{g}^{4}) + a_{p} (R_{I} - 4sT_{p}^{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_{p}^{4})
where A_{s} is the volumetric particle surface area, m^{-1} calculated as follows:
A_{s} = 6r_{2}d_{p}^{-1}
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 emissivity.
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 emissivity 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.
The radiation model is implemented for the combustion of pulverized coal in a wall-fired furnace as described in details here. The composite-radiosity model of the original simulation is replaced by P-1 formulations detailed above.
All model settings have been made by PIL commands and PLANT settings of PHOENICS 3.4