Multi-Fluid Modelling of Turbulent Combustion

by

S.V.Zhubrin

Edited by D.B.Spalding, Feb 13-18, 2003

Abstract

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:



Main text

Contents

  1. Definition of the model
  2. MFM for Two-Step Combustion Reactions
  3. MFM for Thermal Radiation
  4. MFM for NOX Formation in a Wall-Fired Boiler

  5. Conclusions
  6. References
  7. Appendix 1 A radiation example 1
  8. Appendix 2 A radiation example 2


1. Definition of the model

1.1 General features

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:

  1. The flowing and combusting mixture is treated as a population of interspersed fluids characterised by:
    1. their mass fractions and
    2. their other relevant attributes

  2. A fluid property used for distinguishing one fluid from another is known as a PDA (i.e. Population Distinguishing Attribute).
    Such properties are 'discretised', in the sense that one fluid is distinguished from another by having a PDA value lying between fluid-specific upper and lower values.

  3. All non-discretised properties are known as CVAs (i.e. continuously-varying attributes).

  4. Both mass fractions and CVAs obey their own conservation principles, expressible as diffferential 'transport equations'.

  5. Also formulated and solved are equations for the corresponding population-average properties.
    These may be the averages of either PDAs or CVAs.

  6. The fluids interact with each other through the exchange of mass and energy.

  7. Solution of the mass-fraction equations leads to 'fluid-population distributions' (FPDs) for each point in space and time.

  8. The FPDs can equally well be called 'probability-density functions' (PDFs), because a mass fraction can be interpreted as a 'presence probability' for the fluid in question.

  9. The transport equations which are solved for CVAs contain terms representing exchanges of mass and other entities resulting from inter-fluid encounters.


1.2 Particular features

Assumptions which characterise the particular MFM used here are as follows:

Summary of the particular assumptions


1.3 Population-Defining and Continuously-Varying Attributes

The PDA is uniformly discretized; therefore each fluid of the population has its own average value of mixture fraction uniquely pre-defined as follows:

fmix,k = ( k-1 )/( Nfluid - 1 )

where k is a current fluid number and Nfluid 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, fmix,1 = 0, while the largest-number one represents the pure fuel, fmix,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.


1.4 Fluid-Mass-Fraction Conservation

The mass fraction of each fluid, or its "presence probability", mk, 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(rmk)/dt+ div(rVmk- Gt grad mk ) = Rm,k

Rm,k, the net rate of k-fluid generation, is the balance of micromixing rate, Rmix,k, and interphase transfer, Sp,k:

Rm,k = Rmix,k + Sp,k

The source term, Sp,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, Rmix,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:

Rmix,k = r Si Sj Fk,i,j mimj Ti,j

wherein:

For all the computations reported below, Ti,j is assumed to be independent of i and j; and calculated as inversely proportional to the eddy-break-up time scale:

Ti,j = Cmixe/K

with K standing for the kinetic energy of turbulence, e for its dissipation rate and Cmix for an empirical constant.

The fractional loss of mass is computed by following rules:

Fk,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:

Rmix,1 = -0.5(m3+ m4+ m5)m1
Rmix,2 = m1m3+ m1m4/2+ m1m5/3- 0.5(m4+ m5)m2
Rmix,3 = m4(m1/2+m2)+ m5(m1/3+m2/2)- 0.5(m1+ m5)m3
Rmix,4 = m3m5+ m2m5/2+ m1m5/3- 0.5(m1+ m2)m4
Rmix,5 = -0.5(m1+ m2+ m3)m5


1.5 Transport Equations for CVAs

Let Ck be the value of a continuously-varying attribute of fluid k. The conservation equation for Ck takes the following general form:

d(rCk)/dt+ div(rVCk- Gt grad Ck ) = Rck + Sck,p

where Rck is within-fluid mass rate of creation and depletion of Ck and Sck,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, Rcmix,k, and the source of CVA due to its in-fluid generation and/or dissipation, Rcgen,k:

Rck = Rcmix,k + Rcgen,k

The contributions resulting from micro-mixing by fluid encounters are written as:

Rcmix,k= SiMk,i(Ci-Ck)

where Mk,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:

SiMk,i= rSi Sj Fck,i,j mimjTi,j

where
Fck,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:

Rcmix,1 = 0 ;
Rcmix,2 =  (m1m3/2+m1m4/4+m1m5/6) (C1-C2)
                                         +m1m3/2  (C3-C2)
                                         +m1m4/4  (C4-C2)
                                         +m1m5/6  (C5-C2) ;
Rcmix,3 =                (m1m4/4+m1m5/6) (C1-C3)
                           +(m2m4/2+m2m5/4) (C2-C3)
                           +(m1m4/4+m2m4/2) (C4-C3)
                           +(m1m5/6+m2m5/4) (C5-C3) ;
Rcmix,4 =                               m1m5/6  (C1-C4)
                                          +m2m5/4  (C2-C4)
                                          +m3m5/2  (C3-C4)
              +(m3m5/2+m2m5/4+m1m5/6) (C5-C4) ;

Rcmix,5 = 0

The generation/dissipation rates, Rcgen,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:

Rcgen,k = Rfu,k = - Are/K min( Cfu,k, Cox,k/s ) mk

where

Here, the rate of reaction has been taken as proportional to the "presence probability", mk, so as to preclude the fuel depletion in non-existent fluid.

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, H1k, of the fluid.

The source term of the within-fluid heat losses to wall can be expressed as:

Rcgen,k = RH1,k = SH1,kmk

where net rates of heat transfer to the wall, SH1,k, can be readily accounted for via near-wall heat transfer coefficients, either explicitly or in terms of wall functions employed.

More examples of Rcgen,k formulations will be shown in later sections.


1.6 Boundary Conditions for mk and CVA

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:

The appropriate values of CVAs are specified as consistent with their nature.


1.7 Derived Attributes

Often the fluid attributes, either PDA or CVA, require further processing in order to derive other interesting quantities, for example:

The population-average value of a discretized variable, such as mixture fraction, is given by:

mix= Skmkfmix,k

where mk are used as the weighting coefficients to determine the fluid-averaged mean values.

The population-average value of a CVA, such as Ck, is given by a similar relation:

k= SkmkCk

The root-mean-square of the fluctuations for any attribute, f can be written as

RMSf = (Sk (fªmix-f mix,k)²mk

An example of auxiliary attribute is a fluid temperature, Tk. It can be calculated via two within-fluid CVAs, namely, total enthalpy, H1k, and mass fraction of fuel, Cfu,k as follows:

Tk = H1k - Cfu,k

wherein H° represents the specific heat of fuel combustion.

In some cases, e.g. non-adiabatic "mixed-is-burned" combustion, Tk manifests itself as a joint function of both PDA, fmix,k, and CVA, H1k, whereas, for the adiabatic infinitely-fast-reaction, combustion system it becomes a function of mixture fraction only.


2. MFM for Two-step Reaction of Combustion

2.1 Introduction

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:

  1. to provide an MFM formulation the reaction phenomena involved
  2. to incorporate this formulation into PHOENICS;
  3. to compare the results with single-fluid predictions; and
  4. to investigate the physical realism of the results obtained.

2.2 Problem description

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.

2.3 Theoretical model; dependent and independent variables

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, Cmix, 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.

2.4 The governing equations

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:

SH1,k = (RcomCH4,kCH4 + RcomCO,kCO)Fk

where

combustion for methane and carbon monoxide.

Within-fluid methane-source terms are calculated from the rate of combustion as:

SCH4,k = RcomCH4Fk

Within-fluid oxygen-source terms are defined via net rates as:

SO2,k = RnetO2,k Fk

Within-fluid carbon-dioxide-source terms are defined via net rates as:

SCO2,k = RnetCO2,kFk

Within-fluid carbon-monoxide-source terms are defined via net rates as:

SCO,k = RnetCO,kFk

Within-fluid water-vapour-source terms are defined via net rates as follows:

SH2O,k = RnetH2O,kFk

Within-fluid nitric-oxide-source terms are defined via formation rates as:

SNO,k = RNO,kFk

2.5 Combustion model

Combustion taking place within k-fluid 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 limiting blend of Arrhenius-kinetics and eddy-dissipation rates:

RcomCH4,k = - min ( RaCH4,k, ReCH4,k ) and

RcomCO,k = - min ( RaCO,k, ReCO,k ),

where Rak and Rek, in kg/m3/s, are the kinetic and eddy-dissipation rates:

ReCH4,k = 4 RHO1 EP/KE min( FUk, OXk/3)
ReCO,k = 4 RHO1 EP/KE min( COk, OXk/0.57)
RaCH,k4 = 1.15 109RHO12e-24444/Tk FUk-0.3OXk1.3
RaCO,k = 5.42 109RHO12e-15152/Tk OXk0.25HOk0.5COk

The remaining rates are defined through associated stoichiometric coefficients:

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

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

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

RnetCH4,k = RcomCH4,k
RnetCO,k = RcomCO,k + R1CO,k
RnetO2,k = R1O2,k + R2O2,k
RnetCO2,k = R2CO2,k
RnetH2O,k = R1H2O,k

2.6 NOx formation chemistry

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/m3/s, is given by:

RNO,k = 2 RHO1 K1,k NN2k [O]k MNO/MN2

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

N2 + O ==> NO +H

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

[O]k = 3.97 105Tk-0.5 (RHO1 OXk/MO2) 0.5e-31090/Tk

wherein MO2 =32 is the molecular mass of oxygen.

2.7 Thermodynamic properties and auxiliary attributes

The gas mixture density is computed from the ideal-gas law.

The fluid-specific enthalpies are related to fluid temperatures, Tk, and fluid specific heat:

H1k = CP,kTk

The k-fluid specific heat, CP,k, is assumed to be equal for all gas species and is a function of fluid temperature as follows:

CP,k = 1059 + 0.25( Tk - 300 )

2.8 Population-average attributes

The population-average values of general fluid attribute, Fª, is computed as:

Fª= SFk Fk

The RMS of the attribute fluctuations is calculated as

RMSF = (S (Fª-F k)²Fk

2.9 Boundary conditions

The following boundary conditions are applied:

2.10 Results

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 :

2.11 Convergence and computer requirements

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%.

2.12 Conclusions

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.

2.13 The Q1 file

The q1 file from which the results were generated can be seen by clicking here.

It contains PLANT statements, and therefore requires the recompilable version of PHOENICS.

3. Multi-Fluid Model of Thermal Radiation

3.1 Introduction

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:

  1. to formulate a multi-fluid approach to the widely-used method P-1 radiation model of radiative transfer;
  2. to incorporate it into PHOENICS; and
  3. to assess the physical realism of the simulation results by comparison with conventional, single-fluid, predictions.

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.

3.2 P-1-MFM: model equations

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, RI,k in Wm-2 is written as:

div ( Grad,kgradRI,k ) + ak( 4sTk4 - RI,k)mk = 0

where s = 5.68 10-8, is the Stefan-Boltzman constant, Wm-2K-4, and mk is the fluid mass fraction, representing its "presence probability".

The within-fluid exchange coefficient, Grad,k, is expressed by:

Grad,k = ( 3(ak+sk)- Cg,ksk )-1

where

The volumetric source term for the mixture enthalpy of k-fluid, H1k, due to radiation, is then given by:

SH1,rad,k = ak (RI,k - 4sTk4)

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,k,wall = 0.5 ew (4sTw4 - RI,k)(2 - ew)-1 mk

where ew 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:

I = SmkRI,k

The use of the model will now be exemplified.

3.3 Case study: Problem description

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.

3.4 The dependent and independent variables

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, Cmix, 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.

3.5 The governing equations

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:

SH1,k = a (CRk - 4 face="symbol">sTk4)

Within-fluid fuel-source terms are calculated from the rate of combustion as:

SFU,k = RFU,kFk

Within-fluid oxidant-source terms are defined as:

SOX,k = ROX,kFk

Within-fluid product-source terms are defined as:

SPR,k = RPR,kFk

Note that the solution for PRk is not strictly necessary, for it may be calculated from within-fluid composition as:

PRk = 1 - FUk - OXk

Here, however, the solution for PRk is retained as a check for within-fluid mass conservation.

The within-fluid radiation exchange coefficient, Grad, is defined as:

Grad = ( 3(a+s) )-1

3.6 Combustion model

Combustion taking place within k-fluid is treated as a one-step irreversible chemical reaction of carbon monoxide oxidation as follows:

CO + 0.5( O2 + 0.86N2 ) = CO2 + 0.43N2

The stoichiometric ratio, Stoic, is 1, and the reaction rates of combustion, in kgm-3s-1, are obtained from the eddy-dissipation rates:

RFU,k = 4 RHO1 min(FUk,OXk) EP/KE

The oxidant is consumed at the same rate:

ROX,k = StoicRFU,k = RFU,k

The rate of product generation is twice as much:

RPR,k = (1+Stoic)RFU,k = 2 RFU,k

3.7 Properties of the medium and auxiliary attributes

The within-fluid absorbtion coefficient is expressed as a linear function of fuel and product mass fractions:

ak = 0.2FUk + 0.1PRk

The scattering is supposed to be isotropic with:

sk = 0.1 and
Cg,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, Tk, fluid specific heat and fuel mass fractions:

H1k = CPTk + HoFUk

where Ho is the heat of combustion, J/kg.

The specific heat, CP, is assumed to be equal for all gas species and constant.

3.8 Population-average attributes

The population-average values of general fluid CVA, Fª, is computed as:

Fª= SFk Fk

wherin Fs stand for CRk, H1k, FUk, OXk and PRk, as appropriate.

3.9 Boundary conditions

The following boundary conditions are applied:

3.10 the results

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 :

3.11 Concluding remarks

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.

3.12 The Q1 file

The q1 file from which the results were generated can be seen by clicking here.

It contains PLANT statements, and therefore requires the recompilable version of PHOENICS.


4. MFM for NOX formation in a Wall-Fired Boiler

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.

4.1 Introduction

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.


4.2 Modelling Considerations

4.2.1 Conservation equations

The independent variables of the problem are the three cartesian coordinates X,Y and Z.

The main dependent (solved-for) variables are:

The main auxiliary 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 main auxiliary (stored) attributes are: 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 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.


4.2.2 Main features of the sub-models

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.


4.2.3 Devolatilisation 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


4.2.4 Devolatilisation kinetics

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

RV = 2.103e-2829/T2mV2MVPV

Here

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:

MVPV = RHO2 YVR2


4.2.5 Volatiles combustion

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: 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 limiting blend of Arrhenius kinetics and eddy-dissipation rates:

RcomCH4,k = - min ( RaCH4,k, ReCH4,k ) and

RcomCO,k = - min ( RaCO,k, ReCO,k ),

where Rak and Rek, in kg/m3/s, are the kinetic and eddy-dissipation rates:

ReCH4,k = 4 RHO1 EP/KE min( FUk, OXk/3)
ReCO,k = 4 RHO1 EP/KE min( COk, OXk/0.57)
RaCH,k4 = 1.15 109RHO12e-24444/Tk FUk-0.3OXk1.3
RaCO,k = 5.42 109RHO12e-15152/Tk OXk0.25HOk0.5COk

The remaining rates are defined through associated stoichiometric coefficients:

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

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

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

RnetCH4,k = RcomCH4,k
RnetCO,k = RcomCO,k + R1CO,k
RnetO2,k = R1O2,k + R2O2,k
RnetCO2,k = R2CO2,k
RnetH2O,k = R1H2O,k


4.2.6 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:

w = MCOYCO2/( MCOYCO2 + MCO2YCO )

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 a harmonic blend of a kinetic- and diffusion-controlled mass-transfer rates:

RC = Ap [ ( KkC)-1 + (KdC )-1 ]-1PO2

where

Ap = 6R2/SIZE, is volumetric particle surface area, 1/m, and PO2 is partial pressure of O2, N/m2.

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 )

where

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:

RC,CO2 = wMCO2/MC

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

4.2.7 Thermal radiation

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, CRk or RI,k in Wm-2, the conservation equation of the P-1 model is written as:

div ( Grad,kgradRI,k ) + ak( 4sTk4 - RI,k)mk = 0

where s = 5.68 10-8, is Stefan-Boltzman constant, Wm-2K-4, and mk is the fluid mass fraction, representing its "presence probability".

The within-fluid exchange coefficient, Grad,k, is expressed by:

Grad,k = ( 3(ak+sk)- Cg,ksk )-1

where


4.2.8 The formation of NOx

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, NO2, and nitrous oxide, N2O.

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/m3/s, is given by:

RNO,k = 2 RHO1 K1,k NNk [O]k MNO/MN2

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

N2 + O ==> NO +H

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

[O]k = 3.97 105Tk-0.5 (RHO1 OXk/MO2) 0.5e-31090/Tk

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,HCN = RCmNCMHCN/MN

    where

    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,HCN = RVmNVMHCN/MN

    where

    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 ApPNO

    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


4.3 Computational Details

4.3.1 Model sources

Only the model sources which required non-standard settings are described below.

Interphase mass transfer, kg/s

CMDOT = Vcell( RV + RC )

Wherein Vcell is a cell volume, m3

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

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

Composite radiosity, W m-3

The source term for the transport equation of RCR is:

SR = 4( apEp + agEg) - 4 ( ap + ag )RCR

At the inlets the following boundary-source term per unit area is also applied:

SR,in = ( ap + ag + sp ) sTin4 -RCR

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:

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

The contribution from the heats of reactions is:

SH1,com = RcomCH4CH4 + RcomCOCO + RCC

Enthalpy-source term due to radiation is:

SH1,rad = 4 ag( RCR - Eg)

The interphase heat transfer contribution is:

SH1,int = HCOF( T2 - T1 )

2nd phase energy, W m-3

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

SH2 = SH2,rad + SH2,int

where

SH2,rad = 4 ap( RCR - Ep) and

SH2,int = HCOF( T1 - T2 )

4.3.2 Temperature calculations

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

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

4.3.3 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.


4.3.4 Interphase heat transfer

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

NUSS = 2. +0.65 REYN0.5

Volumetric heat transfer coefficient, HCOF, in W m-3 K-1, is then evaluated as:

HCOF = ApkgasSIZE-1NUSS

where, kgas is gas heat conductivity taken as constant.



4.4 Case Study

4.4.1 Geometry and situation

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.


4.4.2 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.

4.4.3 Simulation results

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 :


4.5 Concluding Remarks

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.

4.6 The Q1 file

The q1 file from which the results were generated can be seen by clicking here.

It contains PLANT statements, and therefore requires the recompilable version of PHOENICS.



5. Conclusions;

The results presented so far can be summarised as follows:


6. References;

  1. D.B. Spalding, "Multi-Fluid Models for Simulating Turbulent Combustion", - Presentation at CODE Annual SEMINAR in Teraelahti, Finland, 3-4 October, 2001, www.cham.co.uk/phoenics/d_polis/d_lecs/turb2001/mfm_comb.htm

  2. D.B. Spalding and S.V.Zhubrin, "Development of Multi-Fluid Turbulence Model and its Engineering Applications", - Presentation at VIII IPUC, Luxembourg, 17-20 May, 2000,

  3. D.B. Spalding, "Multi-Fluid Model of Turbulent Flow, Mixing and Combustion", - 1996, www.cham.co.uk/phoenics/d_polis/d_lecs/mfm/mfmbas1.htm"


Appendix 1

P-1 radiation model: A stove with significant radiation

BY : Dr S V Zhubrin, CHAM Ltd

DATE : December 2001

FOR : Demonstration case for V3.4

INTRODUCTION

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.

MODELLING CONSIDERATIONS

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:

Turbulence and combustion models

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.

Properties and auxiliary relations

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.

P-1 model of thermal radiation

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 RI, W m-2,conservation equation:

div ( GradgradRI ) + a ( 4sT4 - RI) = 0

where s = 5.68 10-8, is Stefan-Boltzman constant, W m-2 K-4.

The exchange coefficient, Grad, is expressed by:

Grad = ( 3(a+s)- Cgs )-1

where

The symmetry factor 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

The volumetric source term for the mixture enthalpy, due to radiation, is then given by:

SH1,rad = a (RI - 4sT4)

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 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.

CASE STUDY

Geometry

The furnace chamber comprised a rectangular enclosure. The fuel, CH4, 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.

Boundary conditions

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 RESULTS

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 :

THE IMPLEMENTATION

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

Appendix 2

Radiating particles in a flame systems

BY : Dr S V Zhubrin, CHAM Ltd

DATE : January 2002

FOR : Demonstration case for V3.4

INTRODUCTION

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.

MODELLING CONSIDERATIONS

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 ( 4sTg4 - RI) + ap ( 4sTp4 - RI) = 0

where s = 5.68 10-8, is Stefan-Boltzman constant, W m-2 K-4, Tp and Tg are the gas and particle temperatures, K.

The exchange coefficient, Grad, is expressed by:

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

where

The equivalent particle absorption coefficient is defined as:

ap = epAp

where ep is the emissivity 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 - 4sTg4) + ap (RI - 4sTp4)

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

SH2,rad = epAs (0.25RI - sTp4)

where As is the volumetric particle surface area, m-1 calculated as follows:

As = 6r2dp-1

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 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 IMPLEMENTATION

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