Computational Study on Co-Combustion of Coal and Blast Furnace Gas in an Industrial Fluidized Bed Boiler

Fossil fuels, which are formed over millions of years from organic matter, continue to dominate global energy consumption. However, rapid industrialization and the growing demand for energy have accelerated the depletion of these resources. This situation has created an urgent need for efficient fuel utilization strategies and the development of alternative energy approaches. Among these approaches, fuel co-combustion has emerged as an effective technique for improving fuel flexibility while reducing dependence on conventional fossil fuels. Biomass-based granular solid fuels, such as sawdust, wood chips, and recycled wood waste, are increasingly being utilized for thermal energy production due to their high energy density and wide availability. However, their combustion can lead to environmental concerns, including air pollution and greenhouse gas emissions [1, 2].

Blast Furnace Gas (BFG), a low-calorific-value byproduct generated during the ironmaking process, mainly consists of CO, CO₂, and N₂, with smaller fractions of H₂ and CH₄. The efficient utilization of BFG, particularly through co-combustion with coal, has attracted significant attention due to its potential to improve overall energy efficiency in industrial applications. Previous investigations have demonstrated the feasibility of using BFG in reheating furnaces with advanced burner configurations, enabling the achievement of high operating temperatures without the need for supplementary fuels [3]. Additionally, combustion characteristics and flame stability of BFG-based fuel mixtures have been investigated under industrial conditions, highlighting the influence of hydrogen content and operational parameters on combustion dynamics [4].

Fluidized bed combustion (FBC) is widely recognized as a suitable technology for co-combustion due to its excellent mixing, fuel adaptability, and lower emissions. Experimental studies have provided insights into particle dynamics, circulation behavior, and energy efficiency in fluidized beds [5-8]. However, such experiments are often expensive, time-consuming, and limited in capturing detailed flow and reaction phenomena [9, 10].

To overcome these limitations, Computational Fluid Dynamics (CFD) has emerged as a powerful tool for analyzing multiphase reacting systems. Both Eulerian–Lagrangian and Eulerian–Eulerian approaches have been extensively applied in fluidized bed simulations [11-15]. Comparative studies indicate that while the Eulerian–Lagrangian method provides detailed particle-scale resolution, the Eulerian–Eulerian approach offers better computational efficiency and scalability for industrial systems, with comparable accuracy in predicting hydrodynamic behavior [16, 17]. Furthermore, incorporation of the Kinetic Theory of Granular Flow (KTGF) enables accurate representation of solid-phase interactions and bed dynamics [18]. Advanced CFD studies have also successfully modeled complex reacting systems such as fluid catalytic cracking regenerators and chemical looping combustion processes [4, 5].

Despite these advancements, most existing research is limited to single-fuel combustion systems or laboratory-scale configurations. Comprehensive studies addressing the co-combustion of coal and BFG in industrial scale circulating fluidized bed (CFB) boilers, particularly with validation against real plant data, remain scarce.

To address this gap, the present study investigates a 67 TPH circulating fluidized bed combustion (CFBC) boiler operating in an integrated steel plant. A three-dimensional CFD model based on the Eulerian–Eulerian approach is developed and validated using industrial measurements, including pressure, temperature, and phase distribution along the boiler height. Following validation, the model is applied to analyze the co-combustion of coal and BFG, focusing on temperature distribution, species evolution, emission characteristics, and the interaction of homogeneous and heterogeneous reactions.

The outcomes of this study provide both scientific and practical significance by offering a validated computational framework for industrial-scale co-combustion systems. The proposed approach can support design optimization, improve fuel utilization, and reduce dependence on costly experimental investigations, thereby contributing to more sustainable energy practices in the steel and power sectors.

3.1 Geometry creation

The computational domain of the CFBC boiler includes the distributor plate, coal feed inlets, secondary air-ports, and the riser section as depicted in Figure 1. At the base of the reactor, the distributor plate is equipped with 300 circular nozzles, each with a diameter of 100 mm, ensuring uniform air distribution into the riser. Coal is supplied through two feed inlets mounted on the side walls, each having a cross-sectional area of 500 mm × 500 mm. In addition, seven secondary air-ports measuring 200 mm × 200 mm are positioned along the front and rear riser walls to enhance mixing and combustion. The riser geometry is divided into two distinct regions: a lower section that gradually diverges to a height of 4000 mm and an upper straight section extending 15,000 mm. Under bubbling fluidized bed operating conditions, only a small fraction of solid particles reach the cyclone zone. For this reason, the cyclone section is excluded from the numerical model used in the present study.

Figure 1. Computational domain

Table 1. Grid independency test

3.2 Discretization

The domain is divided into a total of 1,457,482 hexahedral cells, a count established through a rigorous grid independence test performed across a range of cell counts (482,642, 792,365, 1,046,862, 1,457,482, and 1,757,832), as shown in Table 1. During this test, results showed convergence and stability once the cell count reached 1,457,482, indicating a sufficient level of accuracy for capturing the flow and combustion characteristics within the CFBC boiler. While 1,457,482 cells were determined to be optimal for computational efficiency without compromising accuracy, additional calculations were performed at 1,457,482 cells. The discretised domain is depicted in Figure 2.

Figure 2. Meshing

The findings from this denser grid are detailed in the results and discussion section to provide a thorough validation and ensure the robustness of the model outcomes.

3.3 Governing equation

Mass balance equations for the gas and solid phases:

$\frac{\partial \left( {{\alpha }_{g}}{{\rho }_{g}} \right)}{\partial t}+\nabla \left( {{\alpha }_{g}}{{\rho }_{g}}{{v}_{g}} \right)=0$    (1)

$\frac{\partial \left( {{\alpha }_{c}}{{\rho }_{c}} \right)}{\partial t}+\nabla \left( {{\alpha }_{c}}{{\rho }_{c}}{{v}_{c}} \right)=0$    (2)

$\frac{\partial \left( {{\alpha }_{s}}{{\rho }_{s}} \right)}{\partial t}+\nabla \left( {{\alpha }_{s}}{{\rho }_{s}}{{v}_{s}} \right)=0$    (3)

Momentum balance equations for the gas and particle phases:

$\begin{aligned} & \frac{\partial\left(\alpha_g \rho_g v_g\right)}{\partial t}+\nabla\left(\alpha_g \rho_g v_g v_g\right) = -\alpha_g \nabla P+\nabla \tau_g+\alpha_g \rho_g \vec{g}+F_{\text {ext }}\end{aligned}$    (4)

$\begin{gathered}\frac{\partial\left(\alpha_c \rho_c v_c\right)}{\partial t}+\nabla\left(\alpha_c \rho_c v_c v_c\right)=-\alpha_c \nabla P-\nabla P_c+\nabla\vec \tau_c+\alpha_c \rho_c \vec{g}-F_{e x t}\end{gathered}$    (5)

$\begin{gathered}\frac{\partial\left(\alpha_s \rho_s v_s\right)}{\partial t}+\nabla\left(\alpha_s \rho_s v_s v_s\right)=-\alpha_s \nabla P-\nabla P_s+\nabla \vec\tau_s+\alpha_s \rho_s \vec{g}-F_{e x t}\end{gathered}$    (6)

In the Eulerian-Eulerian granular model ${{P}_{c}}$ and ${{P}_{s}}~$is calculated as follows [19]:

${{P}_{c}}=\text{ }\!\!~\!\!\text{ }{{\alpha }_{c}}{{\rho }_{c}}{{\theta }_{c}}+2{{\rho }_{c}}\left( 1+{{e}_{cc}} \right)\alpha _{c}^{2}{{g}_{o}}{{\theta }_{c}}$    (7)

${{P}_{s}}=\text{ }\!\!~\!\!\text{ }{{\alpha }_{s}}{{\rho }_{s}}{{\theta }_{s}}+2{{\rho }_{s}}\left( 1+{{e}_{ss}} \right)\alpha _{s}^{2}{{g}_{o}}{{\theta }_{s}}$     (8)

${{g}_{o}}={{\left[ 1-{{\left( \frac{{{\alpha }_{s}}}{{{\alpha }_{{{s}_{max}}}}} \right)}^{{}^{1}/{}_{3}}} \right]}^{-1}}$   (9)

and, ${{\tau }_{g}},$ ${{\tau }_{c}}$ and ${{\tau }_{s}}$ are the stress tensor for the gas, coal, and sand phases.

$\tau_g=\alpha_g \mu_g\left(\nabla v_g+\nabla v_g^t\right)-\frac{2}{3} \alpha_g \mu_g\left(\nabla \cdot v_g\right) \overline{\bar{I}}_g$    (10)

${{\tau }_{c}}={{\alpha }_{c}}{{\mu }_{c}}\left( \nabla {{v}_{c}}+\nabla v_{c}^{t} \right)-{{\alpha }_{c}}\left( {{\varepsilon }_{c}}-\frac{2}{3}{{\mu }_{c}} \right).{{\vec{v}}_{c}}$    (11)

${{\tau }_{s}}={{\alpha }_{s}}{{\mu }_{s}}\left( \nabla {{v}_{s}}+\nabla v_{s}^{t} \right)-{{\alpha }_{s}}\left( {{\varepsilon }_{s}}-\frac{2}{3}{{\mu }_{s}} \right).{{\vec{v}}_{s}}$    (12)

Mathematical formulae for granular temperature and solid viscosity, particularly for the coal phase, are shown in the next section. We will also apply an analogous set of equations to the sand phase:

${{\mu }_{c}}={{\mu }_{Collisional}}+{{\mu }_{kin}}+{{\mu }_{friction}}$

Ding and Gidaspow [20] and Lun [21] provide the formulae for shear, kinetic, and collisional viscosities, which are given in terms of particle diameter and the restitution coefficient.

3.3.1 Collisional viscosity

${{\mu }_{Col}}=\frac{4}{5}{{\alpha }_{c}}{{\rho }_{c}}{{d}_{c}}{{g}_{o}}\left( 1+{{e}_{cc}} \right){{\left( \frac{{{\theta }_{c}}}{\pi } \right)}^{{}^{1}/{}_{2}}}$    (13)

Viscosity due to fluctuation in kinetic energy:

${{\mu }_{kin}}=\frac{{{\alpha }_{c}}{{d}_{c}}{{\rho }_{c}}\sqrt{{{\theta }_{c}}\pi }}{6\left( 3-{{e}_{cc}} \right)}\left[ 1+\frac{2}{5}\left( 1+{{e}_{cc}} \right)\left( 3{{e}_{cc}}-1 \right){{\alpha }_{c}}{{g}_{o}} \right]$     (14)

3.3.2 Frictional viscosity

${{\mu }_{f}}=\frac{{{P}_{c}}sin\varnothing }{2\sqrt{{{I}_{2D}}}}$    (15)

The mixture technique is used to describe turbulence interactions between the solid and gas phases. By examining the mixture's combined characteristics and velocities, this model illustrates important aspects of turbulent flow. Below is the model's governing equation for the turbulent mixture.

$\begin{aligned} & \frac{\partial}{\partial t}\left(\rho_m k\right)+\nabla \cdot\left(\rho_m v_m k\right)= \nabla \cdot\left(\frac{\mu_{t, m}}{\sigma_k} \nabla k\right)+G_{k, m}-\rho_m \varepsilon\end{aligned}$    (16)

$\begin{gathered}\frac{\partial}{\partial t}\left(\rho_m \varepsilon\right)+\nabla \cdot\left(\rho_m v_m \varepsilon\right)=\nabla \cdot\left(\frac{\mu_{t, m}}{\sigma_{\varepsilon}} \nabla \varepsilon\right)+\frac{\varepsilon}{k}\left(C_{1 e} G_{k, m}-C_{2 e} \rho_m \varepsilon\right)\end{gathered}$   (17)

${{\rho }_{m}}=\mathop{\sum }_{i=1}^{N}{{\alpha }_{i}}{{\rho }_{i}}$   (18)

${{\vec{v}}_{m}}=\frac{\mathop{\sum }_{i=1}^{N}{{\alpha }_{i}}{{\rho }_{i}}{{v}_{i}}}{\mathop{\sum }_{i=1}^{N}{{\alpha }_{i}}{{\rho }_{i}}}$    (19)

where, ${{\mu }_{t,m}}~$is turbulent viscosity, calculated as:

${{\mu }_{t,m}}={{\rho }_{m}}{{C}_{\mu }}\frac{{{k}^{2}}}{\varepsilon }$    (20)

where, ${{C}_{\mu }}$ = 0.09 and ${{G}_{k,m}}$is turbulent kinetic energy production given below:

${{G}_{k,m}}={{\mu }_{t,m}}\left( \nabla {{v}_{m}}+{{\left( \nabla {{v}_{m}} \right)}^{T}} \right):\nabla {{v}_{m}}$    (21)

KTGF is applied to represent the fluctuating motion of particles and their associated kinetic energy. Due to the random movement of particles, their fluctuating energy component is expressed through a parameter commonly referred to as the granular temperature [22, 23]. The governing transport equation derived from kinetic theory is presented below:

$\begin{gathered}\frac{3}{2}\left[\frac{\partial}{\partial t}\left(\alpha_c \rho_c \theta_c\right)+\nabla \cdot\left(\alpha_c \rho_c v_c \theta_c\right)\right]=\left(-P_c \cdot \overline{\bar{I}}+\tau_c\right) \cdot \nabla v_c-\gamma_{\theta_c}+\nabla \cdot\left(k_{\theta_c} \cdot \nabla \theta_c\right) +3 K_{g c} \theta_c\end{gathered}$   (22)

Within this formulation, the first term on the right-hand side represents energy production caused by the solid stress tensor. The second term corresponds to energy diffusion, while the third accounts for collisional dissipation. The final contribution describes interphase energy exchange. Each term in the equation is expressed as a function of granular temperature, solids pressure, solids viscosity, and bulk viscosity. The coefficient ${{k}_{\theta c}}$ denotes the diffusive constant for the coal phase and is defined in Eq. (20):

${{k}_{{{\theta }_{c}}}}=\frac{150\text{ }\!\!~\!\!\text{ }{{\rho }_{c}}{{d}_{c}}\sqrt{{{\theta }_{c}}\pi }}{384\left( 1+{{e}_{cc}} \right){{g}_{o}}}{{\left[ 1+\frac{6}{5}{{\alpha }_{c}}{{g}_{o}}\left( 1+{{e}_{cc}} \right) \right]}^{2}}+2\alpha _{c}^{2}{{\rho }_{c}}{{d}_{c}}\left( 1+{{e}_{cc}} \right){{g}_{o}}\sqrt{\frac{{{\theta }_{c}}}{\pi }}$      (23)

${{\gamma }_{{{\theta }_{c}}}}=\frac{12{{\left( 1-{{e}_{cc}} \right)}^{2}}{{g}_{o}}}{{{d}_{c}}\sqrt{\pi }}\alpha _{c}^{2}{{\rho }_{c}}\theta _{c}^{{}^{3}/{}_{2}}$     (24)

For highly dense FB, the convection and diffusion variables may be disregarded since the algebraic equation, which is the second version of the granular temperature equation, shows that the granular energy generation and dissipation rates are in equilibrium [24].

$0=\left(-P_c \cdot \overline{\bar{I}}+\tau_c\right) \cdot \nabla v_c-\gamma_{\theta_c}+3 K_{g c} \theta_c$    (25)

3.3.3 Interfacial force characterization

Interactions between gas and solid phases involve several interfacial forces, including drag, lift, and virtual mass effects. Among these, the drag force between coal particles and the gas phase is the primary focus of this section. The sand phase is modeled using an equivalent set of governing relations.

${{F}_{ext}}={{F}_{drag}}+{{F}_{lift}}+{{F}_{vm}}$

${{F}_{drag}}={{K}_{cg}}\left( {{v}_{g}}-{{v}_{c}} \right)$

Numerous empirical and semi-empirical expressions for the drag coefficient have been reported in the literature for gas–solid fluidized systems. The correlations most frequently applied include the Gidaspow model [25], the Wen and Yu formulation [26], and the Syamlal–O’Brien model. The predicted distributions of gas and solids volume fractions, as well as overall flow structures, are strongly dependent on the selected drag closure. For bubbling fluidized beds, the Syamlal–O’Brien correlation is widely recognized as one of the most reliable options and is therefore adopted in the present work. The corresponding expression for the drag interaction is given below:

${{K}_{cg}}=\frac{3{{\alpha }_{c}}{{\alpha }_{g}}{{\rho }_{g}}}{4v_{rc}^{2}{{d}_{c}}}{{C}_{d}}\left( \frac{{{R}_{ec}}}{{{v}_{rc}}} \right)\left( {{v}_{g}}-{{v}_{c}} \right)$    (26)

The interphase momentum exchange coefficient is adjusted using a terminal settling velocity relationship, where ${{v}_{rc}}~$represents the ratio of a particle’s terminal settling velocity to that of an equivalent spherical particle.

${{v}_{rc}}=0.5\left( A-0.06R{{e}_{c}}+\sqrt{{{\left( 0.06R{{e}_{c}} \right)}^{2}}+0.12R{{e}_{c}}\left( 2B-A \right)+{{A}^{2}}} \right)$     (27)

$A=\alpha _{g}^{4.14}$

$B=0.8\alpha _{g}^{1.28}$ For ${{\alpha }_{g}}$ < 0.85

$B=\alpha _{g}^{2.65}$ For ${{\alpha }_{g}}$ > 0.85

${{C}_{d}}=\left( 0.63+\frac{4.8}{\sqrt{\frac{R{{e}_{s}}}{{{v}_{rc}}}}} \right)$     (28)

$R{{e}_{c}}=\frac{{{\rho }_{g}}\left( {{v}_{g}}-{{v}_{c}} \right){{d}_{c}}}{{{\text{ }\!\!\mu\!\!\text{ }}_{g}}}$

3.3.4 Formulation for species mass fraction transport

A finite-rate/eddy-dissipation formulation is adopted to describe the kinetics of homogeneous reactions. In this approach, the overall reaction rate is determined by selecting the smaller value predicted by either the chemical finite-rate expression or the turbulence-controlled eddy dissipation mechanism [27, 28].

${{\hat{R}}_{i,r}}=\left( v_{i,r}^{\prime\prime }-v_{i,r}^{'} \right)\left( {{k}_{fr}}\mathop{\prod }_{j=1}^{N}{{\left[ {{c}_{j,r}} \right]}^{\left( n_{i,r}^{'}+n_{i,r}^{\prime\prime } \right)}} \right)$     (29)

The stoichiometric coefficients of reactants and products are represented by $v_{ij}^{\text{ }\!\!'\!\!\text{ }}~$and $v_{ij}^{\text{ }\!\!'\!\!\text{  }\!\!'\!\!\text{ }}$, respectively. Here, ${{k}_{fr}}~$denotes the forward reaction rate constant, and ${{C}_{j}}~$corresponds to the molar concentration of species $j$. The parameter $N~$indicates the total number of chemical species involved in the reaction system. The subsequent relation defines the reaction rate contribution associated with turbulent eddy dissipation.

${{\hat{R}}_{i,r}}=v_{i,r}^{'}{{M}_{w,i}}AB\rho \frac{\varepsilon }{k}\frac{\mathop{\sum }_{p}{{Y}_{p}}}{\mathop{\sum }_{j}^{N}v_{j,r}^{N}{{M}_{w,i}}}$      (30)

The turbulence-generated kinetic energy is represented by $k$, while its dissipation rate is indicated by $\varepsilon $. Unless otherwise specified, the empirical model constants take default values of $A=4$ and $B=0.5$. The molecular weight and mass fraction of the ${{i}^{th}}$chemical species are expressed as ${{M}_{w,i}}$ and ${{Y}_{i}}$, respectively.

3.4 Rate equation for solid gas reaction system

Coal is composed primarily of ash, moisture, fixed carbon, and volatile matter. In contrast, BFG is a gaseous fuel mixture containing roughly 25% CO, 20% CO₂, 3% H₂, and 52% N₂ [6]. When introduced into the boiler, coal and BFG undergo a sequence of both heterogeneous and homogeneous reactions with the surrounding gas environment. Because these fuels differ substantially in their physical characteristics and chemical behavior, separate kinetic parameters are required to accurately represent their combustion processes.

The literature has plenty of information on the kinetic parameters and reaction rates for coal's heterogeneous processes, including char oxidation, tar oxidation, moisture evaporation, volatile release, and char interactions with CO₂, H₂O, and NO [29, 30]. In this study, relevant kinetic parameters from various sources have been used to model these coal reactions. C programming is used for creating the intrinsic reaction rates for the heterogeneous processes, which are then integrated into the Fluent solver via a user-defined function (UDF). Below provides a summary of the kinetics and stoichiometry of the heterogeneous reactions between various phases.

3.4.1 Coal moisture release

${{\left( {{H}_{2}}O \right)}_{Coal}}\to {{\left( {{H}_{2}}O \right)}_{Vapour}}$

$R_{m r}=A_{\text {drying }} \times \exp \left(-\frac{E_{\text {drying }}}{R_{\text {gas }} \cdot T_c}\right) \alpha_c \rho_c Y_{\left(H_2 O\right)_{\text {coal }}}$    (R1)

Here, ${{R}_{gas}}\text{ }\!\!~\!\!\text{ }$represents the universal gas constant, and ${{T}_{c}}$ denotes the coal particle temperature. The term ${{Y}_{{{H}_{2}}O,coal}}\text{ }\!\!~\!\!\text{ }$corresponds to the moisture mass fraction within the coal. The parameters ${{A}_{drying}}\text{ }\!\!~\!\!\text{ }$and ${{E}_{drying}}$ are the pre-exponential factor and activation energy associated with the drying process of the coal phase.

3.4.2 Coal devolatilization

During coal devolatilization, various gaseous hydrocarbons are released, including CH₄, H₂O (vapor), CO₂, CO, H₂S, NH₃, and tar. The MGAS model [31] serves as the basis for the devolatilization rate in this investigation. Mass balance calculations are used to calculate the stoichiometric coefficients for the devolatilization reaction as well as the molecular formulae for tar and volatile matter.

                                        $Volatile\to {{\alpha }_{1}}Tar+{{\alpha }_{2}}CO+{{\alpha }_{3}}C{{O}_{2}}+{{\alpha }_{4}}C{{H}_{4}}+{{\alpha }_{5}}{{H}_{2}}O+{{\alpha }_{6}}{{H}_{2}}S+{{\alpha }_{7}}N{{H}_{3}}+{{\alpha }_{8}}{{H}_{2}}$

$R_d=\left\{\begin{array}{cc}A_{\text {cod_devoletilization}} \times \exp \left(\frac{E_{\text {cod_devoletilization}}}{R_{gas} T_c}\right)\left(1-\alpha_c\right) \rho_c\left(Y_{VM_{\text {cod }}}-Y^*\right) & Y_{M_{\text {cod }}} \geq Y^* \\ 0 & Y_{M_{\text {cod }}}<Y^*\end{array}\right.$    (R2)

${{Y}^{*}}=\frac{{{\rho }_{{{c}_{o}}}}\left( {{Y}_{F{{C}_{o}}}}+{{Y}_{V{{M}_{o}}}} \right)Y_{o}^{*}}{{{\rho }_{c}}}$

$Y_o^*=\left\{\begin{array}{cc}\left(\frac{867.2}{\left(T_c-273\right) / 100}\right)^{3.914} & T_c<1223 \\ 0 & T_c<1223\end{array}\right.$

Here, ${{Y}_{VM}}$ represents the instantaneous mass fraction of volatile matter remaining in the coal. The parameters ${{Y}_{FC,0}}~$and ${{Y}_{VM,0}}$ denote the initial mass fractions of fixed carbon and volatile matter, respectively. The quantity ${{Y}^{\text{*}}}$ defines the remaining fraction of volatile material in the solid phase.

3.4.3 Tar oxidation

When coal devolatilizes, tar—a viscous, black, and sticky liquid—is created. It is made up of intricate chemical molecules that include a lot of carbon. Tar oxidizes when it meets oxygen, producing CO₂ and H₂O. The tar oxidation kinetics and stoichiometry used in this investigation are taken from the study [29].

$Tar+13.557{{O}_{2}}\to 10.79C{{O}_{2}}+5.715{{H}_{2}}O$

$\begin{aligned} R_{\text {tar_oxidation }}= & A_{\text {tar_oxidation }} \times \exp \left(-\frac{E_{\text {tar_oxidation }}}{R_{\text {gas }} \cdot T_c}\right) \times \alpha_g \times\left(\frac{\rho_g Y_{O_2}}{M_{W, O_2}}\right)^{1.5} \times\left(\frac{\rho_g Y_{\text {tar }}}{M_{W, O_2}}\right)^{0.25}\end{aligned}$   (R3)

3.4.4 Char oxidation

Char refers to the solid residue remaining in the furnace after the devolatilization stage of coal conversion. It is primarily composed of carbon, along with smaller quantities of hydrogen, nitrogen, oxygen, and sulphur. Although the detailed mechanisms governing char oxidation are complex and involve multiple interacting processes, the present study adopts a simplified modelling framework [32]. In this representation, char is treated mainly as a carbon-rich solid, while other species released to the gas phase during devolatilization, including hydrogen, nitrogen, and oxygen, are accounted for separately.

Char oxidation is controlled by three dominant resistances: intrinsic surface reaction kinetics, diffusion through the surrounding gas film, and diffusion within the ash layer. The oxidation model proposed in the study [30] is implemented in the current work.

$Cha{{r}_{coal}}+0.511{{O}_{2}}\to CO+0.0075NO+0.008S{{O}_{2}}$

$R_{\text {char_oxidation }}=\frac{6 \alpha_c p_{\mathrm{O}_2}}{d_c\left(\frac{1}{k_{\text {film }}}+\frac{1}{k_{\text {ash }}}+\frac{1}{k_r}\right) \times M_{W, O_2}}$    (R4)

${{k}_{film}}=\frac{{{D}_{{{O}_{2}}}}Sh}{{{d}_{c}}\times \frac{R}{{{M}_{W,{{O}_{2}}}}}\times {{T}_{g}}}$

${{D}_{{{O}_{2}}}}=4.26{{\left( \frac{{{T}_{g}}}{1800} \right)}^{1.75}}$

$Sh=\left( 7-10{{\alpha }_{g}}+5\alpha _{g}^{2} \right)\times \left( 1+0.7\operatorname{Re}_{c}^{0.2}S{{c}^{1/3}} \right)+\left( 1.33-2.4{{\alpha }_{g}}+1.2\alpha _{g}^{2} \right)\times \operatorname{Re}_{c}^{0.2}S{{c}^{1/3}}$

${{k}_{ash}}=\frac{2{{r}_{d}}{{D}_{e}}}{\left( 1-{{r}_{d}} \right){{d}_{c}}\times \frac{R}{{{M}_{W,{{O}_{2}}}}}\times {{T}_{c}}}$

The current and original carbon mass fractions in the coal particle are used to determine ${{r}_{d}}$, defined as the ratio between the shrinking char core diameter and the overall particle diameter. The ash layer’s effective diffusivity, denoted by ${{D}_{a}}$, is expressed as follows:

$D_e=D_{O_2}\left(0.25+0.75\left(1-Y_{A_o}\right)\right)^{2.5}$

$k_r=8710 \times \exp \left(\frac{-27000}{R_{g a s} T_c}\right) r_d^2$

3.4.5 Char reaction with CO₂

By interacting with CO₂, char undergoes gasification, producing CO, NO, and SO₂. The study's char gasification reaction kinetics with CO₂ were taken from the study [33].

$Cha{{r}_{coal}}+1.023C{{O}_{2}}\to CO+0.008NO+0.0075S{{O}_{2}}$

${{R}_{char\_c{{o}_{2}}\_gasification}}=4.4{{T}_{c}}\times \exp \left( \frac{-1.77\times {{10}^{4}}}{{{T}_{c}}} \right){{C}_{C{{O}_{2}}}}$   (R5)

3.4.6 Char reaction with H₂O

In addition to oxidation, char reacts with steam to generate gaseous products such as H₂, CO, NO, and SO₂. The gasification behavior of char in the presence of water vapor is modeled using the reaction framework reported in the study [34].

$Cha{{r}_{coal}}+1.023{{H}_{2}}O\to 1.023{{H}_{2}}+CO+0.008NO+0.0075S{{O}_{2}}$

${{R}_{char\_{{H}_{2}}O\_gasification}}=1.33{{T}_{c}}\times \exp \left( \frac{-1.77\times {{10}^{4}}}{{{T}_{c}}} \right){{C}_{{{H}_{2}}O}}$   (R6)

3.4.7 Homogeneous reaction kinetics

Various gas-phase species—such as CO, CH₄, H₂, H₂S, and NH₃—are released through heterogeneous reactions. These species then participate in homogeneous oxidation reactions with O₂. The set of homogeneous reactions incorporated in the present study is listed in Table 2, with rate parameters adopted from the study [34].

Table 2. Homogeneous reactions kinetics

3.4.8 Predictive modelling of emission formation

Modelling the development of pollutants is intrinsically difficult, because trace chemical species are involved. The CFD simulation employs a number of simplifications to improve computational performance. During devolatilization, sulphur and nitrogen are released primarily as H₂S and NH₃. These species then react with O₂ through gas-phase homogeneous reactions to form SO₂ and NO, as presented in Table 3. Additionally, the char residue, which retains some sulphur and nitrogen, further contributes to SO₂ and NO formation during its oxidation. For this study, it is assumed that NO represents the total NOₓ emissions. Table 3 lists the pertinent homogeneous reactions that contribute to the creation and degradation of contaminants.

$\begin{aligned} & \text { Char }_{\text {coal }}+1.015 \mathrm{NO} \ \quad \rightarrow \text { CO }+0.511 N_2+0.0075 S_2 \\ & \quad R_{\text {char_NO_gasification }} \quad=6 \times \alpha_c \rho_s Y_{\text {charcoal }} \alpha_g \rho_g Y_{N O} \times 1.3\quad \times10^5 \exp \left(\frac{-1.77 \times 10^4}{T_g}\right)\left(\rho_c d_p\right)\end{aligned}$   (R10)

Table 3. Kinetic parameters for the homogeneous gas-phase reactions

3.5 Boundary conditions

To accurately capture hydrodynamic and combustion behavior, appropriate initial and boundary conditions are established. Preheated primary air enters through bottom nozzles under a uniform mass flow rate, while a coal and BFG mixture is injected from side wall ports. Secondary air is introduced via seven preheated inlets on each side wall. The gas mixture's domain walls are subjected to no-slip conditions, while the solid phases' Specularity coefficient remains constant. The domain exit is defined by a pressure outlet. An initial sand bed occupies the region immediately above the distributor plate, filling the first 700 mm with a solids volume fraction of 0.48. The model considers three Eulerian phases: the gas phase, coal particles acting as fuel, and sand serving as the inert bed material, within an Eulerian–Eulerian multiphase framework. Coal composition is detailed in Table 4, while species involved and mathematical inputs are listed in Tables 5 and 6 respectively. The simulation uses a pressure-based segregated solver with the SIMPLE algorithm for pressure-velocity coupling, second-order upwind discretization, and a transient time step of 0.0001 s with 80 inner iterations for convergence.

Table 4. Characterization of coal via proximate and ultimate techniques

Table 5. Identification of Eulerian phases with respective species

Table 6. Simulation-specific mathematical inputs

Three-dimensional CFD simulations are performed on an industrial-scale CFBC boiler operating in an integrated steel facility to investigate hydrodynamic behavior and combustion performance under two conditions: pure coal firing and co-firing with a blend of 10% BFG and 90% coal. The system is modeled as a multiphase flow consisting of four interacting Eulerian phases — sand, coal, and the gas mixture — using an Eulerian–Eulerian framework coupled with the KTGF. Field measurements collected from plant operation, including temperature profiles, pressure distribution, bed height, suspension density, and void fraction, are employed to validate the numerical model. After validation, the same CFD approach is applied to evaluate boiler performance under blended fuel operation. The principal findings derived from the hydrodynamic and combustion analysis of the coal–BFG mixture are summarized below.

• For both coal and blended fuels, a greater pressure drop is seen close to the fluidized bed's bottom and diminishes with height. When compared to pure coal, blended fuel exhibits a marginally smaller pressure drop (8%, 12% reduction), suggesting better fluid movement. In industrial processes, this decrease in pressure drop can improve fluidized bed efficiency and reduce energy usage.
• A relatively high sand volume fraction is observed under both pure coal firing and blended fuel operation. Increasing the concentration of bed material near the wall region can enhance heat transfer between the hot flue gases and the furnace walls. However, excessive accumulation of solids in this region may promote particle agglomeration. Therefore, maintaining an optimal sand concentration near the walls is important to ensure stable and trouble-free operation.
• The circulation pattern of solids and gas phases within the CFBC boiler promotes strong intermixing, with sand particles rising through the central core and descending along the near-wall region. This behavior, supported by axial velocity and vector field results, enhances internal recirculation. Such continuous solids circulation is critical for achieving uniform and efficient fuel combustion inside the reactor.
• When burning mixed fuels, the mass proportion of SO₂ and NO decreases. Decreasing SO₂ and NO emissions improves the air quality and supplements the CFBC boiler in achieving and fulfilling environmental rules and emission restrictions. Higher SO₂ emissions induce corrosion in the flue gas route and downstream components. Reduction in the SO₂ mass fraction by using BFG mixture fuel has the potential to increase the longevity of components and mitigate expenses associated with maintenance of the components.
• Using the Eulerian-Eulerian multiphase model, the simulation effectively creates a thorough CFBC coal combustion model by considering coal, sand, and gas in various stages. A comprehensive physics picture of a large-scale CFBC boiler is provided by mathematical modelling, which is also useful for modelling the combustion of coal and other biofuels.