A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling

ARTICLE IN PRESS

JID: CNF

[m5G;October 14, 2015;22:25]

Combustion and Flame 000 (2015) 1–14

Contents lists available at ScienceDirect

Combustion and Flame journal homepage: www.elsevier.com/locate/combustflame

A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling ˛ a, Roman Weber b,∗ Jakub Bibrzycki a,b, Marco Mancini b, Andrzej Szlek a b

Institute of Thermal Technology, Silesian University of Technology, Konarskiego 22, Gliwice 44-100, Poland Institute for Energy Process Engineering and Fuel Technology, Clausthal University of Technology, Agricolastrasse 4, Clausthal-Zellerfeld 38678, Germany

a r t i c l e

i n f o

Article history: Received 30 November 2014 Revised 20 September 2015 Accepted 21 September 2015 Available online xxx Keywords: Char combustion Fluidized bed CFD Grain-Thermo-Balance

a b s t r a c t The Grain-Thermo-Balance (GTB) rig was built for generation of the char, at high heating rates, and for experimental determination of its oxidation rates. A zero-dimensional mathematical model for calculating oxidation rates of millimeter-size char particles was developed to serve as a char particle sub-model in a CFD-based software for predicting performance of fluidized bed boilers. When morphology of the char was determined using mercury porosimetry and the char oxidation rates in the kinetic regime were measured using TGA, the zero-dimensional Shrinking Particle Model was able to reproduce the GTB measurements well - with exception of the last 2% burnout where due to a model singularity the calculated temperatures exceeded the measured values. When the char was generated in TGA, at low heating rates, its intrinsic reactivity was four times lower and the reactivity decrease was attributed to alterations to the char morphology (40% slow-down) and annealing (60%). © 2015 Published by Elsevier Inc. on behalf of The Combustion Institute.

1. Introduction and objectives In industrial coal-fired fluidized beds, the dense bed is formed using inert (sand) particles typically not larger than 1 mm whereas fuel particles can be as large as 3–8 mm. Bed temperatures are typically in the 1150–1230 K range. Current research efforts, aiming at development of CFD-based mathematical models for performance predictions of fluidized bed combustors [1], concentrate on handling the complex fluid dynamics of two-phase flows. As reviewed by Singh et al. [2], there are three numerical methods available and these are (a) Eulerian–Lagrangian methods which use single particle, particle parcel or a group of particles description of the solid-phase, (b) Eulerian–Eulerian Two Fluid Model methods and (c) DiscreteElement Method which also uses the Eulerian–Lagrangian framework. The three zones of fluidized bed combustion namely, dense bed, splash zone and free board, require different CFD-methods and fluid dynamics of the first two zones is the most difficult to simulate. CFD-simulations of such a complex fluid dynamic structure require a

List of acronyms: BCURA, British Coal Utilisation Research Association; FB, Fluidized Bed; GTB, Grain Thermo-Balance; IR, InfraRed; RPM, Random Pore Model; SCM, Shrinking Core Model; SPM, Shrinking Particle Model; TGA, ThermoGravimetric Analysis/Analyzer. ∗ Corresponding author. E-mail addresses: [email protected] (J. Bibrzycki), [email protected]. ˛ de (M. Mancini), [email protected] (A. Szlek), [email protected], [email protected] (R. Weber).

constant intervention of the modeler to force and secure convergence and computational times are very long. Adamczyk et al. [3] have recently modeled an industrial circulating fluidized bed boiler using a hybrid Euler–Lagrange approach and, since it was not possible to achieve a steady-state convergence, the model was run as timedependent until oscillations around a mean value were observed. One may question such an approach but this is the state of CFD-modeling in fluidized bed technology today. It is then perhaps not too surprising that at this stage of the CFD-modeling, BCURA-type models [4] are used [2,3,5] to compute oxidation rates of char particles. In such models, the oxidation rates are computed by taking into account both oxygen diffusion through the particle boundary layer and an Arrhenius expression accounting for kinetics. The objective of our work is to develop a zero-dimensional mathematical model for calculating oxidation rates of millimeter-size char particles which is to serve as a sub-model in CFD-based software for predicting performance of fluidized beds. The goal is to develop and validate a char oxidation model which goes beyond BCURAformulation [4] by taking into account both reactants transport and reactions inside the porous structure. There exist many models that just do that, as exemplified by works of Simons [6] and Haugen et al. [7]; however, they are too complex to be incorporated into a CFD software handling fluidized bed combustion. In other words, zerodimensionality is the constraint we must adhere to. In reviews by Laurendeau [8] and Smith [9] (see also Weber and Mancini [10]) the framework of the mathematical modeling has been summarized and since then a proliferation of models have been developed which can

http://dx.doi.org/10.1016/j.combustflame.2015.09.024 0010-2180/© 2015 Published by Elsevier Inc. on behalf of The Combustion Institute.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

ARTICLE IN PRESS

JID: CNF 2

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

List of symbols A Aint cash cC cchar CCO2 CO2 cp air (Tp ) cp cr av Deff Dg Dgb Dk dpore av dpore av0 dp dp0 dp ∗ Ea int gash gC Tp

Ash content in coal, % Pre-exponential factor for intrinsic kinetics, m/s Specific heat capacity of ash, kJ/(kg · K) Specific heat capacity of carbon, kJ/(kg · K) Specific heat capacity of solid material in char, kJ/ (kg · K) Molar concentration of carbon dioxide, kmol CO2 / m3g Molar concentration of oxygen, kmol O2 /m3g Specific heat capacity of air at constant pressure, at Tp , kJ/(kg · K) Specific heat capacity of char particle, kJ/(kg · K) Average char conversion rate, kg/s Effective diffusivity, m2 /s Gas diffusivity, m2 /s Gas diffusivity at bulk temperature, m2 /s Knudsen diffusivity, m2 /s Average pore diameter at time, τ m Initial average pore diameter, m Particle diameter, m Initial particle diameter, m Equivalent particle diameter, m Activation energy, kJ/kmol Ash mass fraction in char, kg ash/kg char Carbon mass fraction in char, kg C/kg char

hCO

Specific enthalpy of carbon dioxide at Tp , kJ/kg

h298.15 CO2

Specific enthalpy of carbon dioxide at 298.15, K kJ/kg

Tp hCO h298.15 CO

Specific enthalpy of carbon monoxide at Tp kJ/kg

hTO2b

Specific enthalpy of oxygen at Tb , kJ/kg

h298.15 O

Specific enthalpy of oxygen at 298.15, K kJ/kg

2

2

Specific enthalpy of carbon monoxide at 298.15, K kJ/ kg

hr Reaction enthalpy of carbon oxidation, kJ/kg C hr (CO) Reaction enthalpy of carbon oxidation to CO, kJ/kg C hr (CO2 ) Reaction enthalpy of carbon oxidation to CO2 , kJ/kg

kint kpseudo LV m M m0 mash MC MCO MCO2 MO2 p Sext Sm0 SV SV0 TA Tb Tp Twall t V

C Intrinsic kinetics constant, m/s Pseudo-kinetics constant, m/s Volumetric pore length, m/m3 Mass of char sample at time t, kg Moisture content in coal, % Initial mass of char sample, kg Mass of ash, kg Molar mass of carbon, kg C/kmol C Molar mass of carbon monoxide, kg CO/ kmol CO Molar mass of carbon dioxide, kg CO2 /kmol CO2 Molar mass of oxygen, kg O2 /kmol O2 Absolute pressure, Pa External particle surface area, m2 Initial char specific internal surface area, m2 /g Char volumetric internal surface area at time t, m2 / m3 Initial char volumetric internal surface area, m2 /m3 Activation temperature, K Gas temperature at bulk conditions, K Particle temperature, K Pipe internal surface temperature, K Time, s Volatiles content in coal, %

Particle volume, m3 Total intrusion volume of mercury per gram of char, m3 /g Volume of char sample, m3 Vsample w Average velocity from outflow of the ceramic pipe, m/s X Char burnout Char burnout in 20–95% range X∗ Molar content of N2 in the gas mixture zN2 zO2 Molar content of O2 in the gas mixture DTG Normalized sample mass loss rate, 1/s HHV Higher heating value of coal, kJ/kg LHV Lower heating value of coal, kJ/kg c carbon mass fraction in coal, % h hydrogen mass fraction in coal, % s sulfur mass fraction in coal, % n nitrogen mass fraction in coal, % o oxygen mass fraction in coal, % Nu Nusselt number Re Reynolds number Pr Prandtl number Sc Schmidt number Sh Sherwood number α Convective heat transfer coefficient, kW/(m2 · K) β eff Effective mass transfer coefficient of O2 , m/s δm−mod Relative difference of average conversion rates, %  Char porosity at time t, m3g /m3 0 Initial char porosity, m3g /m3 p Emissivity of the char particle surface  p−pw Mutual emissivity of the char particle-wall arrangement η Effectiveness factor μb Gas dynamic viscosity at bulk temperature, Pa · s λb Gas thermal conductivity at bulk temperature, kW/ (m · K) γ Stoichiometric coefficient in reaction (3), kmol O2 / kmol C γCO2 = 1 Stoichiometric coefficient of carbon/CO2 reaction kmol, CO2 / kmol C ρ app Apparent char density at t, kg/m3 ρ app0 Initial apparent char density, kg/ m3 φ Thiele modulus ψ Pore structural parameter ρ air (Tp , p) Density of air calculated at p = 1bar and Tp , kg/m3 ρ true True char density, kg/m3solid σ Stefan-Boltzmann constant, kW/ m2 · K4 τ pore Pore tortuosity Vp Vpore

Subscripts 0 Initial app Apparent b Bulk g Gas p Particle Superscripts ad Air dried conditions daf Dry ash-free basis exp Experimental mod Modeling be classified into shrinking particle models (SPM) and shrinking core models (SCM). In the latter models, the particle diameter remains constant during oxidation, and the reactions take place at the sharp interface between the unreacted core and the ash layer. The ash layer

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

ARTICLE IN PRESS

JID: CNF

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14 Table 1 Proximate and ultimate analysis of coal “Janina.” Parameter Moisture content Ash content Volatiles content Higher heating value Lower heating value Carbon Hydrogen Nitrogen Sulfur Oxygen

Symbol ad

M Aad Vad (Vdaf ) HHVad LHVad cad had nad sad oad

Unit

Value

% % % kJ/kg kJ/kg % % % % %

11.0 12.3 30.4 (39.7) 22 910 21 831 58.2 4.0 0.9 1.5 12.1

ad – air dried; daf – dry and ash free.

remains on the particle outskirts causing an extra diffusion resistance for the reactants to reach the reaction front. Such a SCM has recently been used by Liu et al. [5] in 15 s simulations of a fluidized bed combustion of 1.5 mm char particles using DEM-CFD approach. In fluidized bed boilers the ash layer formed on the particle outskirts is continuously removed by attrition and in our opinion the shrinking particle models are more suitable than shrinking core models. Thus, in this work we develop a Shrinking Particle Model with combustion reactions occurring both on the external surface and inside the porous structure; the ash formed on the particle outskirts is being permanently removed causing the particle diameter to decrease. In order to validate our model we have developed an experimental rig called Grain-Thermo-Balance (GTB). The rig is used both for generation of the char at 650–800 K/min heating rates and for experimental determination of oxidation rates of 5.2–6.3 mm char particles. We also generate char samples at 10 K/min heating rates using Thermogravimetry (TGA) in order to compare both char morphology and reactivity of chars generated at different heating (devolatilization) rates. 2. Experimental

coal mine, which is located in Libia˛ z˙ in Lesser Poland Voivodeship and lies on the Upper Silesian Coal Basin. Industrial reserves of this coal field are the biggest among all the running coal mines in Poland and were estimated in 2012 to amount to 357 mln Mg [11]. Table 1 shows proximate and ultimate analysis of coal “Janina” as given by the mine. In reality the ash content of the Janina coal varies in the 6–15% range and therefore the ash content of the char particles was determined to be 14% as an arithmetic average of ten char samples, see below. 2.2. Grain Thermo-Balance experiments Experimental investigations on combustion of millimeter-size particles are numerous [12–21] and usually the center particle temperature [12,16–19,21] has been measured using thermocouples. The particle mass loss has been directly measured in a few experiments only [13–16,20]. Instead, the mass loss has been calculated using the measured CO and CO2 concentrations [17–19]. In our work, the particle mass is determined directly using an analytical balance whilst particle surface temperature is measured using an infrared camera. The test rig (see Figs. 1 and 2), named as Grain Thermo-Balance (GTB), is composed of the balance, electric tube furnace, infrared camera (FLIR ThermaCAM SC2000), rotameter, computer and a special stand, which holds the fuel particle. The electrical heater is used to heat up the reacting gas to a predefined temperature up to circa 795 °C. The GTB experiments have been carried out in two steps: as char generation and char combustion studies. For char generation, coal particles in the 0.2–0.3 g mass range (circa 4–6 mm) have been placed on the particle holder and devolatilized in pure nitrogen at 790 °C. The average devolatilization rate, calculated for the 10 –85% volatiles release range, equals to circa 1%/s. Some of such generated char particles have been manually shaped to increase their sphericity (see Table 2) and their diameter has been measured using a caliper. Table 2 lists the particles considered; the diameter of an equivalent sphere (the fifth column of Table 2) is calculated as:

d∗p = (6 · m p /(π · ρapp ))1/3

2.1. Fuel Experimental investigations described in this work have been performed using Polish bituminous hard coal extracted from the “Janina”

3

(1)

with the char apparent density of 1163 kg/m3 (see Table 3). Another ten of such generated char particles have been analyzed and their averaged composition is: ash content 0.128, moisture content 0.015,

Fig. 1. The test rig; 1 - Electric tube furnace, 2 - Analytical balance, 3 - Infrared camera, 4 - Particle holder, 5 - Rotameter, 6 - Computer, 7 - Outlet pipe.

Fig. 2. The scheme of the outlet pipe and the particle holder.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

ARTICLE IN PRESS

JID: CNF 4

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14 Table 2 Mass and diameter of char particles (named as GTB-char). Particle number3

Mass g

Shape

Diameter by caliper, mm

Equivalent diameter, mm

Time1 s

Remains2 %

1 2 3 4 5 6 7 8 9 10

0.085 0.068 0.151 0.126 0.086 0.150 0.112 0.129 0.119 0.121

Spherical Spherical Spherical Spherical Spherical Spherical Original Original Original Original

5.2 4.8 6.7 6.2 5.3 – – – – –

5.2 4.8 6.2 5.9 5.15 6.3 5.7 6.0 5.8 5.8

169 96 171 134 121 226 126 174 132 132

100 15 17 36 5 78 39 95 10 47

1

Time at 85% conversion. Indicative estimate of the % of the the initial char ash which remains on the GTB stand at the end of the experiment – the initial ash content for all particles was assumed to be 14%. 3 Also Experimental Run Number. 2

Table 3 Values obtained directly from porosimeter software. Parameter

GTB-char

TGA-char

Vpore 0 , cm3 /g Sm0 , m2 / g dpore av 0 , nm ρ app 0 at ≈ 1 bar, kg/m3  0 , m3g /m3 msample , g

0.270 26.8 40.2 1163 0.286 0.459

0.166 15.8 41.9 1257 0.191 0.762

volatiles 0.077 so that the average ash content of the generated char particles is 0.14 on the moisture and volatiles free basis. Thus, in the forthcoming discussion the figure of 14% is used for ash content of the GTB char. In the char combustion experiments, air has been used as an oxidizer and its volumetric flow rate and temperature have been maintained constant at 500 dm3 /h (p = 101325 Pa, t = 20 ◦ C) and 795 ◦ C, respectively, which gives the average air velocity of 3.3 m/s. The particle mass is recorded every second while particle surface temperature every three seconds. Measuring surface temperature using infrared camera requires surface emissivity which in this work is assumed to be 0.8 [22–26]. One should realize that the measurements of particle surface temperature are imprecise because the particle emissivity changes with both temperature [27] and char structure alterations during combustion. The 1400 K surface temperature obtained for 0.8 emissivity would change to 1521 K and 1304 K if the particle surface emissivity were to be 0.7 and 0.9, respectively. The precise positioning of the particle inside the outlet pipe is shown in Fig. 2. Figure 3 shows the measured mass loss for the nearly-spherical particles (particles 1–5) listed in Table 2. The smaller the particle, the shorter the combustion time, as seen in many other works [13,17,19,21,28]. Figure 3 also shows the maximum particle surface temperature as a function of time. The maximum particle surface temperature is the highest (pixel) temperature registered by the IR camera focused on the burning particle assuming 0.8 particle surface emissivity. As shown in Fig. 3 the particles are heated from room temperature to the gas temperature. During this time moisture and residual volatiles are released and then char combustion is initialized which is indicated by a rapid temperature increase until a temperature of 1400–1500 K is reached. Then char combustion proceeds during which a slight increase of the maximum recorded surface temperature is observed. Table 2 also lists particles No. 6–10 of original shape and they are non-spherical. Comparison of combustion times, for nearly-spherical and non-spherical char particles, shows that initial shape of the char particle has negligible effect on the combustion time of particles with

Fig. 3. Char mass (dashed lines) and maximum temperature of particle surface (dotted lines) plotted against time of combustion; Numbers 1–5 indicate particles listed in Table 2.

0.07–0.15 g initial mass. For instance, combustion times of particles No. 4, 9 and 10, having similar initial mass (Table 2) despite different shapes, vary only by circa 2 s when taking 85% char burnout as a reference. Table 2 lists also percentage of the particle initial ash which remains on the particle holder at the end of each experimental run. For all but particles No. 1, 6 and 8, the mass of the remains constitutes not more than 47% of the ash content of the char. This is so since the ash formed on the particle surface has been continuously removed by the air flow throughout the GTB experiments, as it has been confirmed by video sequences taken [29]. In this paper, we focus on char combustion conditions when the outer ash layer is being removed so that Shrinking Particle Model is appropriate. Table 2 includes also particles No. 1, 6 and 8 which in terms of their initial mass and diameter are similar to particles No. 5, 3 and 4, respectively. During the GTB experiments of particles No. 1, 6 and 8, it has been observed that the ash layer remained on the surface and at the end of each experimental run, around 100%, 78% and 95% of the initial char ash formed mechanically stable ash-remains of a size similar to the initial char particle diameter. Indeed the ash layer, which remained on the surface throughout the experiments, increased the heat and mass transfer resistance so that the burnout time increased to around 169 s, 226 s and 174 s, if compared to 121 s, 171 s and 134 s for particles No. 5, 3 and 4, respectively (see also mass loss of particles 1 and 5 in Fig. 3). In this paper we do not consider particles No. 1, 6 and 8 any further, although these are interesting cases for examining the extra resistances occurring due to the ash layer. Among the ten particles examined (see

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF

ARTICLE IN PRESS

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

5

Table 2) seven have burned in the Shrinking Particle Mode and it has been impossible to know a priori whether the ash would remain on the particle outskirts or it would be blown away by the gas flow. 2.3. Thermogravimetry The TGA measurements have been carried out in two steps as char generation and char oxidation/gasification studies. 2.3.1. Char generation For char generation in the TGA, coal particles of original shape and mass in the range of 0.15–0.3 g (circa 4–6 mm) have been used. The char generation is performed in pure nitrogen at 79 ml/min volumetric flow rate and with 10 K/min heating rate. The temperature control program starts with sample heating up to 105 °C with 10 K/min, then fifteen minutes isotherm at 105 °C is applied. The next step is sample heating up to 750 °C with 10 K/min heating rate. For the last step a thirty minutes isotherm at 750 °C is used. An average devolatilization rate is 0.05%/s which has been calculated for the 10–85% volatiles release interval. An average composition of the TGA-generated char has been determined to be: 0.017 moisture, volatiles 0.024, ash 0.111 so that 0.12 ash content (moisture and volatiles free basis) is applicable. One of the goals is to examine how devolatilization conditions affect the char combustion rates. To this end, two char types are produced; one in the GTB test rig at high heating rates of 590 K/min (see Section 2.2) and the second one in the TGA, thus at low heating rates of 10 K/min . These two chars are named as GTB-char and TGA-char. 2.3.2. Char combustion in O2 /N2 atmosphere To determine intrinsic kinetics as well as the conditions of transition from intrinsic to apparent kinetics, special attention is needed in carrying out TGA experiments. The TGA has been operated in isothermal mode to eliminate both the complexity and ambiguity in deriving the kinetic parameters at non-isothermal operation [30]. For obtaining reliable data, free of scattering, attention must be paid to the crucible choice [24,31,32], sample size [24,31–35], sample granulation [35], oxidizer partial pressures and flow rate [32], and temperature [33,36–38]. Discussing all these factors is beyond the scope of this paper and can be found in [29]. Both chars (TGA-char and GTB-char) have been milled to the size less than 0.25 mm and 8 mg sample has been placed in the platter crucible. The temperature program starts with a sample heating up with 50 K/min rate to the required temperature (in N2 ) where, for temperature stabilization, a ten minute isotherm is applied. Then the flow of oxidizer, a mixture of 79% of N2 and 21% of O2 , is applied. Each TGA-experimental run has been performed twice and very good repeatability has been obtained [29]. The GTB-char oxidation measurements have been carried out at 350 °C, 400 °C, 450 °C, 600 °C, 750 °C, 850 °C and 950 °C while the TGA-char has been oxidized at 400 °C, 450 °C, 500 °C, 600 °C, 750 °C, 850 °C and 950 °C. The first three temperatures have been chosen in order to perform the oxidation experiments in kinetic-controlled conditions for both chars (the mass loss figures and the char conversion rates can be found in [29]). As shown in Fig. 4 , at 450 °C in the intrinsic kinetic regime, the combustion rates for the TGA-char and the GTB-char are substantially different; the TGA-char burns with a much lower rate. Indeed, Fig. 5 shows that for temperatures lower than around 500 °C the GTB-char burns four times faster than the TGA-char. On the left-vertical axis of Fig. 5 the average char conversion rate (DTGCav ) is plotted which is calculated as an arithmetic average of the char conversion rates in the 10–85% conversion interval. At lower temperatures each increase of the temperature causes the DTGCav to increase significantly until around 600 °C (Fig. 5). From this temperature onwards the conversion rate increases but at a lower rate. At 750–950 °C temperatures the combustion rates of the two chars do not differ significantly (Fig. 5).

Fig. 4. The TGA-char and the GTB-char combustion rates measured in TGA at 450 °C; Black lines indicate normalized char mass and blue lines show the char conversion rates of GTB-char and TGA-char; oxidizer composition: zO2 = 0.21, zN2 = 0.79.

Fig. 5. Average conversion rates of both chars and the conversion rate ratio of GTBchar and TGA-char plotted against temperature; oxidizer composition: zO2 = 0.21, zN2 = 0.79.

The GTB-char to TGA-char conversion rates ratio of 1:1 has been observed which indicates that oxygen diffusion to the sample surface is a rate controlling mechanism. 2.3.3. Char gasification in CO2 /N2 atmosphere The aim of the char gasification study in N2 /CO2 atmosphere is determination of the gasification rates at various temperatures. The transition from kinetically controlled to diffusion/kinetic controlled regime typically occurs in the 900–1100 °C range [32,36,39,40] and it is affected by CO2 partial pressure and particle size. Three temperature points have been chosen for gasification rates determination; those are: 800 °C, 850 °C and 950 °C. The partial pressure of carbon dioxide, similarly to the char combustion case, is set to be 0.21 bar. A 8 mg sample of GTB-char is used and the temperature program is similar to the one used for char oxidation in N2 /O2 atmosphere. The gasification time is roughly equal to 6.5 h for 8 mg char sample at 850 °C temperature and 85 min at 950 °C temperature. 3. Arrhenius constants determination The procedure for determination of intrinsic rates of char combustion with O2 and gasification with CO2 using the isothermal-TGA is presented in this Section. The char oxidation rates have been determined considering the 10–85% char conversion interval. In some

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF 6

ARTICLE IN PRESS

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

intrinsic reactivity, internal surface area development and oxygen concentration.

dm/dt = −Vsample ·  · 1/γ · MC · kint · SV · CO2

(9)

The sample volume (Vsample ) can be calculated knowing the char apparent density according to

Vsample = m0 /ρapp0

(10)

whilst char porosity can be related to the apparent density using

 = 1 − ρapp /ρtrue

(11)

where the true particle density is calculated as

ρtrue = ρapp0 /(1 − 0 )

(12)

The intrinsic reaction rate constant can be expressed using the Arrhenius equation Fig. 6. CO/CO2 ratio and γ according to Arthur [43] and Tognotti et al. [44] (see also Fig. 1 in Shaddix et al. [45]).

works 50% of char burnout was used [31,35,41] while Karlstrom et al. [42] used 75%. The char burnout is calculated as:

X = 1 − m/m0

(2)

3.1. Equations Assuming that char consists of ash and carbon only, char oxidation can be described by the following stoichiometric equation:

C + γ O2 ⇒ (2γ − 1) CO2 + (2–2γ ) CO

(3)

where carbon is oxidized to both carbon dioxide and carbon monoxide, see Arthur [43]. The CO/CO2 ratio depends on the temperature and oxygen concentration; for oxygen mole fractions not larger than around 21%, the CO/CO2 ratio is [43]:

CO /CO2 = 2512 · exp(−6244/Tp )

(4)

The stoichiometric coefficient (γ ) of global reaction (3) is, for a given temperature, calculated as follows:

γ = (1 + 0.5 · (CO/ CO2 ))/(1 + (CO/CO2 ))

(5)

The coefficient γ varies between one, when only carbon dioxide is produced, and 0.5, at sufficiently high temperatures, when CO is produced only. Arthur relationship (4) was determined experimentally using both artificial graphite and coal char made from Warwickshire coal (the Ryder Seam) and 1–2.83 mm particles were examined in the 753– 1173 K temperature range [43]. Experiments of Tognotti et al. [44] were carried out using 180–250 μm Spherocarb char particles in the 670-1670 K temperature range and, for 20% oxygen content in nitrogen, the following relationship was obtained

CO /CO2 = 69.39 · exp(−3070/Tp )

(6)

Shaddix et al. [45] have recently re-examined Tognotti et al. [44] experiments and concluded that the above relationship is accurate if the particle temperature does not exceed 1250 K. We use both Arthur [43] and Tognotti et al. [44] expressions, see Fig. 6. Thus, the Arrhenius constants have been determined for reaction:

C + O2 ⇒ CO2

(7)

and for reaction:

C + 0.5 O2 ⇒ CO

(8)

as well as for reaction (3) with CO/CO2 dependence given either by Eq. (4) or by Eq. (6). If the reaction of char with oxygen is first order [7,38,42,46– 48], Eq. (9) can be then derived, which relates the mass loss to the

kint = Aint · exp(−Ea int /(R · T ))

(13)

The char specific (volumetric) internal surface area changes during oxidation. According to the Random Pore Model [49] the internal surface development is a function of char burnout and is depended on the pore structural parameter (ψ ). The internal surface area is then described using Eq. (14) [50]

SV = SV 0 · (1 − X ) · (1 − ψ · ln(1 − X ))0.5

(14)

while the pore structural parameter can be calculated according to [49] using Eq. (15)

ψ = 4 · π · LV · (1 − 0 )/SV2 0

(15)

The volumetric pore length, assuming the cylindrical pore shape, can be calculated from the total intrusion volume of mercury per gram of char (Vpore ) and from initial average pore diameter (dpore av0 ) according to Eq. (16) [51]. 2 LV = 4 · Vpore · ρapp0 /(dpore av0 · π )

(16)

For char gasification experiments with CO2 the Boudouard reaction is assumed:

C + CO2 ⇒ 2CO

(17)

According to Batchelder et al. [47], this reaction can be assumed to be first order when the temperature is in the 900–1060 ◦ C range. Thus, Eq. (9) has to be only slightly modified to yield:

dm/dt = −Vsample ·  · 1/γCO2 · MC · kint · SV · CCO2

(18)

Other parameters, which are required for intrinsic kinetics of charCO2 reaction, are calculated using Eqs. (10)–(16). 3.2. Mercury porosimetry In order to obtain initial apparent char density (ρ app0 ), initial porosity ( 0 ), initial internal surface area (SV0 ), total mercury intrusion volume of pores (Vpore ) and initial average pore diameter (dpore av0 ), which are necessary in calculations of intrinsic kinetics, the measurements have been carried out using the mercury porosimeter for both chars. The porosimeter allows for pore volume determination up to 3 nm which covers all macropores and most mesopores (2–50 nm). The porosimeter program determines the following parameters: total specific intrusion volume, total specific pore area, median pore diameter, apparent and skeletal densities, average pore diameter calculated as 4 · Vpore0 /Sm0 and total porosity (see Table 3). These parameters (see Eqs. (12), (15), (16), (19) and (20)) allow for calculation of char true density, volumetric pore length, pore structural parameter, initial volumetric internal surface area and pore tortuosity. The initial volumetric internal surface area is calculated as follows:

SV 0 = Sm0 · ρapp0

(19)

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

ARTICLE IN PRESS

JID: CNF

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

7

Table 4 The calculated parameters for GTB-char and TGA-char. Parameter

ρ true , kg/ m LV , m/m3 ψ, − SV0 , m2 / m3 τ pore , −

3

GTB-char

TGA-char

1628 2.47E + 14 2.28 3.11E + 7 1.20

1554 1.51E + 14 3.88 1.99E + 7 1.19

Table 5 Arrhenius constants for GTB-char and TGA-char oxidation using C + O2 ⇒ CO2 reaction. Temperature range and char type

Aint , m/s

Ta int , K

R2 , –

350−450 ◦ C (GTB-char) 400−500 ◦ C (TGA-char) 750−950 ◦ C (GTB-char) 750−950 ◦ C (TGA-char)

9.76E + 03 1.64E + 03 1.84E − 041 1.81E − 04

15970 15320 20001 15701

0.9998 0.9991 0.9984 0.9838

1

Fig. 7. Arrhenius constants determination for C + O2 ⇒ CO2 reaction (see Table 5).

Pseudo-intrinsic constants. Table 6 Arrhenius constants for GTB-char and TGA-char oxidation.

Pore tortuosity, which is defined as the ratio of the real path length of gas diffusion to the straight path from hypothetical points A to B [13], is evaluated as [9]:

τpore = (dpore av0 · SV 0 /(4 · 0 ))2

(20)

Table 3 shows values obtained directly from the porosimeter software, while Table 4 lists the calculated parameters for both chars. 3.3. Calculation results Using both the TGA and the mercury porosimetry data the intrinsic kinetics constants can be calculated. The 10–85% char conversion interval is divided into tens sub-intervals for which kint values are determined using either Eq. (9) (char oxidation with oxygen) or Eq. (18) (char oxidation with CO2 ). Then, an arithmetic average over the subintervals is calculated. Since two measurements have been performed for each temperature, the final kint value is calculated as an average of these two runs. To obtain values of the pre-exponential factor and the activation energy (or activation temperature = Ea int /R) Eq. (13) is linearized:

ln(Aint ) − (Ea int /R) · (1/T ) = ln(kint )

(21)

Plots of ln(kint ) as a function of reciprocal of sample temperature, for two temperature regions (the first one from 350 to 500 ◦ C and the second one from 750 to 950 ◦ C), are shown for reaction (7). Table 5 shows the Arrhenius parameters together with the correlation factor for char-O2 reaction to form CO2 . In the low temperature region the conversion rate increases strongly with the temperature. In the higher temperature region (750–950 ◦ C), the conversion rates are weakly dependent on temperature (Fig. 7) since the process is controlled by diffusion. Similar trends were observed by Jess and Andresen [48] who found five times smaller activation energy in the lower temperature region (700–850 K) in comparison to the high temperature region (950–1100 K) [48]. The ratio of the activation energies calculated for the lower temperature region and the higher temperature region reaches 8.0 and 9.8 for the GTB-char and the TGA-char, respectively (Table 5). The correlation factors are larger than 0.98 which indicates good accuracy of the TGA data. The Arrhenius constants have also been determined for reactions (3) and (8) for the intrinsic temperature region as shown in Table 6. The activation energy for both TGA-char and GTB-char for C + 0.5 O2 ⇒ CO (see Table 6) and C + O2 ⇒ CO2 (see Table 5) global reactions is the same, however the pre-exponential factor is two times smaller

Reaction

Temperature range and char type

Aint , m/s

Ta int , K

R2 , –

C + 0.5 O2 ⇒ CO

350–450 ◦ C (GTB-char) 400–500 ◦ C (TGA-char) 350–450 ◦ C (GTB-char)

4.88E + 03

15, 970

0.9998

8.21E + 02

15, 320

0.9991

4.04E + 03

15, 450

0.9996

400–500 C (TGA-char)

4.79E + 02

14, 560

0.9993

350 − 450 ◦ C (GTB-char)

3.85E + 03

15, 500

0.9997

6.26E + 02

14, 835

0.9991

C + 0.5 O2 ⇒ CO C + γ O2 ⇒ (2 γ −1) CO2 + (2–2 γ ) CO γ based on Arthur, Eq. (4) C + γ O2 ⇒ (2 γ −1) CO2 + (2–2 γ ) CO γ based on Arthur, Eq. (4) C + γ O2 ⇒ (2 γ −1) CO2 + (2–2 γ ) CO γ based on Tognotti, Eq. (6) C + γ O2 ⇒ (2 γ − 1) CO2 + (2–2 γ ) CO γ based on Tognotti, Eq. (6)

◦

◦

400–500 C (TGA-char)

for the first reaction, which is the result of taking γ = 0.5. For reaction (3), both the activation energy and the pre-exponential factor are slightly reduced in comparison to the results shown in Table 5. In the derivation of the Arrhenius constants listed in Tables 5 and 6 we have assumed (see Eq. (9)) that reactions (3), (8) and (7)are first order. Almost perfect correlations in the low-temperature (kinetic) regime indicate that the reaction order remains constant over the 10–85% char conversion range at 0.21 bar oxygen pressure. Discussing the complex issue of the reaction order is beyond the scope of this paper (see for example [52–54]) and cannot be done without a systematic variation of the oxygen partial pressure. The same calculation procedure has been used for char gasification, as shown in Fig. 8. The three points are almost perfectly positioned on the line, therefore, the attained correlation coefficient is almost one. The activation temperature is 26177 K while the preexponential factor is equal to 2775 m/s. Table 7 shows literature values of activation temperature for char reactions with O2 and CO2 and they vary significantly. The differences are, among other factors, consequences of different experimental procedures used. Obtaining the kinetic constants from non-isothermal TGA experiments requires special mathematical procedures seldom used [30]. Often experiments are performed not in the intrinsic region and the obtained values of activation energy are substantially lower. The pseudo-kinetic constants are then obtained as they have been obtained in this work in the high temperature region (Table 5).

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

ARTICLE IN PRESS

JID: CNF 8

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14 Table 8 Inputs into the shrinking particle model. Parameter

Unit

GTB char

TGA char

gash gC hr (CO) hr (CO2 ) MC MCO2 MCO MO2 p SV0 Tb Tp0 Twall w

kg ash/kg char kg C/kg char kJ/kg C kJ/kg C kg C/kmol C kg CO2 /kmol CO2 kg CO/kmol CO kg O2 / kmol O2 bar m2 /m3 K K K m/s m3 gas / m3 − kg/m3 kg/m3 solid −

0.14 0.86 9.212E + 3 [62] 32.796E + 3 [62] 12 44 28 32 1 3.11E + 7 1069 300 600 3.3 0.286 0.8 1163 1628 1.2

0.12 0.88 9.212E + 3 [62] 32.796E + 3 [62] 12 44 28 32 1 1.99E + 7 1069 300 600 3.3 0.191 0.8 1257 1554 1.19

0 p ρ app 0 ρ true τ pore

Fig. 8. Arrhenius constants determination for char gasification reaction (C + CO2 ⇒ 2 CO).

4. Shrinking particle model The Shrinking particle model described in the Appendix is used in the following four versions: (a) (b) (c) (d)

γ γ γ γ

= 1 (CO2 is produced only) = 0.5 (CO is produced only) calculated using Eq. (4) as in Arthur [43] calculated using Eq. (6) as in Tognotti et al. [44].

Table 8 shows input variables used in the SPM.

4.1. Results The model has been run to simulate the experiments performed for four nearly-spherical char particles listed in Table 2 (see Fig. 3).

The initial particle diameters in the SPM are: 4.75 mm, 6.20 mm, 5.90 mm and 5.15 mm for experimental runs No. 2, 3, 4 and 5, respectively. Figure 9 shows measured (Runs no. 4 and 5) and calculated mass loss; the char conversion progresses faster when combustion of carbon to CO is considered in comparison to combustion to CO2 . This is associated with the fact that two times less oxygen is consumed per kmol of carbon if carbon monoxide is the final product. Nevertheless, two times faster char conversion is not observed in the model predictions. This is because the particle temperature is higher in the case of combustion to CO2 (Fig. 10) since the reaction enthalpy of carbon oxidation to carbon dioxide is higher (see Table 8). When CO and CO2 are the oxidation products, the combustion progresses faster than when CO is the only product.

Table 7 Literature values of activation temperature for char oxidation. Reference

Ta int , K (reaction with O2 )

[22]

8000–10, 000

[55]

4900–6400

[56]

circa 18, 600

[20]

16, 900

[33]

18, 760 and circa 13, 200

[13] [36] [57]

15400–16840 16840–18400 15, 640

[58] [59] [60]

13, 750 15, 640–16, 720 15, 560 and 16, 720

Ta int , K (reaction with CO2 )

25, 140–26, 800

26, 930

[39]

16, 960, 19, 970 and 27, 180

[61]

24, 740 and 12, 420

This work

14, 560–15, 970

26, 177

Description

Reactive coals; Assumed model value Non-isothermal TGA, Elbistan lignite Non-isothermal TGA, Semi-anthracite and Low-volatile bituminous chars Single char particle test rig, Brown coal char Non-isothermal TGA, Commercial graphite and high Volatile bituminous coal char TGA, bituminous coal chars TGA, Australian coal chars Standard value used for calculation of char thermal deactivation FBB model assumptions Non-isothermal TGA, coal chars Non-isothermal and isothermal TGA, Bituminous coal char Czech coke, Polish coke and Graphite, first two values are Apparent parameters, technical TGA Small test reactor, Polish bituminous and lignite chars Isothermal TGA in intrinsic Region, char produced from Bituminous coal “Janina”

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF

ARTICLE IN PRESS

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

9

Table 9 Relative difference of average conversion rates between model predictions and measurements; δm−mod in %. Particle number

Reaction products of SPM: CO2 and CO Eq. (4)/Eq. (6)

Reaction product of SPM: only CO2

Reaction product of SPM: only CO

2 3 4 5

43.9/38.0 65.3/58.8 35.8/30.7 50.9/45.1

9.2 5.0 13.7 4.7

38.3 58.6 31.3 45.8

Average

49.0/43.2

8.2

43.5

Fig. 9. Measured and the SPM calculated mass loss; Squares (Run no. 4) and triangles (Run no. 5) indicate measured data (initial mass and diameters shown in Table 2), while lines indicate calculation results (legend indicates the final product of carbon oxidation).

Fig. 11. Recalculated char burnout in the 20% to 95% range; squares (Run no. 4) and triangles (Run No. 5) indicate measured data (initial mass and diameters shown in Table 2), while lines indicate calculation results (legend on the top of chart indicate the final product of carbon oxidation).

Fig. 10. Measured and SPM calculated particle temperatures and burnout for Run no. 5 (initial mass and diameters shown in Table 2); Squares show measured maximal surface temperature; Diamonds show the measured burnout while lines indicate calculation results (legends indicate the product of carbon oxidation).

The particles mass losses measured for Runs no. 4 and 5 fall between the SPM results obtained for carbon oxidation to CO2 and to CO alone (Fig. 9). During the initial twenty seconds the experimental mass loss rate is larger than the model predictions. The initial experimental mass loss is not only related to carbon combustion, but also to moisture (0.015) and residual volatiles (0.077) release from the particle. At the initial stage of GTB-char combustion, the external surface is larger in comparison to the surface of the modeled sphere. The larger external surface area is associated with cracks found on the surface of manually-rounded char particles. The irregularities on the particle surface are burned first and the particle becomes more and more spherical when combustion progresses. At the initial stage of combustion, the model assumption that only carbon is the combustible part is inaccurate, since residual hydrogen is also present. All these model assumptions make the SPM not accurate for the initial stage of char combustion. Passing through the first stage of char combustion for particles No. 4 and 5, at time around 40–50 s, the mass loss rate gradually decreases (see Fig. 9). The rate decrease is comparable with mass loss rate obtained from the model which assumes carbon combustion to CO2 only. To quantify the models, the relative difference of average conversion rates (δm−mod ) have been calculated according to:

  mod  exp δm−mod = (crexp av − cr av ) /cr av

(22)

Average conversion rates (cr av = m/τ ) are calculated using particle mass at 20% and at 95% of carbon burnout. Table 9 shows clearly that the lowest relative difference (8.2% on average) of average conversion rates is found for the model where carbon dioxide is the only product of carbon oxidation. The models, in which carbon oxidizes according to C + γ O2 ⇒ (2 γ −1) CO2 + (2−2γ ) CO reaction or C + 0.5 O2 → CO reaction, overpredict the conversion rates by 40% or more (see Table 9). Marginally different mass loss has been calculated when Arthur expression (Eq. (4)) is replaced with Tognotti et al. (Eq. (6)) relationship, as shown in Figs. 9 and 10. The char conversion rates predicted by the model and the experimentally-obtained values for the 20–95% char burnout range are shown in Fig. 11. The char burnout (X∗ ) is recalculated here according to Eq. (23)

X ∗ = [m p (X = 20%) − m p (τ )]/[m p (X = 20%) − m p (X = 95%)] (23) so as to eliminate the ambiguities in the initial oxidation stages. The experimentally-obtained conversion rates for Run no. 5 (and also for Run no. 3, not shown here) are within 4.2% margin from the model predictions in which carbon dioxide is the only product of carbon oxidation. For Run no. 4 (and Run no. 2, not shown here) the model predictions underestimate the degree of burnout (X∗ ). For example at 100 s, the model predicts 0.69 burnout instead of 0.8 as observed in the experiments. As shown in Table 9, the char conversion rates obtained using the two other models, where oxidation proceeds either to CO or CO and CO2 , are too high in comparison with the experimental data. Figure 10 shows both the particle temperature obtained using the SPM and the measured values (Run no. 5). Good correlation between maximal registered temperature of the particle surface and modelcalculated particle temperature has been found from around 30 s

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF 10

ARTICLE IN PRESS

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

Fig. 12. Run 5 – Terms appearing in the energy balance-(Eq. (24)). Carbon conversion to CO2 only.

onwards when the carbon oxidation to CO2 is used in the SPM. Not only is the curve tendency preserved but the temperature level is also similar. When carbon oxidation in the SPM progresses to either CO and CO2 or to CO only, the maximum measured value of surface particle temperature is around 200 K larger than the SPM calculated temperatures. It is instructive to examine both the energy balance equation (Eq. (26) in Appendix) and its individual terms in order to comprehend the measured and model calculated temperatures shown in Fig. 10. The energy balance (Eq. (26)) can be simplistically written in the following form:

Reactions Enthalpy Release + Enthalpy Transport + Convection d Tp + Radiation = ρ · c p · Vp · (24) dt where the terms: reactions enthalpy release, enthalpy transport, convection and radiation are easily identifiable in Eq. (26) of the Appendix. These four terms together with their sum, which is called as total (net) energy rate, are shown in Fig. 12 (for carbon conversion to CO2 only) which we analyze together with Fig. 10. During the first 15 s, the particle is heated up mainly by convection. The char oxidation becomes gradually faster causing the particle temperature to reach around 1400 K. From that instant onwards the char oxidation, in the 27–98% burnout range, proceeds at almost constant temperature. The reactions enthalpy release is counterbalanced by heat losses; the total energy rate to the particle is low but always remains positive. As shown in Fig. 12 radiation is the dominant mechanism of heat removal. Both the reactions’ enthalpy release and the radiation term decrease, in absolute values, almost linearly with time. This is because, for burnout larger than around 27%, the combustion process is controlled by oxygen supply rate to the particle external surface; the pseudo-kinetic term (Eq. (31)), being in the 3.4–8.2 m/s range, is substantially larger than the mass transfer rate (Eq. (36)) which is around 0.3 m/s. In other words the d2 -law for the combustion rate holds. Since the particle temperature remains constant at the 1400 K level, the radiation term decreases also linearly with time. A careful inspection of Fig. 10 indicates that although the SPM calculated particle temperature remains at around 1400 K level, it increases very slowly with time. In other words, the total (net) energy supply rate to the particle is very low and positive (see Fig. 12), invoking that dTp /dt remains also low and stays always positive. When the char burnout approaches 100%, the SPM calculated temperature increases rapidly suggesting that a singularity occurs when approaching complete combustion. In Fig. 13 the four energy terms appearing in Eq. (24), each divided by the particle volume, are plotted for

Fig. 13. Run 5 – Energy rate per particle volume. Carbon conversion to CO2 only.

the last 2% burnout. It is apparent that, in absolute values, the reactions enthalpy term increases faster than the heat losses all together (radiation+convection+enthalpy transport). The reactions enthalpy release divided by the particle volume varies with particle diameter as 1/d2p since

Reactions Enthalpy Release (−dm/dt ) · hr ∼ Particle Volume d3p ∼

d2p · d−1 p · hr · d3p

∼

hr d2p

(25)

while the remaining terms (radiation, convection, enthalpy transport) show 1/dp behavior. Thus when the particle diameter decreases, the reaction enthalpy release term increases always faster than the remaining three (negative) terms governing the heat removal rate. The corollary of this observation is that derivative dTp /dt never gets negative for high enough reactions’ enthalpies and it rapidly increases when approaching 100% burnout. This is the SPM singularity occurring when the particle diameter approaches zero. Thus, at the last 2% burnout, the SPM calculated particle temperatures do not agree with the measured values, as shown in Fig. 10. On the other hand, at the end of the GTB-experiments there is always ash remaining on the experimental stand and the particle diameter never goes to zero. One may argue then that the SPM correctly describes the GTBexperiments with exception of the last few percent of burnout. 5. Rate controlling mechanisms In the above presented analysis, we have focused on the GTBexperiments with emphasis on experimental Runs 4 and 5 which, like all other experiments (Table 2), have been carried out in oxidizer (air) containing 21% oxygen. With exception of the initial 20 s of the experiments, where heating up and pseudo-kinetics are rate controlling, the char combustion rate remains controlled by oxygen transfer rate. Thus, a question is whether all our consideration concerning the intrinsic kinetics and inter-particle diffusion is needed or, in other words, whether calculations based on mass transfer only would not be sufficient. Fig. 14 shows the SPM model predictions of Run 5 obtained by ignoring the pseudo-kinetic term in Eq. (30). It is apparent that the mass loss curve for such predictions runs almost parallel to the full (kinetics and mass transfer) SPM calculations. The difference is obviously in the offset of char combustion. In the full SPM model predictions the char combustion is initiated after reaching a particle temperature which provides appreciable rates of pseudo-kinetics. In the SPM predictions, where the oxygen diffusion is considered only, the char combustion proceeds immediately with appreciable rates.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF

ARTICLE IN PRESS J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

Fig. 14. The char mass loss calculations for 21% (Run no. 5), 10%, 5% and 2% oxygen mole fraction. Dashed lines show predictions obtained using the full (oxygen diffusion + pseudo kinetics) SPM while solid lines show SPM predictions obtained with kinetics switched off (diffusion only). Carbon conversion to CO2 only.

This deficiency could be easily fixed by introducing a threshold temperature (for Run no. 5 it would be around 1000K) for the char oxidation initiation and then the mass loss curve would be indeed very similar to the full SPM model predictions. Thus, for a mathematical description of experimental Runs 2,3,4 and 5, a SPM model, taking into account only oxygen transport to the particle surface, would work well. Figure 14 shows also the SPM predictions for 10%, 5% and 2% oxygen mole fractions. For each oxygen mole fraction two predictions are presented; one obtained using the full SPM (oxygen diffusion plus pseudo-kinetics) and the other predictions obtained using the SPM with kinetics switch-off (oxygen diffusion only). Substantial differences are visible and it is obvious that the kinetics cannot be ignored; the lower the oxygen concentration the more important the kinetics. The particle temperature is determined by the enthalpy release rate and at lower oxygen concentrations the temperatures at the plateau level are 1220K, 1000K and 910K for 10%, 5% and 2% O2 , respectively. At such temperatures the pseudo-kinetic rates are becoming comparable to the oxygen diffusion rate. Thus, a SPM model, which is to handle variable oxygen mole fractions (lower than 21%), has to have both mechanism incorporated. We wish to complete this Section with a remark concerning SPM model predictions for experimental Runs no. 4 and 5 in which Arthur [43] and Tognotti et al. [44] expressions for CO/CO 2 ratios have been used. As shown in Figs. 9 and 10 as well as in Fig. 11 both options provide almost identical results. For particle temperatures exceeding 900 K, Tognotti expression gives a slightly larger γ than Arthur formula, as shown in Fig. 6. Thus for the same char oxidation rate (-dm/dt), Tognotti at al. formulation gives just a slightly larger reactions enthalpy. Such a difference results in a very small alteration to the overall combustion rate since the process is controlled by oxygen diffusion to the particle surface and pseudo-kinetics plays a role, for a short period of time, in the initial stages only. Consequently, the mass loss and the particle temperatures are indeed very similar regardless whether Arthur or Tognotti expressions are used. 6. Oxidation rates of the GTB-char and the TGA-char Two chars have been made as described above; the GTB-char generated at a high heating (devolatilization) rate and the TGA-char at a low heating rate. The measured intrinsic reactivity for both chars is listed in Table 5 while their morphology data is given in Tables 3 and 4. The ratio of SV0 for the GTB-char and the TGA-char reaches 1.64 and it partially explains the difference in the TGA observed combustion

[m5G;October 14, 2015;22:25] 11

Fig. 15. Recalculated char burnout in the 20% to 95% range; triangles (Run no. 4) show the measured data for the GTB-char (initial mass and diameters shown in Table 2) while squares indicate measured data for the 5.0 mmTGA-char particle (not listed in Table 2); lines indicate calculation results (legend indicate the final product of carbon oxidation).

Fig. 16. The SPM calculated mass loss obtained for two initial particle mass (Run no. 4 – blue line; Run no. 5 – black line; see Table 2); Solid lines indicate results obtained for GTB-char, dotted lines for TGA-char, while dashed lines show models predictions for TGA-char, in which GTB-char morphology is implemented; final product of carbon oxidation is CO2 only. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

rates (see Section 2.3), as shown in Figs. 5 and 4. Thermal annealing is also likely to slow down the conversion rates of TGA-char in comparison with GTB-char (see Fig. 5). It is then instructive to examine the oxidation rates of both chars using the GTB rig. Figure 15 shows the measured char burnout (X∗ ) for the 5.15 mm GTB-char particle (Run no. 5 in Table 2) and for the 5.0 mm TGA-char particle (not listed in Table 2) and it is apparent that the TGA-char burns slower than the GTB-char. The SPM predictions, obtained for γ = 1 and the model inputs as listed in Table 8, also show the same behavior. To deconvolute the influence of morphology from the annealing effect, the SPM has been employed to calculate conversion rates of the TGA-char using the TGA-char intrinsic reactivity and the GTBchar morphology (see Tables 3 and 4). New Arrhenius constants for the TGA-char have been calculated, following the procedure described in Section 3, and they are only slightly different from the values listed in Table 5 for 400–500 °C range. The subsequent new ∗ pre-exponential factor (ATintGA = 1035 m/s) and actication tempera∗ T GA ture (Ta int = 15339 K) have been entered into the SPM together with the morphology data of the GTB-char (SV0 ,  , ρ app0 , ρ true , τ pore as listed in Table 8) and these simulations are marked in Fig. 16 as TGAchar with GTB-char morphology.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF 12

ARTICLE IN PRESS

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

The zero-dimensional Shrinking Particle Model (SPM) has been developed to reproduce, in the simplest possible way, the experimentally observed combustion rates at 21% oxygen mole fraction. Three versions of the model have been examined. In the first one, carbon is oxidized directly to CO2 , in the second one to CO while in the third option the oxidation produces a CO/CO2 mixture in ratios determined by the expression of either Arthur [43] or Tognotti et al. [44]. Only the first option provides the predictions which are in agreement with the measured data. In general the SPM model predictions agree well with the measured data with exception of the last 2% burnout where excessive particle temperature is calculated due to the model singularity occurring when the particle diameter goes to zero. Acknowledgments

Fig. 17. The SPM calculated particle temperature obtained for Run no. 5 measured maximal particle surface temperature; Solid lines indicate results obtained for GTBchar, dotted lines for TGA-char, dashed lines show models predictions for TGA-char, in which GTB-char morphology is implemented, while squares indicate measured data; final product of carbon oxidation is CO2 only.

Figure 16 shows the SPM predicted particle mass loss for the GTBchar, the TGA-char, and for TGA-char in which GTB-char morphology is implemented. Initial particle masses are as for Runs no. 4 and 5 (see Table 2) and GTB-char mass loss calculations are the same as presented in Fig. 9. The mass loss rate predicted by the SPM for TGAchar is lower than for GTB-char. The difference in combustion time reaches around 75 s and 58 s for Runs no. 4 and 5, respectively. The SPM predictions for the case where both chars have the same morphology are positioned between calculations of the TGA-char particle and the GTB-char particle, as expected. Eq. (22) is used to quantify these differences; the difference between average conversion rates obtained for the GTB-char and the TGA-char is 21.8%, while between calculations obtained for GTB-char and TGA-char with morphology from GTB-char is 13.1%. This means that annealing is responsible for around 60% of the slow-down while morphology for 40%. Figure 17 shows the char particle temperature calculated for the considered cases. The temperature for the TGA-char particle combustion (dotted lines) as well as for the TGA-char particle combustion in which morphology is taken as for the GTB-char, is around 50 K lower than the measured maximal value or calculated for the GTBchar. The lower temperature is related to the lower conversion rates for the TGA-char particle.

Authors of this work would like to thank Anna Pajdak of the Insti˛ tute of Advanced Energy Technologies, the Czestochowa University of Technology, who carried out the porosimetry measurements. Jakub Bibrzycki is a beneficiary of the project “DoktoRIS – Scholarship program for innovative Silesia” co-financed by the European Union from the European Social Fund. This work has been carried out with financial assistance within the frame of the statute research of the Institute of Thermal Technology, the Silesian University of Technology. Appendix. The Shrinking Particle Model Figure 18 shows the energy balance diagram for the SPM. The energy conservation equation (Eq. 26), which is used to calculate the particle temperature, includes convective and radiative heat transfer, energy generation inside the particle due to heterogeneous char oxidation reactions as well as enthalpy transport due to species transport.

(−dm/dt ) · hr  T p  298.15 − hCO2 − hCO2 +

 

hTCOp

 298.15

− hCO

· (2 − 2 · γ ) · (−dm/dt ) · (MCO /MC )

 298.15

− hTOb2 2 − hO2

· (2 · γ − 1) · (−dm/dt ) · (MCO2 /MC )

· (−dm/dt ) · γ · (MO2 /MC )



4 + α · Sext · (Tp − Tb ) +  p−pw · σ · Sext · Tp4 − Twall

=



dTp · (Vp · (1 − 0 ) · cchar · ρtrue + Vp · 0 · c p air (Tp ) · ρair (Tp , p)) dt (26)

The specific heat capacity of char is calculated as: 7. Conclusions This work presents a study on combustion of millimeter-size char particles. The newly developed Grain Thermo-Balance (GTB) has been used to measure overall combustion rates of 4.8–6.3 mm char particles in air. The measured data is used to develop a zero-dimensional mathematical model which is to serve as a char particle sub-model in a CFD-based software for predicting performance of fluidized bed boilers. To generate the char samples, coal particles have been devolatilized at 10 K/min heating rate in the TGA and at around 590 K/min rate in the GTB. At 21% oxygen mole fraction and at temperatures lower than 773 K both chars burn in the kinetic (intrinsic) regime and the GTB-char reactivity is around four times larger than that for the TGA-char. At temperatures larger than 1023 K both chars burn in the diffusion regime with comparable overall combustion rates. It has been found that the TGA-char has lower porosity and internal specific surface area in comparison to the GTB-char. These differences in the particle morphology are responsible for around 40% slow-down of the combustion rate, if compared to the GTB-char, while annealing is responsible for the remaining 60%.

cchar = gash · cash + gC · cC

(27)

Fig. 18. Diagram of energy balance of the char particle used in the SPM.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

ARTICLE IN PRESS

JID: CNF

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

13

where the specific heat of carbon (cc ) and enthalpies of CO2 , CO and O2 have been calculated using JANAF tables while cash is estimated using the following equation [63,64]:

where the Nusselt number is calculated as [13,68]:

cash = (594 + 0.586 · Tp )/1000

In order to calculate the particle mass and temperature at an instant, Eqs. (26) and (29) have to be integrated numerically using the initial particle mass and temperature as initial conditions; mass and temperature are then independent variables. A more elegant way is to normalize these variables to the initial particle mass and temperature, respectively. The latter approach allows for an easy control of the integration error. However, this is just a convenience. Since the equations are not stiff, the integration is rather straight forward. We have performed the integration using both the primitive variables (mass and temperature) and the normalized ones obtaining identical results. As an integration method we have used the 5th order Runge– Kutta method with an adaptive step control [70]. The integration accuracy has been set to 10−8 .

(28)

Taking external and pore diffusion as well as heterogeneous reactions into consideration, the char conversion rate can be calculated as:

(−dm/dt ) = Vp · MC · (1/γ ) · η · kint · SV0 · βeff · Sext · CO2 ·(βeff · Sext + η · kint · SV0 · Vp )−1

(29)

The above relationship can be recast into:

(−dm/dt ) =

MC · (1/γ ) · Sext · CO2 1 kint ·(η·SV 0 /(Sext /Vp ))

+ β1 ef f

(30)

with the pseudo-kinetics (intrinsic kinetics plus pore diffusion) defined as

kpseudo = kint ·

η · SV 0 Sext /Vp

(31)

The effectiveness is the ratio of the internal surface area which participates in heterogeneous reactions, to the internal surface area available for reactions in entire particle volume. This parameter is an indicator of pore diffusion limitation of the conversion rate [36] and is calculated using Thiele modulus (φ ) as [65,66]:

η = (3/φ 2 ) · (φ /tanh(φ) − 1)

(32)

while the modulus is evaluated as [67]:

φ = (d p /2) · (SV 0 · kint /De f f )0.5

(33)

The effective diffusivity coefficient is calculated as:

De f f = (0 /τ pore ) · (1/Dg + 1/Dk )−1

(34)

where the gas diffusivity is evaluated as for a binary O2 -air system at bulk (oxidizer) temperature. The Knudsen diffusivity, describing the diffusion rates inside pores, is calculated as:

Dk = (2/3) · (d pore av /2) · ((8 · R · Tp )/(π · MO2 ))0.5

(35)

while the effective mass transfer coefficient (β eff ), used in Eq. (29), is [38]:

βe f f = Sh · Dgb /d p

(36)

where the Sherwood number [13,68]: Sh = 2 + 0.6 · Reb 0.5 · Sc1/3

(37)

and the Schmidt number are evaluated for the bulk temperature according to [13,69]: Sc = μb /Dgb · ρb

(38)

with the Reynolds number being: Reb = w · d p · ρb /μb

(39)

In this SPM it is assumed that the particle shrinks immediately when the portion of carbon is consumed so that the actual particle diameter is related to burnout through:

d p = d p0 · (1 − X )1/3

(40)

where dp0 is the initial particle diameter. The reaction enthalpy of carbon oxidation, which appears in Eq. (26), is:

hr = hr (CO) · (2 − 2 · γ ) + hr (CO2 ) · (2 · γ − 1)

(41)

Since the particle is small, if compared to the surrounding walls, the mutual emissivity of the particle-wall system is taken to be equal to the particle emissivity ( p−pw =  p ). The convective heat transfer coefficient is estimated using:

α = Nu · λb /(d p )

(42)

Nu = 2 + 0.6 · Reb 0.5 · Pr0.33 b

(43)

References [1] K. Myohanen, T. Hyppanen, Int. J. Chem. Reactor Eng. 9 (2011) 1–55. [2] R.I. Singh, A. Brink, M. Hupa, Appl. Thermal Eng. 52 (2013) 585–614. [3] W. Adamczyk, G. Wecel, M. Klajny, P. Kozolub, A. Klimanek, R. Bialecki, Particuology 16 (2014) 29–40. [4] M.A. Field, D. Gill, B. Morgan, P. Hawksley, Combusti. Pulverized Coal, BCURA, 1967. [5] D. Liu, X. Chen, W. Thou, C. Zhao, Proc. Combust. Inst. 33 (2011) 2701–2708. [6] G.A. Simons, Proc. Combust. Inst. 19 (1982) 1067–1076. [7] N.E.L. Haugen, M.B. Tilghman, R.E. Mitchell, Combust. Flame 161 (2) (2014) 612– 619. [8] N.M. Laurendeau, Progress Energy Combust. Sci. 4 (1978) 221–270. [9] I.W. Smith, Proc. Combust. Inst. 19 (1982) 1045–1065. [10] R. Weber, M. Mancini, Z. Phys. Chem. 229 (5) (2015) 619–641. ´ ´ ˛ A. Malon, [11] R. Bonda, D. Brzezinski, M. Czapigo-Czapla, G. Czapkowski, J. Dylag, S.Z. Mikulski, W. Mi´skiewicz, S. Oszczepalski, A. Piotrowska, D. Siekierka, ´ L. Skrzypczyk, J. Sokołowski, W. Szczygielski, M. Szuflicki, M. Tyminski, A. Wałkuska, Bilans zasobów złóz˙ kopalin w Polsce wg stanu na 31 XII 2012 ´ ´ r., Technical Report, Panstwowy Instytut Geologiczny - Panstwowy Instytut Badawczy, 2013. [12] V. Manovic, M. Komatina, S. Oka, Fuel 87 (2008) 905–914. [13] I. Guedea, D. Pallares, L.I. Diez, F. Johnsson, Fuel Process. Technol. 112 (2013) 118– 128. [14] B. Remiarova, J. Markos, R. Zajdlik, L. Jelemensky, Fuel Process. Technol. 85 (2004) 303–321. [15] A.K. Sadhukhan, P. Gupta, R.K. Saha, Fuel Process. Technol. 90 (2009) 692–700. [16] K. Holikova, R. Zajdlik, J. Markos, L. Jelemensky, Chem. Papers 59(6a) (2005) 413– 420. [17] A.N. Hayhurst, M.S. Parmar, Combust. Flame 130 (2002) 361–375. [18] F. Scala, Energy Fuels 25 (2011) 1051–1059. [19] F. Scala, R. Chirone, Chem. Eng. J. 165 (2010) 902–906. [20] R. Zajdlik, L. Jelemensky, B. Remiarova, J. Markos, Chem. Eng. Sci. 56 (2001) 1355– 1361. [21] M. Komatina, V. Manovic, D. Dakic, Energy Fuels 20 (2006) 114–119. [22] K. Svoboda, M. Hartman, M. Pohorely, O. Tranka, Acta Geodyn. Geomater. 1 (2004) 261–274. [23] M.J. Biggs, P.K. Agarwa, Chem. Eng. Sci. 52 (1997) 941–952. [24] A. Gomez-Barea, P. Ollero, R. Arjona, Fuel 84 (2005) 1695–1704. [25] T.H. Fletcher, L.L. Baxter, D.K. Ottesen, Spectral emission characteristics of sizegraded coal particles, Technical Report, Sandia National Laboratories, 1987. [26] H. Gao, J. H., M. Nomura, Energy Fuels 24 (2010) 18–28. [27] C.A. Heidenreich, H.M. Yan, D.K. Zhang, Fuel 78 (1999) 557–566. [28] F. Winter, M.E. Prach, H. Hofbauer, Combust. Flame 108 (1997) 302–314. [29] J. Bibrzycki, Investigations of coal particle combustion and gasification Ph.D. thesis, Silesian University of Technology and Clausthal University of Technology, 2014. [30] R. Weber, J. Energy Inst. 81 (2008) 226–233. [31] P. Ollero, A. Serrera, R. Arjona, S. Alcantarilla, Fuel 81 (2002) 1989–2000. [32] B. Nowak, O. Karlstrom, P. Backman, A. Brink, M. Zevenhoven, S. Voglsam, F. Winter, M. Hupa, J. Thermal Anal. Calorim. 111 (2012) 183–192. [33] A. Zolin, A.D. Jensen, P.A. Jensen, K. Dam-Johansen, Fuel 81 (2002) 1065–1075. [34] L. Gasparovic, Z. Korenova, L. Jelemensky, Chem. Papers 64 (2010) 174–181. [35] A. Gomez-Barea, P. Ollero, C. Fernandez-Baco, Energy Fuels 20 (2006) 2202–2210. [36] D.G. Roberts, D.J. Harris, Energy Fuels 14 (2000) 483–489. [37] G. Hakvoort, J.C. Schouten, P.J.M. Valkenburg, J. Thermal Anal. 35 (1989) 335–346. [38] S. Brunello, I. Flour, P. Maissa, B. Bruyet, Fuel 75 (1996) 536–544. [39] S. Ulzama, A Theoretical Analysis of Single Coal Particle Behavior during Spontaneous Devolatilization and Combustion Ph.D. thesis, Otto von Guericke University Magdeburg, 2007. [40] D.G. Roberts, E.M. Hodge, D.J. Harris, J.F. Stubington, Energy Fuels 24 (2010) 5300– 5308. [41] G. Chen, Q. Yu, K. Sjostrom, J. Anal. Appl. Pyrolysis 40-41 (1997) 491–499.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024

JID: CNF 14

ARTICLE IN PRESS

[m5G;October 14, 2015;22:25]

J. Bibrzycki et al. / Combustion and Flame 000 (2015) 1–14

[42] O. Karlstrom, A. Brink, E. Biagini, M. Hupa, L. Tognotti, Proc. Combust. Inst. 34 (2013) 2427–2434. [43] J. Arthur, Trans. Faraday Soc. 47 (1951) 164–178. [44] L. Tognotti, J. Longwell, A. Sarofim, Proc. Combust. Inst. (1990) 1207–1213. [45] C.R. Shaddix, F. Holtzleithner, M. Geier, B.S. Haynes, Combust. Flame 160 (2013) 1827–1834. [46] C. Chen, T. Kojima, Fuel Process. Technol. 47 (1996) 215–232. [47] H.R. Batchelder, R.M. Busche, W.P. Armstrong, Ind. Eng. Chem. 45 (1953) 1856– 1878. [48] A. Jess, A.-K. Andresen, Fuel 89 (2010) 1541–1548. [49] S.K. Bhatia, D.D. Perlmutter, AIChE J. 26 (1980) 379–386. [50] B. Feng, S.K. Bhatia, Carbon 41 (2003) 507–523. [51] P.A. Webb, An Introduction to the Physical Characterization of Materials by Mercury Intrusion Porosimetry with Emphasis on Reduction and Presentation of Experimental Data, Technical Report, Micromeritics Instrument Corp., 2001. [52] R.H. Hurt, J.M. Calo, Combust. Flame 125 (2001) 1138–2145. [53] R.H. Hurt, B. Haynes, Proc. Combust. Inst. 30 (2005) 2161–2168. [54] M. Geier, C. Shaddix, F. Holtzleither, Proc. Combust. Inst. 34 (2013) 1411–1418. [55] H. Sis, J. Thermal Anal. Calorim. 88 (2007) 863–870. [56] B. Arias, C. Pevida, F. Rubiera, J.J. Pis, J. Thermal Anal. Calorim. 90 (2007) 859–863. [57] N.V. Russell, J.R. Gibbins, C.K. Man, J. Williamson, Energy Fuels 14 (2000) 883–888.

[58] [59] [60] [61] [62]

[63] [64] [65] [66] [67] [68] [69] [70]

F. Chejne, J.P. Hernandez, Fuel 81 (2002) 1687–1702. R.C. Everson, H.W.J.P. Neomagus, D. Njapha, Fuel 85 (2006) 418–422. E. Sima-Ella, G. Yuan, T. Mays, Fuel 84 (2005) 1920–1925. A. Mianowski, Z. Robak, M. Tomaszewicz, S. Stelmach, J. Thermal Anal. Calorim. 110 (2012) 93–102. H.M. Tolvanen, L.I. Kokko, R. Raiko, Challenges in Combustion Technology today, The National Committees of Sweden and Finland of the International Flame Research Foundation and the Scandinavian-Nordic Section of the Combustion Institute, 2011, pp. 1–14. Swedish-Finnish Flame Days, Pitea, Sweden D. Merrick, Fuel 62 (5) (1983) 540–546. B. Le´sniak, L. Słupik, G. Jakubina, Chemik 67 (2013) 560–571. L. Kelebopile, R. Sun, H. Wang, X. Zhang, S. Wu, Fuel Process. Technol. 111 (2013) 42–54. E.W. Thiele, Ind. Eng. Chem. Res. 31 (1939) 916–920. B. Peters, A. Dziugys, R. Navakas, Mechanika 18 (2012) 177–185. W.E. Ranz, W.R. Marshall, Chem. Eng. Progress 48 (1952) 141–146. F.P. Incropera, D.P. Dewitt, T.L. Bergman, A.S. Lavine, Fundamentals of Heat and Mass Transfer, Sixth edition, John Wiley & Sons, 2007. W.H. Press, S.A. Teukolsky, W.T. Vetting, B.P. Flannery, Numerical Receipes in C. The Art of Scientific Computing., Cambridge University Press, 1992.

Please cite this article as: J. Bibrzycki et al., A char combustion sub-model for CFD-predictions of fluidized bed combustion - experiments and mathematical modeling, Combustion and Flame (2015), http://dx.doi.org/10.1016/j.combustflame.2015.09.024