A Comprehensive Three-Dimensional Analysis of a Large-Scale Multi-Fuel CFB Boiler Burning Coal and Syngas. Part 1. The CFD Model of a Large-Scale Multi-Fuel CFB Combustion

Abstract

The paper is focused on the idea of multi-fuel combustion in a large-scale circulating fluidized bed (CFB) boiler. The article discusses the concept of simultaneous coal and syngas combustion. A comprehensive three-dimensional computational fluid dynamics (CFD) model is developed, which allows us to describe complex phenomena that occur in the combustion chamber of the CFB boiler burning coal and syngas produced from coal sludge.

1. Introduction

Awareness of the growing environmental pollution and climate change is leading to the emergence of new methods to reduce emissions from fossil fuel combustion [1,2,3]. One of the most favored steam generation technology in recent times is circulating fluidized bed (CFB) technology. The main advantages of CFB boilers are stable operation, low acid gas emissions, and fuel flexibility with the potential for reliable operation with difficult-to-burn fuels [4,5,6].

Since the complexity of fuel combustion in the CFB boilers, especially in the large-scale units, is still not sufficiently recognized, the model research of such objects operation is of practical significance. Different modeling approaches to modeling CFB facilities can be distinguished, including artificial intelligence (AI) and the programmed computing approach [7,8,9,10]. The AI approach constitutes the use of the methods, which can reproduce a process from training samples [11,12,13,14]. The programmed computing approach base on writing algorithms to solve elaborate mathematical descriptions of the considered processes [5,6,15,16,17,18,19,20,21]. Basu et al. [5,6] provided a wide-range review and comparison of circulating fluidized bed combustor CFBC models. A similar qualification shows Gungor and Eskin et al. [22]. The authors selected three groups of CFB models, corresponding to three levels of sophistication. The Group I encompasses the 1D, mostly, plug flow/stirred tanks or axial solids distribution models, using simplified mass and energy balance. Group II contains 2D (1.5D), core/annulus models with a broad consideration of combustion and other related processes, Finally, Group III constitute 3D models, mainly CFD (computational fluid dynamics) models.

Models built using the CFD approach are the most general and numerically sophisticated, advanced, based on detailed consideration of chemical kinetics and individual physical processes, gas, and solids continuity equations, momentum balances, and appropriate constitutive equations [23]. Many various 3D models of CFB boilers exist in literature.

An interesting application of a hybrid Euler–Lagrange approach to model the dense gas–solid flow combined with a combustion process in a large-scale industrial CFB boiler was developed in [24]. Modifications of the existing code via the user-defined function UDF mechanism as well as the selection of a proper combination of submodels, numerical techniques, and the careful control of the solution convergence implemented in the ANSYS Fluent allowed the successful validation of the developed model against the measured data sets [24].

A model for the simulation of reactive gas–solids flows in large, industrial CFB combustors is developed in [25]. The 105 MW_e Duisburg and the 235 MW_e Turów combustors are discussed at work. A semiempirical approach was employed to describe the 3D combustion of coal in [26]. The calculations were performed on the 12 MW_th CFB boiler of the Chalmers University of Technology is modeled. A three-dimensional numerical model has was successfully established to simulate olive cake combustion in a CFB combustor in [27]. The calculations were conducted for a laboratory-scale CFD model unit.

A comprehensive CFD combustion model for a large-scale supercritical CFB boiler was presented in [28]. The authors considered a large-scale 350 MW supercritical CFB boiler.

The fuel flexibility of CFB technology has caused a discussion on the feasibility of different fuels co-combustion in a CFB furnace. Thus the concept of multi-fuel CFB units has emerged. As a general rule, the kind of fuel dedicated to an individual boiler depends on various factors, including fuel properties and design parameters [29]. Since the CFB boilers allow using wide size distribution feedstock [30], many works are devoted to biomass and coal co-combustion [31,32,33], considering that biomass belongs to the environmentally friendly, renewable fuels. However, the authors underlined that the use of biofuels tends to operational problems during the combustion process due to, e.g., defluidization as the effect of the bed sintering, superheater fouling, and high-temperature corrosion. The main reason for such biofuels’ behavior is their high content of sodium and potassium, leading to lowering the ash softening point [4,15,34]. Two combustion approaches with wood-based biomass, coal, peat, and solid recovered fuel (SRF) utilization were discussed in [29]. The authors underlined the feasibility of using 100% SRF, and 100% demolition wood.

The co-combustion results of sewage sludge with coal and wood in laboratory-scale and pilot-scale 12 MW_th CFB boilers were considered in [35].

The interaction of fuels in co-firing in FBC was also thoroughly discussed by Huppa [36]. The author emphasized that some fuel mixtures show surprising interactions as there are no sufficient data about the behavior of the fuels when they are burned in the mixtures. Since the change in operational conditions occurs while the fuel type changes, it is difficult to draw out the process behavior. The author gives a review of some research on the fuel interaction in fluidized bed combustors. The cases with a gas co-combustion were also listed but not discussed [36].

The use of gaseous fuel, with a properly organized gas supply system, can be beneficial for power unit flexibility. A key issue is to provide the minimum technically justified load to avoid the instability of combustion [37]. However, an essential advantage of converting solid fuel into process gas is enabling stable operation of the boiler furnace with a substantial reduction of capacity. Additional positive aspects of the proposed concept are economic, social, and ecological benefits from the utilization of coal sludge [33,35,36].

On the other hand, the acceptable parameters may be exceeded, which may result in the destruction or, at best, reduction of the service life of individual components. Such effects may be a consequence of a change in the flue gas velocity and temperature profiles [4,18,37,38,39].

The above literature review reveals that the latest bibliography lacks such comprehensive 3D-CFD modeling of a large-scale, multi-fuel, CFB boiler, burning coal, and syngas.

In this paper, the multi-fuel coal and syngas combustion process in the large-scale OFz-425 CFB boiler installed in Poland was simulated. The commercial CFD ANSYS Workbench package, including the ANSYS Fluent Solver, was employed to describe flow and thermal conditions in the furnace of the CFB boiler.

To the best of our knowledge, the present work is the first one in the open in the literature, such as a comprehensive 3D CFD model of a large-scale multi-fuel CFB boiler burning coal and syngas.

2. Description of the OFZ-425 CFB Boiler

The simulations are conducted for the OFz-425 large-scale CFB boiler, produced by RAFAKO S.A., Poland. It is a two-pass CFB unit burning bituminous coal combustion (Figure 1). The main parts of the boiler are the riser with the combustion chamber made of membrane-walls, superheater II (SH II), and reheater II (RH II), the flue gas discharge with cyclones, and the second pass in the convection cage consisting of built-in superheaters and reheaters as well as an economizer and a rotary air heater.

The technical data of the boiler are listed in

Table 1. Technical data of the large-scale OFz-425 CFB boiler.

Boiler Type	OFz-425
Nominal thermal capacity, MWth	335.8
Live steam temperature, °C	560
Live steam pressure, MPa	16.1
Secondary steam temperature, °C	560
Secondary steam pressure, MPa	4.0/3.8
Feedwater temperature, °C	250
The efficiency of the boiler, %	91
Maximum continuous rating, t/h	425
Live steam temperature, °C	560
Secondary steam flow rate, t/h	371.2–382.0
Feedwater pressure, MPa	17.9
Excess air in the furnace,	1.2
The temperature of the fluidized bed, °C	850–870
Fluidizing gas velocity, m/s	5–8

.

Bituminous coal is the primary fuel dedicated to the CFB boiler. Its properties are listed in

Table 2. Properties of coal (as received).

Lower heating value (LHV), MJ/kg	16.7
Proximate analysis ar, wt %
Fixed Carbon (by difference) FC	21.1
Volatile, V	36.6
Ash, A	20.6
Moisture, W	21.7
Ultimate analysis ad, wt %
Carbon, C	44.5
Hydrogen, H	3.62
Oxygen, O	12.89
Nitrogen, N	1.17
Sulfur, S	1.8

.

One of the options for processing coal sludge is its gasification in gasification reactors and then the use of the obtained syngas as an auxiliary fuel in power boilers [2,40]. The feasibility of coal and syngas co-combustion in the OFz-425 boiler is considered during the study. One of the syngas production technologies is the use of gasification technology proposed by Synthesis Energy Systems, Inc. (SES), Huston USA. During the gasification mixture of 70% coal sludge (8% moisture) and 30% fine coal, the gas product with the properties given in

Table 3. Properties of syngas.

Lower heating value (LHV), MJ/kg	6.09
Temperature, °C	128
Pressure, bar	9.6
Mole mass, kg/kmol	21.45
Density, kg/m3	6.228
Viscosity, Pas	1.84 × 10−6
Specific heat, kJ/(kg K)	1.623
Thermal conductivity, W/(m K)	0.0568
Composition, mol%
H	24.70
CO	18.84
CO2	17.55
COS	0.01
CH4	3.39
H2O	26.13
H2S	0.28
NH3	0.15
N	8.44
Ar	0.51

is obtained.

The amount of batch sand is equal to 64,500 kg, of which 50,000 kg fills the furnace chamber, and the rest is contained in the seals of the boiler.

3. Methods

The defined research problem concerns the analysis of CFD modeling of the multi-fuel OFz-425 CFB boiler supplied by bituminous coal and syngas from coal sludge. The purpose of the model research is to recognize thermal-flow conditions in the CFB boiler in a nominal operating configuration and under the modified mode, i.e., with the supply of the auxiliary syngas.

ANSYS software was used during the model research, as the world’s leading simulation package, enabling comprehensive calculations in almost every field of science and industry. The mathematical model of the OFz 425 CFB boiler base on the algorithms of numerical flow simulations with chemical reactions that occur in the furnace. The two-phase model of homogenous and heterogeneous combustion employs the computational mechanics of CFD fluids and transport equations for the flow of the reactive mixture. The equations for mass, enthalpy, momentum, and selected gaseous components relevant to the combustion process are applied using the finite volume method.

The practical purpose of the study is to determine the increase in temperatures inside the combustion chamber as a result of the syngas supply and to compare the results with the critical operating temperatures of the fluidized bed. The mass balance and solid-phase transport DPM (discrete phase modeling) equations are considered in the model.

The Eddy dissipation model is employed with the radiation heat transfer equations, according to the discrete ordinates model.

In order to accurately simulate flow phenomena in non-laminar flows, Reynolds averaged Navier–Stokes equations (RANS) are supplemented with the k-omega BSL turbulence model, based on the Boussinesq hypothesis.

The equations are solved in the ANSYS Fluent solver based on the finite volume method. A detailed description of the physical models used in the paper and the methods of numerical solutions of the equations are described in the consecutive subsections.

3.1. Transport Equations with Combustion (Eddy Dissipation Method)

The mass conservation equation (the continuity equation) can be written as (1):

$$\frac{\partial \mathsf{\rho }}{\partial \mathrm{t}}+\nabla ·\left(\mathsf{\rho }\overline{\mathrm{v}}\right)={\mathrm{S}}_{\mathrm{m}}$$

(1)

Equation (1) is a general form of mass conservation equation for both compressible and incompressible flows, where the S_m is the mass source.

Equation (1) for 2D axial-symmetric issues can be written, as follows:

$$\frac{\partial \mathsf{\rho }}{\partial \mathrm{t}}+\frac{\partial }{\partial \mathrm{x}}\left({\mathsf{\rho }\mathrm{v}}_{\mathrm{x}}\right)+\frac{\partial }{\partial \mathrm{r}}\left({\mathsf{\rho }\mathrm{v}}_{\mathrm{r}}\right)+\frac{{\mathsf{\rho }\mathrm{v}}_{\mathrm{r}}}{\mathrm{r}}={\mathrm{S}}_{\mathrm{m}}$$

(2)

In an inertial system (without acceleration), the momentum conservation equation can be described by the Equation (3):

$$\frac{\partial }{\partial \mathrm{t}}\left(\mathsf{\rho }\overline{\mathrm{v}}\right)+\nabla ·\left(\mathsf{\rho }\overline{\mathrm{v}}\overline{\mathrm{v}}\right)=-\nabla \mathrm{p}+\nabla ·(\overline{\overline{\mathsf{\tau }}})+\mathsf{\rho }\overline{\mathrm{g}}+\overline{\mathrm{F}}$$

(3)

where: p is static pressure, $\overline{\overline{\mathsf{\tau }}}$ is the stress tensor (described below), $\mathsf{\rho }\overline{\mathrm{g}}$ is the gravity, and $\overline{\mathrm{F}}$ expresses the external mass forces and other source terms present in the model.

The stress tensor mentioned above ($\overline{\overline{\mathsf{\tau }}}$) is defined as (4):

$$\overline{\overline{\mathsf{\tau }}}=\mathsf{\mu }\left[\left(\nabla \overline{\mathrm{v}}+\nabla {\overline{\mathrm{v}}}^{\mathrm{T}}\right)-\frac{2}{3}\nabla \overline{\mathrm{v}}\mathrm{I}\right]$$

(4)

where μ is the dynamic viscosity, and $\mathrm{I}$ stands for the unit tensor.

Three mechanisms of heat transport, i.e., conduction, convection, and radiation, can be considered in the energy equation of the following form (5):

$$\frac{\partial }{\partial \mathrm{t}}\left(\mathsf{\rho }\mathrm{E}\right)+\nabla ·\left(\overline{\mathrm{v}}\left(\mathsf{\rho }\mathrm{E}+\mathrm{p}\right)\right)=\nabla ·\left({\mathrm{k}}_{\mathrm{eff}}\nabla \mathrm{T}-{\displaystyle \sum }_{\mathrm{j}}{\mathrm{h}}_{\mathrm{j}}{\overline{\mathrm{J}}}_{\mathrm{j}}+({\overline{\overline{\mathsf{\tau }}}}_{\mathrm{eff}}·\overline{\mathrm{v}})\right)+{\mathrm{S}}_{\mathrm{h}}$$

(5)

where ${\mathrm{k}}_{\mathrm{eff}}$ is an effective conduction coefficient, ${\overline{\mathrm{J}}}_{\mathrm{j}}$ the mass flow rate of a diffusive component $\mathrm{j}$, the next terms of Equation (5) define energy transfer by conduction, diffusion of substances, viscous dissipation, and the term ${\mathrm{S}}_{\mathrm{h}}$ is a user-defined term, representing volume heat source.

$\mathrm{E}$ from Equation (5) is defined as (6):

$$\mathrm{E}=\mathrm{h}-\frac{\mathrm{p}}{\mathsf{\rho }}+\frac{{\mathrm{v}}^2}{2}$$

(6)

where $\mathrm{h}$-enthalpy is expressed for ideal gases as:

$$\mathrm{h}={\displaystyle \sum }_{\mathrm{j}}{\mathrm{Y}}_{\mathrm{j}}{\mathrm{h}}_{\mathrm{j}}$$

(7)

and for incompressible flows (8):

$$\mathrm{h}={\displaystyle \sum }_{\mathrm{j}}{\mathrm{Y}}_{\mathrm{j}}{\mathrm{h}}_{\mathrm{j}}+\frac{\mathrm{p}}{\mathsf{\rho }}$$

(8)

The term ${\mathrm{Y}}_{\mathrm{j}}$ in Equations (7) and (8) stands for the mass fraction of component j and ${\mathrm{h}}_{\mathrm{j}}$ is expressed by Equation (9):

$${\mathrm{h}}_{\mathrm{j}}=\underset{{\mathrm{T}}_{\mathrm{ref}}}{\overset{\mathrm{T}}{{\displaystyle \int }}}{\mathrm{c}}_{\mathrm{p},\mathrm{j}}\mathrm{dT}$$

(9)

The reference temperature ${\mathrm{T}}_{\mathrm{ref}}$, employed to calculate the enthalpy, depends on the solver and the models used.

In the combustion model, without premixing, the energy equation in the form of the total enthalpy Equation (10) can be written:

$$\frac{\partial }{\partial \mathrm{t}}\left(\mathsf{\rho }\mathrm{H}\right)+\nabla ·\left(\mathsf{\rho }\overline{\mathrm{v}}\mathrm{H}\right)=\nabla ·\left(\frac{{\mathrm{k}}_{\mathrm{eff}}}{{\mathrm{C}}_{\mathrm{p}}}\nabla \mathrm{H}\right)+{\mathrm{S}}_{\mathrm{h}}$$

(10)

Assuming that the Lewis number (Le) is equal to 1, the conduction and the diffusion term of the substance combine forming the first term in Equation (10), and the viscous dissipation is included in the second term of the above equation. Total $\mathrm{H}$ is defined as follows:

$$\mathrm{H}={\displaystyle \sum }_{\mathrm{j}}{\mathrm{Y}}_{\mathrm{j}}{\mathrm{H}}_{\mathrm{j}}$$

(11)

where ${\mathrm{Y}}_{\mathrm{j}}$ is the mass fraction of the component $\mathrm{j}$ and

$${\mathrm{H}}_{\mathrm{j}}=\underset{{\mathrm{T}}_{\mathrm{ref},\mathrm{j}}}{\overset{\mathrm{T}}{{\displaystyle \int }}}{\mathrm{c}}_{\mathrm{p},\mathrm{j}}\mathrm{dT}+{\mathrm{h}}_{\mathrm{j}}^0\left({\mathrm{T}}_{\mathrm{ref},\mathrm{j}}\right)$$

(12)

where ${\mathrm{h}}_{\mathrm{j}}^0\left({\mathrm{T}}_{\mathrm{ref},\mathrm{j}}\right)$ is the enthalpy of the substance i formation at the reference temperature.

When heat is supplied to a fluid, its density changes with temperature, generating convection. The effect of buoyancy forces on mixed convection can be determined by the ratio of Grashof and Reynolds numbers (3.13):

$$\frac{\mathrm{Gr}}{{\mathrm{Re}}^2}=\frac{\mathrm{g}\mathsf{\beta }\mathsf{\Delta }\mathrm{TL}}{{\mathsf{\nu }}^2}$$

(13)

For the pure form of natural convection, the importance of buoyancy-induced flow can be determined by the Rayleigh number (14):

$$\mathrm{Ra}=\frac{{\mathrm{g}\mathsf{\beta }\mathsf{\Delta }\mathrm{TL}}^3\mathsf{\rho }}{\mathsf{\mu }\mathsf{\alpha }}$$

(14)

where: $\mathsf{\beta }$ is the thermal expansion coefficient (15):

$$\mathsf{\beta }=-\frac{1}{\mathsf{\rho }}{\left(\frac{\partial \mathsf{\rho }}{\partial \mathrm{T}}\right)}_{\mathrm{p}}$$

(15)

The value of the thermal diffusivity factor is expressed as (16):

$$\mathsf{\alpha }=\frac{\mathrm{k}}{{\mathsf{\rho }\mathrm{C}}_{\mathrm{p}}}$$

(16)

When modeling of combustion processes, the transport and mixing components, as well as chemical reactions, should also be considered. For this purpose, ANSYS Fluent allows introducing introduces an additional Equation (17):

$$\frac{\partial }{\partial \mathrm{t}}\left({\mathsf{\rho }\mathrm{Y}}_{\mathrm{i}}\right)+\nabla ·\left({\mathsf{\rho }\overline{\mathrm{v}}\mathrm{Y}}_{\mathrm{i}}\right)=-\nabla ·{\overline{\mathrm{J}}}_{\mathrm{i}}+{\mathrm{R}}_{\mathrm{i}}+{\mathrm{S}}_{\mathrm{i}}$$

(17)

where: ${\mathrm{R}}_{\mathrm{i}}$ is the index of net production of the substance i as a result of chemical reactions, and ${\mathrm{S}}_{\mathrm{i}}$ means the source term. Equation (17) is solved for N − 1 substances, where N is the total number of fluid phase substances present in the system.

The term ${\mathrm{R}}_{\mathrm{i}}$ expressing the formation rate of the substance i via chemical reactions for turbulent flows can be calculated using the Eddy dissipation model. This model of chemical–turbulent interaction assumes that the reaction rate is controlled by turbulence ignoring chemical time scales. Such an approach allows avoiding calculating the chemical kinetics of Arrhenius, which computationally are very demanding [41].

On the basis of this model, the production rate ${\mathrm{R}}_{\mathrm{i},\mathrm{r}}$ of a substance $\mathrm{i}$ as a result of the chemical reaction r. The ${\mathrm{R}}_{\mathrm{i},\mathrm{r}}$ is the lower value of the following two (18) and (19):

$${\mathrm{R}}_{\mathrm{i},\mathrm{r}}={\mathsf{\nu }}_{\mathrm{i},\mathrm{r}}^{\prime }{\mathrm{M}}_{\mathrm{w},\mathrm{i}}\mathrm{A}\mathsf{\rho }\frac{\mathsf{\epsilon }}{\mathrm{k}}{\min }_{ℛ}\left(\frac{{\mathrm{Y}}_{ℛ}}{{\mathsf{\nu }}_{ℛ,\mathrm{r}}^{\prime }{\mathrm{M}}_{\mathrm{w},ℛ}}\right)$$

(18)

$${\mathrm{R}}_{\mathrm{i},\mathrm{r}}={\mathsf{\nu }}_{\mathrm{i},\mathrm{r}}^{\prime }{\mathrm{M}}_{\mathrm{w},\mathrm{i}}\mathrm{AB}\mathsf{\rho }\frac{\mathsf{\epsilon }}{\mathrm{k}}\left(\frac{{{\displaystyle \sum }}_{\mathrm{P}}{\mathrm{Y}}_{\mathrm{P}}}{{{\displaystyle \sum }}_{\mathrm{j}}^{\mathrm{N}}{\mathsf{\nu }}_{\mathrm{j},\mathrm{r}}^{″}{\mathrm{M}}_{\mathrm{w},\mathrm{j}}}\right)$$

(19)

where: ${\mathrm{Y}}_{\mathrm{P}}$—the mass fraction of the product, $\mathrm{P}$, ${\mathrm{Y}}_{ℛ}$—the mass fraction of a considered reagent, $ℛ$, $\mathrm{A}$—an empirical constant of 4.0, and $\mathrm{B}$—an empirical constant of 0.5.

3.2. K-Omega BSL Turbulence Model

For an average turbulent fluid movement, each flow parameter at a given time can be expressed as the sum of the averaged value and its fluctuation. For a given velocity component, we can write:

$${\mathrm{u}}_{\mathrm{i}}=\overline{{\mathrm{u}}_{\mathrm{i}}}+{\mathrm{u}}_{\mathrm{i}}^{\prime }$$

(20)

In the same way, we may write scalars:

$$\mathsf{\varphi }=\overline{\mathsf{\varphi }}+{\mathsf{\varphi }}^{\prime }$$

(21)

where $\mathsf{\varphi }$ indicates a scalar, such as pressure, energy, or concentration.

The following form of transport equations in the tensor notation can be formulated:

$$\frac{\partial \mathsf{\rho }}{\partial \mathrm{t}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left({\mathsf{\rho }\mathrm{u}}_{\mathrm{i}}\right)=0$$

(22)

$$\frac{\partial }{\partial \mathrm{t}}\left({\mathsf{\rho }\mathrm{u}}_{\mathrm{i}}\right)+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left({\mathsf{\rho }\mathrm{u}}_{\mathrm{i}}{\mathrm{u}}_{\mathrm{j}}\right)=-\frac{\partial \mathrm{p}}{\partial {\mathrm{x}}_{\mathrm{i}}}+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left[\mathsf{\mu }\left(\frac{\partial {\mathrm{u}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{j}}}+\frac{\partial {\mathrm{u}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{i}}}-\frac{2}{3}{\mathsf{\delta }}_{\mathrm{ij}}\frac{\partial {\mathrm{u}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)\right]+\frac{\partial {\mathrm{R}}_{\mathrm{ij}}}{\partial {\mathrm{x}}_{\mathrm{j}}}$$

(23)

where R_ij is a tenor of Reynolds’ stress, related to the average rate of momentum change due to turbulent pulsations. In the Cartesian coordinate system, this tenor has a form consistent with Equation (24):

$${\mathrm{R}}_{\mathrm{ij}}=-\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{j}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{i}}^{\prime }}=\left[\begin{array}{ccc}-\mathsf{\rho }\overline{{\left({\mathrm{u}}_{\mathrm{x}}^{\prime }\right)}^2}& -\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{x}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{y}}^{\prime }}& -\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{x}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{z}}^{\prime }}\\ -\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{y}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{x}}^{\prime }}& -\mathsf{\rho }\overline{{\left({\mathrm{u}}_{\mathrm{y}}^{\prime }\right)}^2}& -\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{y}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{z}}^{\prime }}\\ -\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{z}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{x}}^{\prime }}& -\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{j}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{i}}^{\prime }}& -\mathsf{\rho }\overline{{\left({\mathrm{u}}_{\mathrm{z}}^{\prime }\right)}^2}\end{array}\right]$$

(24)

The problem of missing six equations for each component in the Reynolds tensor can be solved by solving the influence of turbulence, according to the Boussinesq approach [42]:

$${\mathrm{R}}_{\mathrm{ij}}=-\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{j}}^{\prime }}\overline{{\mathrm{u}}_{\mathrm{i}}^{\prime }}={\mathsf{\mu }}_{\mathrm{T}}\left(\frac{\partial {\overline{\mathrm{u}}}_{\mathrm{i}}}{\partial {\mathrm{x}}_{\mathrm{j}}}+\frac{\partial {\overline{\mathrm{u}}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{i}}}\right)-\frac{2}{3}\left(\mathsf{\rho }\mathrm{k}+{\mathsf{\mu }}_{\mathrm{T}}\frac{\partial {\mathrm{u}}_{\mathrm{k}}}{\partial {\mathrm{x}}_{\mathrm{k}}}\right){\mathsf{\delta }}_{\mathrm{ij}}$$

(25)

where ${\mathsf{\mu }}_{\mathrm{T}}$ stands for turbulent viscosity. It is a scalar function (often non-linear) of several variables, such as the physicochemical properties of the fluid or the nature of the flow.

The baseline turbulence (BSL) model k-ω was employed in the paper as it combines in a hybrid way the advantages of the k-ε model for free flow and the k-ω model solved in the wall area. The balance equations for the individual components of the model takes the following form:

$$\frac{\partial }{\partial \mathrm{t}}\left(\mathsf{\rho }\mathrm{k}\right)+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left({\mathsf{\rho }\mathrm{ku}}_{\mathrm{i}}\right)=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left({\mathrm{\Gamma }}_{\mathrm{k}}\frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)+{\mathrm{G}}_{\mathrm{k}}-{\mathrm{Y}}_{\mathrm{k}}+{\mathrm{S}}_{\mathrm{k}}$$

(26)

$$\frac{\partial }{\partial \mathrm{t}}\left(\mathsf{\rho }\mathsf{\omega }\right)+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left({\mathsf{\rho }\mathsf{\omega }\mathrm{u}}_{\mathrm{i}}\right)=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left({\mathrm{\Gamma }}_{\mathsf{\omega }}\frac{\partial \mathsf{\omega }}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)+{\mathrm{G}}_{\mathsf{\omega }}-{\mathrm{Y}}_{\mathsf{\omega }}+{\mathrm{D}}_{\mathsf{\omega }}+{\mathrm{S}}_{\mathsf{\omega }}$$

(27)

where ${\mathrm{G}}_{\mathrm{k}}$ expresses the kinetic energy turbulence production and is defined in the same way as in the standard model k-ω and ${\mathrm{G}}_{\mathsf{\omega }}$ represents the production of $\mathsf{\omega }$. Terms ${\mathrm{\Gamma }}_{\mathrm{k}}$ and ${\mathrm{\Gamma }}_{\mathsf{\omega }}$ are the effective diffusion coefficients of the k and ω variables, while ${\mathrm{Y}}_{\mathrm{k}}$ and ${\mathrm{Y}}_{\mathsf{\omega }}$ are the terms of the variable dissipation of the $\mathrm{k}$ and $\mathsf{\omega }$ due to turbulence, ${\mathrm{D}}_{\mathsf{\omega }}$ is the cross-diffusion term, while ${\mathrm{S}}_{\mathrm{k}}$ and ${\mathrm{S}}_{\mathsf{\omega }}$ stand for user-defined function sources.

The above mentioned effective diffusion coefficients are determined by the following Equations (28) and (29):

$${\mathrm{\Gamma }}_{\mathrm{k}}=\mathsf{\mu }+\frac{{\mathsf{\mu }}_{\mathrm{t}}}{{\mathsf{\sigma }}_{\mathrm{k}}}$$

(28)

$${\mathrm{\Gamma }}_{\mathsf{\omega }}=\mathsf{\mu }+\frac{{\mathsf{\mu }}_{\mathrm{t}}}{{\mathsf{\sigma }}_{\mathsf{\omega }}}$$

(29)

where ${\mathsf{\sigma }}_{\mathrm{k}}$ and ${\mathsf{\sigma }}_{\mathsf{\omega }}$ are mean turbulent Prandtl’s numbers for the $\mathrm{k}$ and $\mathsf{\omega }$ variables, defined by (30) and (31):

$${\mathsf{\sigma }}_{\mathrm{k}}=\frac{1}{\frac{{\mathrm{F}}_1}{{\mathsf{\sigma }}_{\mathrm{k},1}}+\frac{\left(1-{\mathrm{F}}_1\right)}{{\mathsf{\sigma }}_{\mathrm{k},2}}}$$

(30)

$${\mathsf{\sigma }}_{\mathsf{\omega }}=\frac{1}{\frac{{\mathrm{F}}_1}{{\mathsf{\sigma }}_{\mathsf{\omega },1}}+\frac{\left(1-{\mathrm{F}}_1\right)}{{\mathsf{\sigma }}_{\mathsf{\omega },2}}}$$

(31)

Turbulent viscosity ${\mathsf{\mu }}_{\mathrm{t}}$ is expressed by the Equation (32):

$${\mathsf{\mu }}_{\mathrm{t}}={\mathsf{\alpha }}^{*}\frac{\mathsf{\rho }\mathrm{k}}{\mathsf{\omega }}$$

(32)

where ${\mathsf{\alpha }}^{*}$ stands the correction factor for low Reynolds numbers given by Equation (33):

$${\mathsf{\alpha }}^{*}={\mathsf{\alpha }}_{\infty }^{*}\left(\frac{{\mathsf{\alpha }}_0^{*}+{\mathrm{Re}}_{\mathrm{t}}/{\mathrm{R}}_{\mathrm{k}}}{1+{\mathrm{Re}}_{\mathrm{t}}/{\mathrm{R}}_{\mathrm{k}}}\right)$$

(33)

and ${\mathrm{Re}}_{\mathrm{t}}=\frac{\mathsf{\rho }\mathrm{k}}{\mathsf{\mu }\mathsf{\omega }}$, ${\mathrm{R}}_{\mathrm{k}}=6$, ${\mathsf{\alpha }}_0^{*}=\frac{{\mathsf{\beta }}_{\mathrm{i}}}{3}$, and ${\mathsf{\beta }}_{\mathrm{i}}=0.072$.

The mixing function ${\mathrm{F}}_1$ is defined by the Equation (34):

$${\mathrm{F}}_1=\tan \mathrm{h}\left({\mathrm{\Phi }}_1^4\right)$$

(34)

where: ${\mathrm{\Phi }}_1=\min \left[\max \left(\frac{\sqrt{\mathrm{k}}}{0.09\mathsf{\omega }\mathrm{y}},\frac{500\mathsf{\mu }}{{\mathsf{\rho }\mathrm{y}}^2\mathsf{\omega }}\right),\frac{4\mathsf{\rho }\mathrm{k}}{{\mathsf{\sigma }}_{\mathsf{\omega },2}{\mathrm{D}}_{\mathsf{\omega }}^{+}{\mathrm{y}}^2}\right]$, ${\mathrm{D}}_{\mathsf{\omega }}^{+}=\max \left[2\mathsf{\rho }\frac{1}{{\mathsf{\sigma }}_{\mathsf{\omega },2}}\frac{1}{\mathsf{\omega }}\frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{j}}}\frac{\partial \mathsf{\omega }}{\partial {\mathrm{x}}_{\mathrm{j}}},{10}^{-10}\right],$ and $\mathrm{y}$ is the distance to the next surface.

The term ${\mathrm{G}}_{\mathrm{k}}$ in Equation (26) is defined by Equation (35):

$${\mathrm{G}}_{\mathrm{k}}=-\mathsf{\rho }\overline{{\mathrm{u}}_{\mathrm{i}}^{\prime }{\mathrm{u}}_{\mathrm{j}}^{\prime }}\frac{\partial {\mathrm{u}}_{\mathrm{j}}}{\partial {\mathrm{x}}_{\mathrm{i}}}$$

(35)

and the term ${\mathrm{G}}_{\mathsf{\omega }}$ in Equation (27) describes the Equation (36):

$${\mathrm{G}}_{\mathsf{\omega }}=\frac{{\mathsf{\alpha }\mathsf{\alpha }}^{*}}{{\mathsf{\nu }}_{\mathrm{t}}}{\mathrm{G}}_{\mathrm{k}}$$

(36)

The dissipation of kinetic energy of turbulence ${\mathrm{Y}}_{\mathrm{k}}$ is defined similarly to the standard model $\mathrm{k}-\mathsf{\omega }$ and determined by Equation (37):

$${\mathrm{Y}}_{\mathrm{k}}={\mathsf{\rho }\mathsf{\beta }}^{*}\mathrm{k}\mathsf{\omega }$$

(37)

For the k-omega BSL model the term ${\mathrm{f}}_{{\mathsf{\beta }}^{*}}$, is equal to 1.

The dissipation of the $\mathsf{\omega }$ variable is defined in a similar way as in the standard model. The difference constitutes in the way how ${\mathrm{f}}_{\mathsf{\beta }}$ and ${\mathsf{\beta }}_{\mathrm{i}}$ are calculated. In the model k-omega BSL, the value of ${\mathrm{f}}_{\mathsf{\beta }}$ is constant and equals 1, so ${\mathrm{Y}}_{\mathsf{\omega }}$ is ultimately expressed by Equation (38):

$${\mathrm{Y}}_{\mathrm{k}}={\mathsf{\rho }\mathsf{\beta }\mathsf{\omega }}^2$$

(38)

where ${\mathsf{\beta }}_{\mathrm{i}}$ is defined by as follows:

$${\mathsf{\beta }}_{\mathrm{i}}={\mathrm{F}}_1{\mathsf{\beta }}_{\mathrm{i},1}+\left(1-{\mathrm{F}}_1\right){\mathsf{\beta }}_{\mathrm{i},2}$$

(39)

The term of the cross-diffusion ${\mathrm{D}}_{\mathsf{\omega }}$ was formed by combining the equations of the model $\mathrm{k}-\mathsf{\epsilon }$ and $\mathrm{k}-\mathsf{\omega }$. Finally, this quantity is expressed by Equation (40):

$${\mathrm{D}}_{\mathsf{\omega }}=2\left(1-{\mathrm{F}}_1\right)\mathsf{\rho }\frac{1}{{\mathsf{\omega }\mathsf{\sigma }}_{\mathsf{\omega },2}}\frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{j}}}\frac{\partial \mathsf{\omega }}{\partial {\mathrm{x}}_{\mathrm{j}}}$$

(40)

There are many constants in the model equations, which are listed below: ${\mathsf{\sigma }}_{\mathrm{k},1}=2.0$, ${\mathsf{\sigma }}_{\mathsf{\omega },1}=2.0,$ ${\mathsf{\sigma }}_{\mathrm{k},2}=1.0,$ ${\mathsf{\sigma }}_{\mathsf{\omega },2}=1.168,$ ${\mathsf{\beta }}_{\mathrm{i},1}=0.075,$ and ${\mathsf{\beta }}_{\mathrm{i},2}=0.0828$. Additional constants used in the equations of the k-omega BSL model take the same values as for the standard model $\mathrm{k}-\mathsf{\omega }$.

3.3. Discrete Ordinates Radiation Model

A radiation model solving the radiation transport Equation (41) for a finite number of discrete solid angles was introduced into the CFD fluidized bed model. Radiative transfer equation (RTE) for the absorption, emitting, and dispersing medium in the position $\overline{\mathrm{r}}$ and in the direction $\overline{\mathrm{s}}$ has the following form:

$$\frac{\mathrm{dI}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)}{\mathrm{ds}}+\left(\mathrm{a}+{\mathsf{\sigma }}_{\mathrm{s}}\right)\mathrm{I}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)={\mathrm{an}}^2\frac{{\mathsf{\sigma }\mathrm{T}}^4}{\prod }+\frac{{\mathsf{\sigma }}_{\mathrm{s}}}{4\prod }\underset{0}{\overset{4\prod }{{\displaystyle \int }}}\mathrm{I}\left(\overline{\mathrm{r}},{\overline{\mathrm{s}}}^{\prime }\right)\mathrm{\Phi }\left(\overline{\mathrm{s}}·{\overline{\mathrm{s}}}^{\prime }\right)\mathrm{d}{\mathrm{\Omega }}^{\prime }$$

(41)

where: $\overline{\mathrm{r}}$—location vector, $\overline{\mathrm{s}}$—direction vector, ${\overline{\mathrm{s}}}^{\prime }$—vector of scattering direction, $\mathrm{s}$—path length, $\mathrm{a}—$absorption coefficient, $\mathrm{n}$—refraction index, ${\mathsf{\sigma }}_{\mathrm{s}}$—refraction index, $\mathsf{\sigma }$—Stefan-Boltzmann’s constant, $\mathrm{I}—$radiation intensity depending on location ($\overline{\mathrm{r}}$) and direction ($\overline{\mathrm{s}}$), $\mathrm{T}$—local temperature, and ${\mathrm{\Omega }}^{\prime }$—solid angle. Equation (41) is transformed into a transport equation for radiation intensity in the spatial coordinate system (x,y,z).

The discrete ordinates (DO) model uses the radiation transport equation in the direction $\overline{\mathrm{s}}$ as a surface equation:

$$\nabla \left(\mathrm{I}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)\overline{\mathrm{s}}\right)+\left(\mathrm{a}+{\mathsf{\sigma }}_{\mathrm{s}}\right)\mathrm{I}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)={\mathrm{an}}^2\frac{{\mathsf{\sigma }\mathrm{T}}^4}{\prod }+\frac{{\mathsf{\sigma }}_{\mathrm{s}}}{4\prod }\underset{0}{\overset{4\prod }{{\displaystyle \int }}}\mathrm{I}\left(\overline{\mathrm{r}},{\overline{\mathrm{s}}}^{\prime }\right)\mathrm{\Phi }\left(\overline{\mathrm{s}}·{\mathrm{s}¯}^{\prime }\right)\mathrm{d}{\mathrm{\Omega }}^{\prime }$$

(42)

It is also possible to consider non-grey radiation using a grey-band model. In this case, Equation (42) takes the Equation (43):

$$\nabla \left({\mathrm{I}}_{\mathsf{\lambda }}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)\overline{\mathrm{s}}\right)+\left({\mathrm{a}}_{\mathsf{\lambda }}+{\mathsf{\sigma }}_{\mathrm{s}}\right){\mathrm{I}}_{\mathsf{\lambda }}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)={\mathrm{a}}_{\mathsf{\lambda }}{\mathrm{I}}_{\mathrm{b}\mathsf{\lambda }}+\frac{{\mathsf{\sigma }}_{\mathrm{s}}}{4\prod }\underset{0}{\overset{4\prod }{{\displaystyle \int }}}{\mathrm{I}}_{\mathsf{\lambda }}\left(\overline{\mathrm{r}},{\overline{\mathrm{s}}}^{\prime }\right)\mathrm{\Phi }\left(\overline{\mathrm{s}}·{\overline{\mathrm{s}}}^{\prime }\right)\mathrm{d}{\mathrm{\Omega }}^{\prime }$$

(43)

In the above equation, λ is the wavelength, ${\mathrm{a}}_{\mathsf{\lambda }}$ is the absorption coefficient and ${\mathrm{I}}_{\mathrm{b}\mathsf{\lambda }}$ states for the radiation intensity of the black body, determined by Planck’s function.

Such an implementation of the DO model divides the radiation spectrum into N bands with different wavelength ranges. These bands do not have to be adjacent or equal. The user enters the wavelength ranges in the individual bands himself. The RTE equation is then integrated into each defined band, leading to transport equations of ${\mathrm{I}}_{\mathsf{\lambda }}$ Δλ, the radiation energy in each wavelength range Δλ. Then the radiation in each band is considered to be grey.

The emissivity of a black body in a given range of wavelengths per unit solid angle is defined as (44):

$$\left[\mathrm{F}\left(0\to {\mathrm{n}\mathsf{\lambda }}_2\mathrm{T}\right)-\mathrm{F}\left(0\to {\mathrm{n}\mathsf{\lambda }}_1\mathrm{T}\right)\right]{\mathrm{n}}^2\frac{{\mathsf{\sigma }\mathrm{T}}^4}{\prod }$$

(44)

where $\mathrm{F}\left(0\to \mathrm{n}\mathsf{\lambda }\mathrm{T}\right)$ stands for a portion of the radiation energy emitted by the black body.

The total radiation intensity $\mathrm{I}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)$ in each direction $\overline{\mathrm{s}}$ in the position $\overline{\mathrm{r}}$ is calculated by summing up all wavelength ranges according to the Equation (45):

$$\mathrm{I}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right)={\displaystyle \sum }_{\mathrm{k}}{\mathrm{I}}_{{\mathsf{\lambda }}_{\mathrm{k}}}\left(\overline{\mathrm{r}},\overline{\mathrm{s}}\right){\mathsf{\Delta }\mathsf{\lambda }}_{\mathrm{k}}$$

(45)

The radiation DO model in the considered case is complemented by the weighted sum of gray gases model (WSGGM) to calculate the variable value of the absorption coefficient.

The model makes a compromise between the excessive simplification of the gray gas model and the complete model considering absorption in individual bands.

The WSGGM model, based on the assumption that total emissivity on the $\mathrm{s}$ path, can be expressed as:

$$\mathsf{\epsilon }={\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{I}}{\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}\left(\mathrm{T}\right)(1-{\mathrm{e}}^{-{\mathsf{\kappa }}_{\mathrm{i}}\mathrm{ps}})$$

(46)

where ${\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}$ stands for the emissivity factor for the fictitious grey gas i, the value in brackets is the emissivity of the fictitious grey gas i and ${\mathsf{\kappa }}_{\mathrm{i}}$ is the absorption coefficient of the grey gas i. The value $\mathrm{p}$ is the sum of the partial pressures of all the absorbing component gases.

Terms ${\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}$ and ${\mathsf{\kappa }}_{\mathrm{i}}$ are obtained according to [43,44] and depend on the gas composition, but the coefficient ${\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}$ is also a function of temperature. If the total pressure is not equal to 1 atm, the scaling principle for the coefficient ${\mathsf{\kappa }}_{\mathrm{i}}$ is applied.

For ${\mathrm{p}}_{\mathrm{tot}}<0.9\mathrm{atm}$ or ${\mathrm{p}}_{\mathrm{tot}}>1.1\mathrm{atm}$ a correction is introduced, according to the expression:

$${\mathsf{\kappa }}_{\mathrm{i}}\to {\mathsf{\kappa }}_{\mathrm{i}}{\mathrm{p}}_{\mathrm{tot}}{}^{\mathrm{m}}$$

(47)

where: $\mathrm{m}$ is determined according to [45] and depends on partial pressures and temperatures of the absorbing gas components.

The absorption coefficient for $\mathrm{i}=0$ is equal to 0 to include windows in the radiation spectrum between highly absorbent regions $({{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{I}}{\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}<1)$, whereas the weight for $\mathrm{i}=0$ is determined from the Equation (48):

$${\mathrm{a}}_{\mathsf{\epsilon },0}=1-{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{I}}{\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}$$

(48)

The influence of temperature on the coefficient ${\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}$ is usually approximated using the following formula:

$${\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}={\displaystyle \sum }_{\mathrm{j}=1}^{\mathrm{J}}{\mathrm{b}}_{\mathsf{\epsilon },\mathrm{i},\mathrm{j}}{\mathrm{T}}^{\mathrm{j}-1}$$

(49)

where ${\mathrm{b}}_{\mathsf{\epsilon },\mathrm{i},\mathrm{j}}$ stands for polynomial coefficients of temperature emissivity of gases, while ${\mathrm{b}}_{\mathsf{\epsilon },\mathrm{i},\mathrm{j}}$ and ${\mathsf{\kappa }}_{\mathrm{i}}$ are estimated by matching Equation (46) to the total emissivity obtained experimentally [43,44,46].

The coefficients ${\mathrm{b}}_{\mathsf{\epsilon },\mathrm{i},\mathrm{j}}$ and ${\mathsf{\kappa }}_{\mathrm{i}}$ can be considered as constant, as they change slightly with $\mathrm{ps}$ and $\mathrm{T},$ according to [43,44]. If ${\mathsf{\kappa }}_{\mathrm{i}}\mathrm{ps}“1$ for all values $\mathrm{i}$ Equation (46) can be simplified as follows (50):

$$\mathsf{\epsilon }={\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{I}}{\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}{\mathsf{\kappa }}_{\mathrm{i}}\mathrm{p}$$

(50)

When comparing Equation (50) with the grey gas model with absorption coefficient $\mathsf{\alpha }$, it is observed that the radiation intensity change over distance s in the WSGGM model is the same as in the grey gas model with $\mathsf{\alpha }$:

$$\mathsf{\alpha }={\displaystyle \sum }_{\mathrm{i}=0}^{\mathrm{I}}{\mathrm{a}}_{\mathsf{\epsilon },\mathrm{i}}{\mathsf{\kappa }}_{\mathrm{i}}\mathrm{p}$$

(51)

In the general case $\mathsf{\alpha }$ is calculated as (52):

$$\mathsf{\alpha }=-\frac{\ln \left(1-\mathsf{\epsilon }\right)}{\mathrm{s}}$$

(52)

where the emissivity ε for the WSGGM model is calculated from the Equation (36).

The factor $\mathsf{\alpha }$ defined by Equation (52) is dependent on $\mathrm{s}$, which reflects the non-grey nature of radiation absorption in molecular gases. ANSYS Fluent employs the Equation (51) for $\mathrm{s}\le {10}^{-4}\mathrm{m}$ and the Equation (52) for $\mathrm{s}>{10}^{-4}\mathrm{m}$. For the border value $\mathrm{s}$, both equations lead to an almost identical value.

3.4. Discrete Phase Modeling (DPM)

Currently, two approaches to numerical modeling of multiphase flows are most commonly used: Euler–Lagrange or Euler–Euler method. In the fluidized bed model, the Euler–Lagrange method was adapted. In this method, the fluid phase is treated as continuous by solving Navier–Stokes equations, while the dispersed phase is calculated by tracking a large number of particles, bubbles or droplets in the computing domain. The dispersed phase can exchange momentum, mass, and energy with a continuous fluid phase.

This approach becomes much simpler when interactions between particles can be neglected. This assumption requires the dispersed phase to have a small volume share. However, a high flow load of the dispersed phase (${\dot{\mathrm{m}}}_{\mathrm{particles}}\ge {\dot{\mathrm{m}}}_{\mathrm{fluid}}$) is also acceptable. Particle or droplet trajectories are calculated individually at specified intervals when performing fluid phase calculations making the approach suitable for modeling spray dryers, combustion of coal, and liquid fuels as well as two-phase flows.

The force equation can be expressed by (53):

$$\frac{{\mathrm{d}\overline{\mathrm{u}}}_{\mathrm{p}}}{\mathrm{dt}}=\frac{\overline{\mathrm{u}}-{\overline{\mathrm{u}}}_{\mathrm{p}}}{{\mathsf{\tau }}_{\mathrm{r}}}+\frac{\overline{\mathrm{g}}\left({\mathsf{\rho }}_{\mathrm{p}}-\mathsf{\rho }\right)}{{\mathsf{\rho }}_{\mathrm{p}}}+\overline{\mathrm{F}}$$

(53)

where $\overline{\mathrm{F}}$ means a term of additional acceleration (force/unit of mass of the particle), the term $\frac{\overline{\mathrm{u}}-{\overline{\mathrm{u}}}_{\mathrm{p}}}{{\mathsf{\tau }}_{\mathrm{r}}}$ stands for the drag force corresponding to the unit of mass of the particle.

The particle relaxation time ${\mathsf{\tau }}_{\mathrm{r}}$ is defined as follows:

$${\mathsf{\tau }}_{\mathrm{r}}=\frac{{\mathsf{\rho }}_{\mathrm{p}}{\mathrm{d}}_{\mathrm{p}}^2}{18\mathsf{\mu }}\frac{14}{{\mathrm{C}}_{\mathrm{d}}\mathrm{Re}}$$

(54)

where $\overline{\mathrm{u}}$ and ${\overline{\mathrm{u}}}_{\mathrm{p}}$ mean fluid and particles velocities, respectively, $\mathsf{\mu }$ expresses the molecular viscosity, $\mathsf{\rho }$ and ${\mathsf{\rho }}_{\mathrm{p}}$ are fluid and particle density, respectively and ${\mathrm{d}}_{\mathrm{p}}$ stands for the particle diameter. The Re expresses Reynolds number:

$$\mathrm{Re}\equiv \frac{{\mathsf{\rho }\mathrm{d}}_{\mathrm{p}}\left|{\overline{\mathrm{u}}}_{\mathrm{p}}-\overline{\mathrm{u}}\right|}{\mathsf{\mu }}$$

(55)

Particle rotation has a significant impact on the trajectory of a given particle, moving in a fluid. An additional differential equation for the particle momentum (56) was employed to take into account this effect of particle rotation:

$${\mathrm{I}}_{\mathrm{p}}\frac{{\mathrm{d}\overline{\mathsf{\omega }}}_{\mathrm{p}}}{\mathrm{dt}}=\frac{{\mathsf{\rho }}_{\mathrm{f}}}{2}{\left(\frac{{\mathrm{d}}_{\mathrm{p}}}{2}\right)}^5{\mathrm{C}}_{\mathrm{w}}\overline{\mathrm{\Omega }}=\overline{\mathrm{T}}$$

(56)

where ${\mathrm{I}}_{\mathrm{p}}$ is the moment of inertia, ${\overline{\mathsf{\omega }}}_{\mathrm{p}}$ expresses the angular velocity of the particle, ${\mathsf{\rho }}_{\mathrm{f}}$ is the density of the fluid, ${\mathrm{d}}_{\mathrm{p}}$ is the diameter of the particle, ${\mathrm{C}}_{\mathrm{w}}$ is the rotational resistance coefficient, $\overline{\mathrm{T}}$ is the torque applied to the particle, and $\overline{\mathrm{\Omega }}$ means the relative rotational speed between the particle and the fluid. The value $\overline{\mathrm{\Omega }}$ is determined as (57):

$$\overline{\mathrm{\Omega }}=\frac{1}{2}\nabla {\mathrm{x}\overline{\mathrm{u}}}_{\mathrm{f}}-{\overline{\mathsf{\omega }}}_{\mathrm{p}}$$

(57)

The moment of inertia for a spherical particle is calculated from Equation (58):

$${\mathrm{I}}_{\mathrm{p}}=\frac{\prod }{60}{\mathsf{\rho }}_{\mathrm{f}}{\mathrm{d}}_{\mathrm{p}}^5$$

(58)

Equation (53) contains the term $\overline{\mathrm{F}}$, responsible for additional forces acting on the particle. They may be particularly relevant in specific circumstances. The first is the “virtual mass” force, which is the force needed to accelerate the fluid surrounding the particle. This force is determined by Equation (59):

$$\overline{\mathrm{F}}={\mathrm{C}}_{\mathrm{vm}}\frac{\mathsf{\rho }}{{\mathsf{\rho }}_{\mathrm{p}}}\left({\overline{\mathrm{u}}}_{\mathrm{p}}\nabla \overline{\mathrm{u}}-\frac{{\mathrm{d}\overline{\mathrm{u}}}_{\mathrm{p}}}{\mathrm{dt}}\right)$$

(59)

where: ${\mathrm{C}}_{\mathrm{vm}}$ stands for a virtual mass factor with a default value of 0.5.

This additional force is created by the pressure gradient (60):

$$\overline{\mathrm{F}}=\frac{\mathsf{\rho }}{{\mathsf{\rho }}_{\mathrm{p}}}{\overline{\mathrm{u}}}_{\mathrm{p}}\nabla \overline{\mathrm{u}}$$

(60)

Virtual mass and forces related to the pressure gradient are not relevant when the density of a fluid is significantly lower than the density of particles, as in the considered here case with solid particles in gaseous flows $\left(\frac{\mathsf{\rho }}{{\mathsf{\rho }}_{\mathrm{p}}}\text{”}1\right)$. However, it is recommended to take them into account for density ratios higher than 0.1.

When particles are suspended in gas with large temperature gradients, an additional force acts on the particles in the opposite direction to this gradient. This phenomenon is referred to as thermophoresis (or otherwise thermodiffusion, Soret effect) and can be considered by an additional acceleration term in the Formula (53). It is expressed by Equation (61):

$$\overline{\mathrm{F}}=-{\mathrm{D}}_{\mathrm{T},\mathrm{p}}\frac{1}{{\mathrm{m}}_{\mathrm{p}}\mathrm{T}}\nabla \mathrm{T}$$

(61)

where ${\mathrm{D}}_{\mathrm{T},\mathrm{p}}$ expresses the thermophoresis thermoforming factor. This coefficient can be expressed as a constant, a polynomial, user-defined function, or can be expressed in the form proposed by Talbot [47], making that the Equation (61) takes the following form:

$$\overline{\mathrm{F}}=-\frac{6\prod {\mathrm{d}}_{\mathrm{p}}{\mathsf{\mu }}^2{\mathrm{C}}_{\mathrm{s}}\left(\mathrm{K}+{\mathrm{C}}_{\mathrm{t}}\mathrm{Kn}\right)}{\mathsf{\rho }\left(1+3{\mathrm{C}}_{\mathrm{m}}\mathrm{Kn}\right)\left(1+2\mathrm{K}+2{\mathrm{C}}_{\mathrm{t}}\mathrm{Kn}\right)}\frac{1}{{\mathrm{m}}_{\mathrm{p}}\mathrm{T}}\nabla \mathrm{T}$$

(62)

where $\mathrm{Kn}$—Knudsen number $\left(\frac{2\mathsf{\lambda }}{{\mathrm{d}}_{\mathrm{p}}}\right)$, $\mathsf{\lambda }$—mean free fluid path, $\mathrm{K}=\mathrm{k}/{\mathrm{k}}_{\mathrm{p}}$, $\mathrm{k}$—thermal conductivity of the fluid based on translation energy $\left(\frac{15}{4}\mathsf{\mu }\mathrm{R}\right)$, ${\mathrm{k}}_{\mathrm{p}}$—thermal conductivity of the particle, ${\mathrm{C}}_{\mathrm{s}}=1.17$, ${\mathrm{C}}_{\mathrm{t}}=2.18$, ${\mathrm{C}}_{\mathrm{m}}=1.14$, ${\mathrm{m}}_{\mathrm{p}}$—mass of the particle, $\mathrm{T}$—local fluid temperature, and $\mathsf{\mu }$—fluid viscosity.

3.5. Finite Volume Method

For this work, the finite volume method was used. This method consists of dividing the examined area into a finite number of control cells with computational nodes. Then the differential equations describing a phenomenon are integrated into the control volume of each cell. Different types of approximation expressions replace individual parts of the equations. Finally, a system of algebraic equations, usually non-linear, is obtained [48].

The general form of the equation, which describes the transport of a variable in the differential form, is as follows [42]:

$$\frac{\partial \left(\mathsf{\rho }\mathsf{\varphi }\right)}{\partial \mathrm{t}}+\mathrm{div}\left(\mathsf{\rho }\mathsf{\varphi }\overline{\mathrm{u}}\right)=\mathrm{div}\left({\mathrm{\Gamma }}_{\mathsf{\varphi }}\mathrm{grad}\mathsf{\varphi }\right)+{\mathrm{S}}_{\mathsf{\varphi }}$$

(63)

The variable $\mathsf{\varphi }$ expresses a variable in the equations used to describe fluid flow. Replacing $\mathsf{\varphi }$ with an appropriate variable (u, v, and w) and selecting the appropriate diffusion coefficient, the individual transport equations described in Chapter 3.1 may be obtained.

The control volume method uses the general transport Formula (63). This equation is integrated over the entire volume of the control cell, and then transformed using the Ostrogradsky–Gauss theorem [42]:

$$\underset{\mathrm{CV}}{\overset{}{{\displaystyle \int }}}\mathrm{diva}\mathrm{dV}=\underset{\mathrm{A}}{\overset{}{{\displaystyle \int }}}\mathrm{n}·\mathrm{adA}$$

(64)

Finally, the following equation is obtained [42]:

$$\frac{\partial }{\partial \mathrm{t}}\left(\underset{\mathrm{CV}}{\overset{}{{\displaystyle \int }}}\left(\mathsf{\rho }\mathsf{\varphi }\right)\mathrm{dV}\right)+\underset{\mathrm{A}}{\overset{}{{\displaystyle \int }}}\mathrm{n}·\left(\mathsf{\rho }\mathsf{\varphi }\mathrm{u}\right)\mathrm{dA}=\underset{\mathrm{A}}{\overset{}{{\displaystyle \int }}}\mathrm{n}·\left(\mathrm{\Gamma }\mathrm{grad}\mathsf{\varphi }\right)\mathrm{dA}+\underset{\mathrm{CV}}{\overset{}{{\displaystyle \int }}}{\mathrm{S}}_{\mathsf{\varphi }}\mathrm{dV}$$

(65)

The first term of the above formula is the non-stationary term, the second one represents the convective transport of the dependent variable, and the two last ones are the diffusion term and the source term, respectively.

For stationary problems, the first term disappears from Equation (65). However, when considering non-stationary conditions, it is necessary to take into account the change in the value of a dependent variable over time as well. For this purpose, Formula (65) must be integrated over time. The obtained Equation (66) represents the most general form of the transport equation:

$$\underset{\Delta \mathrm{t}}{\overset{}{{\displaystyle \int }}}\frac{\partial }{\partial \mathrm{t}}\left(\underset{\mathrm{CV}}{\overset{}{{\displaystyle \int }}}\left(\mathsf{\rho }\mathsf{\varphi }\right)\mathrm{dV}\right)\mathrm{dt}+\underset{\Delta \mathrm{t}}{\overset{}{{\displaystyle \int }}}\underset{\mathrm{A}}{\overset{}{{\displaystyle \int }}}\mathrm{n}\left(\mathsf{\rho }\mathsf{\varphi }\mathrm{u}\right)\mathrm{dAdt}=\underset{\Delta \mathrm{t}}{\overset{}{{\displaystyle \int }}}\underset{\mathrm{A}}{\overset{}{{\displaystyle \int }}}\mathrm{n}\left(\mathrm{\Gamma }\mathrm{grad}\mathsf{\varphi }\right)\mathrm{dAdt}+\underset{\Delta \mathrm{t}}{\overset{}{{\displaystyle \int }}}\underset{\mathrm{CV}}{\overset{}{{\displaystyle \int }}}{\mathrm{S}}_{\mathsf{\varphi }}\mathrm{dVdt}$$

(66)

This transformed equation is used in the control volume method for each control cell in the considered computational domain. However, to solve it, the equation must be discretized, i.e., presented in a form that can be solved in each cell into which the area was divided. The equation is then obtained in a discretized form [49]:

$$\frac{\partial \mathsf{\rho }\mathsf{\varphi }}{\partial \mathrm{t}}\mathrm{V}+{\displaystyle \sum }_{\mathrm{f}}^{\mathrm{N}\mathrm{faces}}{\mathsf{\rho }}_{\mathrm{f}}{\overline{\mathrm{u}}}_{\mathrm{f}}{\mathsf{\varphi }}_{\mathrm{f}}·{\overline{\mathrm{A}}}_{\mathrm{f}}={\displaystyle \sum }_{\mathrm{f}}^{\mathrm{N}\mathrm{faces}}{\mathrm{\Gamma }}_{\mathsf{\varphi }}\nabla {\mathsf{\varphi }}_{\mathrm{f}}·{\overline{\mathrm{A}}}_{\mathrm{f}}+{\mathrm{S}}_{\mathsf{\varphi }}\mathrm{V}$$

(67)

where the $\mathrm{N}\mathrm{faces}$ term means the number of walls limiting the control cell, ${\mathsf{\varphi }}_{\mathrm{f}}$ represents the value of the $\mathsf{\varphi }$ variable penetrating through the cell wall ($\mathrm{f}$), ${\mathsf{\rho }}_{\mathrm{f}}{\overline{\mathrm{u}}}_{\mathrm{f}}{\overline{\mathrm{A}}}_{\mathrm{f}}$ is the mass flow through the cell wall, $\nabla {\mathsf{\varphi }}_{\mathrm{f}}$ stands for the gradient of the $\mathsf{\varphi }$ variable on the $\mathrm{f}$ wall, and $\mathrm{V}$ is the cell volume.

In practice, most of the equations describing real phenomena are non-linear, so Equation (67) can still have a non-linear character, expressed, e.g., by the non-linear dependencies. Therefore, Equation (67) is transformed using various approximation methods. The source term is also linearized. As a result of all these operations, a system of linear algebraic equations can finally be obtained:

$${\mathrm{a}}_{\mathrm{P}}{\mathsf{\varphi }}_{\mathrm{P}}={\displaystyle \sum }_{\mathrm{nb}}{\mathrm{a}}_{\mathrm{nb}}{\mathsf{\varphi }}_{\mathrm{nb}}+\mathrm{b}$$

(68)

The index “nb” means adjacent control cells relative to a cell with a center P, and the variable b refers to the source term, ${\mathrm{a}}_{\mathrm{P}}$ and ${\mathrm{a}}_{\mathrm{nb}}$ are linearized coefficients for ${\mathsf{\varphi }}_{\mathrm{P}}$ and ${\mathsf{\varphi }}_{\mathrm{nb}}$. The resulting system of linear equations is solved using the iterative method.

A pressure-based solver was used in fluidized bed modeling. In this approach, the velocity field is obtained by solving the pressure equation (pressure correction) from the continuity and momentum equations. This means that if the correct pressure field is used to calculate the velocity field, then the calculated velocity field will satisfy the continuity equation. This requires the use of an iterative process, wherein each successive iteration values of pressures (and pressure corrections) and velocities are calculated until they meet the continuity equation with the required accuracy. According to the Equation (67), it is necessary to know the value of the dependent variable ϕ on the walls of a control cell. It must be obtained on the basis of knowledge of the variable value in computational nodes, located in the centers of control cells. Interpolation procedures are used for this purpose.

The second order upwind interpolation scheme was used in the study. In this method, the value of the dependent variable ${\mathsf{\varphi }}_{\mathrm{f}}$ on the cell’s wall is calculated based on the values in the center of the nearest cell lying in the opposite direction to the fluid movement, according to the following expression [49]:

$${\mathsf{\varphi }}_{\mathrm{f}}={\mathsf{\varphi }}_{\mathrm{upstream}}+\nabla {\mathsf{\varphi }}_{\mathrm{upstream}}·\overline{\mathrm{r}}$$

(69)

where the index ‘upstream’ means the closest neighboring cell situated in the opposite direction to the direction of the fluid flow, and $\overline{\mathrm{r}}$ stands for the vector connecting the center of that cell with the center point on the wall for which the value of the variable $\mathsf{\varphi }$ is calculated.

In the case of a regular grid, it can be written as follows:

$${\mathsf{\varphi }}_{\mathrm{e}}={\mathsf{\varphi }}_{\mathrm{P}}+\nabla {\mathsf{\varphi }}_{\mathrm{P}}·\frac{1}{2}\Delta \mathrm{x}$$

(70)

which will lead to the following expression:

$${\mathsf{\varphi }}_{\mathrm{e}}\cong \frac{3}{2}{\mathsf{\varphi }}_{\mathrm{P}}-\frac{1}{2}{\mathsf{\varphi }}_{\mathrm{W}}$$

(71)

3.6. Interphase Interactions, Particle–Particle Interaction, and Bed Material Recirculation

It has to be noted that in the considered case of the reactor, enough particles are present such that momentum exchange between dispersed and carrier phase interfaces alters the dynamics of the carrier phase. Due to this fact, it was obviously required to pay attention to the interphase effects that are especially important from the point of view of the turbulence modeling.

The source term F for the momentum exchange due to particle flow through the control volume has been calculated based on the general rule expressed by Equation (72) below [49]:

$$F={{\displaystyle \sum }}^{}\left(\frac{18\mu C_{D}Re}{{\rho }_p{d}_p^{2}24}\left(u_p-u\right)+F_{pther}\right){\dot{m}}_p\Delta t$$

(72)

where: $\mu —$viscosity of the fluid phase material, ${\rho }_p—$the density of the particle, $d_p—$the diameter of the particle, $Re—$Reynolds number, $u_p—$the velocity of the solid phase particle, $u—$velocity of the fluid element, $C_D—$drag coefficient, ${\dot{m}}_p—$the mass flow rate of the particles, $\Delta t—$time step, and $F_{other}—$other interaction forces.

The numerical model that has been used to carry out the computation process allowed us to simulate turbulent multiphase flow based on the principles formulated by Faeth [50] and Amsden [51]. The additional turbulence production in the appropriate continuous phase transport equation in this approach was considered when the particle diameter was higher than 10% of the turbulence length scale.

The heat transfer Q (source or sink) between the continuous phase and solid particles was incorporated into the model by applying the following Equation [49]:

$$Q=\left[\frac{{\overline{m}}_p}{m_{p,0}}{C}_p\Delta T_p+\frac{\Delta m_p}{m_{p,0}}\left(-h_{fg}+h_{pyr}+{{\displaystyle \int }}_{T_{ref}}^{T_p}{C}_{p,i}dT\right)\right]{\dot{m}}_{p,0}+Q_{SR}$$

(73)

where: ${\overline{m}}_p—$the average mass of the particle in the control volume, $m_{p,0}—$the initial mass of the particle, $C_p—$heat capacity of the particle, $\Delta T_p—$the temperature change of the particle in the control volume, $\Delta m_p—$change in the mass of the particle in the control volume, $h_{fg}—$latent heat of volatiles, $h_{pyr}—$the heat of pyrolysis, $C_{p,i}—$heat capacity of the volatiles, $T_p—$the temperature of the particle upon exit of the control volume, $T_{ref}—$reference temperature for enthalpy, ${\dot{m}}_{p,0}—$the initial mass flow rate of the particle injection tracked, and $Q_{SR}—$source expressing the thermal effect of surface reaction (combustion).

A modified stochastic collisions approach has been applied to describe particle–particle interactions collisions. The modification assumed the introduction of so-called “parcels” that represent in statistical meaning several individual solid matter particles. Thus one parcel includes many particles that efficiently reduce computation cost. Collisions of the grouped particles have been estimated stochastically based on the O’Rourke algorithm [52]. This approach assumes two kinds of collision effects: coalescence or/and reflection. The probability of each above-mentioned phenomena is calculated by the collision Weber number (We_c) [52]:

$$W{e}_c=\frac{\rho U_{rel}^2\overline{D}}{\sigma }$$

(74)

where: $U_{rel}$—the relative velocity between two parcels, $\overline{D}$—the average diameter of two parcels, and $\mathsf{\sigma }$—particle Surface stress.

The above-given expression is calculated for each couple of parcels for which collision has been found.

The approach to control solid-phase mass conservation consisted of controlling the total instantaneous mass of inert material in the computational domain.

The experimental data obtained from the real full-scale boiler was applied in the developed CFD model. In order to take into consideration bed material recirculation, the mass flow rate of the solids coming back to the bed was set according to the average mass flow rate declared in the technical specification of the boiler. Such an approach ensured the mass balance of the solid phase in the boiler based on the monitoring of the solid phase total mass in the computational domain. Although this simplification of the system may affect the agreement between conditions in the simulation and real case, such an approach is acceptable due to the goal of the study, concerning the evaluation of the parameter difference for different operation variants. Assessment of the relative change of selected flow variables was still possible and provided the answer regarding final recommendations.

The final confirmation of the correctness of the developed model will be confirmed in the validation chapter.

4. Geometry and Computational Domain Discretization of the Combustion Chamber

The numerical model based on computational fluid mechanics was performed using the ANSYS Workbench 17.2 platform. The geometric model of the boiler was developed in the ANSYS SpaceClaim Direct Modeler. The real dimensions of the tested part of the installation were applied to build spatial geometry. Due to the presence of the geometric and physicochemical symmetry plane, it was possible to reduce the geometry to half of the considered object. Such an approach allowed studying the system with a reduced number of grid control cells, which had a positive effect on the computation time (Figure 2).

The geometry was divided into orthogonal components (Figure 3, Figure 4 and Figure 5), with the whole domain defined as topologically coherent, which also allowed to generate a coherent computing grid, not requiring the use of so-called computational coupling interfaces of individual parts forming the spatial geometry of the domain.

Figure 4 shows how the zone between the platen superheaters section was divided into three zones. Such an approach enabled buffer passage of grid control cells from smaller volumes in the zone within the exchangers to larger volumes in the area between the superheaters. The division mentioned above also allowed one to reduce the negative impact of narrow passages within the superheater and reheater on the effectiveness of the grid generation algorithm.

The most complex computing domain was the wind box (Figure 5), having numerous cylindrical outlets for secondary air nozzles, feeders, and burners. Direct contact of arches and all kinds of curves with cuboid parts of the geometry significantly disturbs the grid structure, which may adversely affect the quality of calculations. Therefore, a possibly large number of orthogonal sections were separated from the wind box. In order to avoid deformations of the mesh at the edges of contact of the burners with the plane of the wind box wall, the gas supply pipes were not considered. Such an approach assured that the geometric structure of this part of the domain to be entirely consistent.

After preparing the spatial geometry, the data was exported to the ANSYS meshing module. Preparation for the discretization process started with the task of the required global grid settings shown in

Table 4. Global grid settings.

Defaults
Physics Preference	CFD
Solver Preference	Fluent
Relevance	0
Element Midside Nodes	Dropped
Sizing
Size Function	Proximity and Curvature
Relevance Center	Fine
Initial Size Seed	Active Assembly
Transition	Slow
Span Angle Center	Fine
Curvature Normal Angle	Default (18.0°)
Num Cells Across Gap	3
Proximity Size Function Sources	Faces and Edges
Min Size	Default (5.77020 mm)
Proximity Min Size	Default (5.77020 mm)
Max Face Size	100.0 mm
Max Tet Size	Default (1154.0 mm)
Growth Rate	Default (1.20)
Automatic Mesh-Based Defeaturing	On
Defeature Size	Default (2.88510 mm)
Minimum Edge Length	38.0 mm
Inflation
Use Automatic Inflation	None
Assembly Meshing
Method	None
Advanced
Number of CPUs for Parallel Part Meshing	16
Number of Retries	0
Rigid Body Behavior	Dimensionally Reduced
Mesh Morphing	Disabled
Triangle Surface Mesher	Program Controlled
Topology Checking	No
Pinch Tolerance	Default (5.19320 mm)
Generate Pinch on Refresh	No

.

The maximum allowed radius of curvature between adjacent tops of grid cells on curves was equal to 18°. Since all sections of spatial geometry were described by local dimensioning of cells (so-called “Body Sizing”), global grid settings in the range of acceptable control cell size were omitted in the considerations.

By setting the “Relevance Center” option (“Fine” mode), maximum quality requirements were imposed on the shape of the grid cells created in the computational domain.

In addition, in order to ensure the highest possible mesh quality with the smallest number of cells, additional local settings were used for subsequent separated sections of spatial geometry.

The size of a single grid element (understood as the radius of the sphere described on the element) was between 40 and 60 mm. Such a dense mesh is necessary to model the transport of fluidized bed particles. The defined surfaces, used when working with the Ansys Fluent solver preprocessor, are shown in Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10.

The total number of grid nodes was 21.5 million. The number of mesh elements amounted to 28 million. In order to assess the created grid, the study of critical parameters was carried out, mainly orthogonal quality and skew. Both parameters took values from 0 to 1. In the case of orthogonal quality, the result “1” means the highest possible quality, while the criterion for the grid admission to the calculation is to achieve the result for the worst cell (minimum value of orthogonal quality in the whole grid) not lower than 0.05. In the case of skew, “0” means the highest possible quality, and the grid acceptance criterion is a worst-case cell score of no more than 0.95.

The grid parameters determined when assessing its quality are listed in

Table 5. Selected parameters of the analyzed grid quality indicators.

Orthogonal Quality
Min	0.14
Max	1.0
Mean	0.95
Standard deviation	8.0779 × 10−0.02
Skew
Min	1.3057 × 10−0.10
Max	0.95
Mean	75243 × 10−0.02
Standard deviation	0.12387

.

As a part of the grid quality evaluation process, the distribution of the number of grid cells according to orthogonal quality and skew was developed [53]. The results are presented in Figure 11, Figure 12, Figure 13 and Figure 14. Diagrams in Figure 12 and Figure 14 provide detailed views of the respective diagrams where the maximum values of the “Y” axis have been reduced to enable an analysis of the distributions in the lower quality range as diagram bars.

According to the diagrams, the grid created was dominated by hexagonal cells with a high degree of orthogonality, characteristic of the highest quality structural grids. Quadrilateral cells can also be distinguished, located primarily in the wind box section, which is too complex in shape to be able to use only cubic elements.

The above mentioned tetragonal elements are used in the wind box section, but also between the highest part of the combustion chamber and the domain section reflecting the exit duct to the cyclone. In this case, the use of the tetragonal grid is aimed at creating a kind of buffer layer. It allows the grid to pass from the orthogonal part of the combustion chamber to the diagonal outlet section so that it is possible to use the grid of hexagonal elements at the exit of the computational domain. Such an approach is essential for avoiding possible calculation errors at the outlet (the so-called numerical diffusion). Figure 15, Figure 16 and Figure 17 show selected details of the generated grid.

5. Boundary Conditions

Such geometry performed in the ANSYS meshing module was exported to the ANSYS Fluent solver preprocessor. Before the calculations, all necessary solver settings were defined. The solution to the flow issue was based on the pressure field in the domain (the so-called “Pressure-Based” method), which is recommended for issues concerning relatively low fluid flow velocities. Time-averaged calculations were conducted, which allow significantly reducing the computation time, especially necessary in case of such large and complex cases. The energy equation required for non-isothermal problems was introduced. The k-omega BSL model was applied to the turbulence field solution. A sophisticated model of heat transport through radiation was applied, assuming the effect of the reactive mixture components (carbon dioxide and water) on the radiation (the so-called “Weighted Sum of Gray Gases Model”), based on the concentration of a mixture component.

The discrete phase model (DDPM) was activated to conduct the fluidized bed simulations. Appropriate modifications of transport equations for the reactive mixture (combustion) with reactions were made, taking into account heterogeneous, surface processes (release of volatiles from fuel particles and char burnout) as well as homogeneous processes (combustion in gaseous phase).

Combustion reactions were based on a model assuming close interaction with turbulence in the flow, where the effectiveness of the reaction was dependent on the mixing of reagents (the so-called “Eddy Dissipation Model”).

The following coal properties were taken into account in the calculations: lower heating value 16.7 MJ/kg, volatile 36.6%, ash 20.6%, moisture 21.7, carbon content 85%, hydrogen 10%, oxygen 4%, and nitrogen 1%.

The first stage of conversion of the coal particles was moisture release:

$$H_2{O}_{Fuel}\stackrel{\dot{r_M}}{\to }{H}_2{O}_{Vapor}$$

(75)

$$\dot{r_M}=\{\begin{array}{c}{f}_M\frac{\left(T_P-T_{evap}\right){\rho }_M{C}_{pm}}{\Delta H_M\delta t}, T_P\ge T_{evap}\\ 0, T_P<T_{evap}\end{array}$$

(76)

where: $T_{evap}—$defined evaporation temperature, ${\rho }_M$—the density of the moisture in fuel, $C_p$—specific heat of the water in fuel, $\Delta H_M$—evaporation heat of the water in fuel, $\delta t$—time step, $T_P$—the local temperature of the solid fuel, and $f_M$—mass fraction of the moisture in the fuel.

The devolatilization rate in Equations (75) and (76) was determined based on the interphase reaction kinetics. It described the process of releasing volatiles from the fuel particles to the gas phase in which homogenous combustion reactions were occurring.

$$Volatile{s}_{Fuel}\stackrel{R_{dev}}{\to }Volatile{s}_{Gas phase}$$

(77)

$$R_{dev}=A_{dev}{e}^{\left(\frac{-E_{a, dev}}{RT}\right)}·Y_{vol\_fuel}·{\rho }_{fuel},\left[\frac{kg}{m^{3}s}\right]$$

(78)

where: $Y_{vol}$—mass fraction of volatiles and ${\rho }_{pal}$—the density of the solid fuel.

Volatiles were treated as an individual, contractual molecule whose chemical formula was determined based on fuel parameters, and especially it’s composition.

The chemical formula of the volatile molecule can be found on the left-hand side of the expression (79). It was a substrate in only one homogenous reaction in the reference-boiler-state model.

$$C_{1.90}{H}_{4.69}{O}_{0.11}{N}_{0.0337}+3.02O_2\to 1.90C{O}_2+2.34H_{2}O+0.0168N_2$$

(79)

Heterogeneous (surface) combustion of the char remaining after devolatilization was described according to Equations (80) and (81):

$$Char+O_2\stackrel{R_c}{\to }C{O}_2$$

(80)

$$R_c=A_c{e}^{\left(\frac{-E_{a,c}}{RT}\right)}·Y_{c\_solid}·{\rho }_{pal},\left[\frac{kg}{m^{3}s}\right]$$

(81)

where $Y_{c\_solid}$—char mass fraction.

The domain material is defined as a reaction mixture characterized by a variable concentration of components: volatile matter released from solid fuel, nitrogen, vapor, carbon dioxide, and oxygen. The following properties of the combustible fuel particle material are determined: density 1600 kg/m³, specific heat 1800 J/(kg K), moisture evaporation temperature 400 °C, the heat of the chemical reaction absorbed by the solid phase 30%, and particle diameter 200 μm.

The following properties of the inert material were applied for calculations: density 2600 kg/m³, specific heat 880 J/(kg K), and particle diameter 200 μm.

The heat transfer between the domain and the combustion chamber walls was considered using convection and radiation models.

The planar symmetry condition allowed only half of the system to be simulated. According to this, some other parameters related to symmetry were modified (mass flow and hydraulic diameter of primary air).

The boundary and initial conditions were as follows:

• Primary air: mass flow rate 100 kg/s, temperature 243 °C, hydraulic inlet diameter 5.1337 m, and turbulence intensity 10%,
• Four secondary air inlets (I stage): hydraulic diameter 0.285 m, mass flow rate of 2.5 kg/s, temperature 251 °C, and turbulence intensity 5%,
• Four secondary air inlets (II stage, start-up burners): hydraulic diameter 0.56 m, mass flow rate 2.3 kg/s, temperature 251 °C, and turbulence intensity 5%,
• Seven secondary air inlets (III stage): hydraulic diameter 0.245 m, mass flow rate: 1 kg/s, temperature 251 °C, and turbulence intensity 5%,
• Exit from the domain to the cyclone: hydraulic diameter 2.9981 m, pressure 300 Pa, and turbulence intensity 5%,
• Two feeders with a total fuel flow rate of 72 t/h.

The total mass of the fluidized bed is 50 tons. No limestone or any other additions were supplied.

In order to improve the convergence of the energy equation calculations, the heat transfer through external surfaces of cylindrical elements air nozzles’ outlets outside the wind box was omitted.

Due to the inclusion of gravity in the model and strong mass interactions in the system, a parallel scheme (“Coupled”) was applied. All the governing equations in the model were solved based on the “Upwind-Second Order” spatial discretization scheme.

Calculations were performed for a boiler operating at maximum capacity. After solving the reference issue described by the above settings, a variant analysis was carried out consisting of a gradual replacement of the secondary air streams with syngas (Table 3).

The previously described model has been extended with additional components of the reactive mixture, to consider the new cases with the syngas employment.

The following additional reactions have also been introduced:

$$H_2+0.5O_2\to H_{2}O$$

(82)

$$CO+O_2\to C{O}_2$$

(83)

$$2C{H}_4+3O_2\to 2CO+4H_{2}O$$

(84)

All the homogenous gas-phase reactions have been considered as the fast chemistry processes. Reaction rates were calculated based on the Eddy Dissipation approach.

The mass flow rate of the syngas was the same as the previously supplied secondary air mass flow rate. At the same time, the mass flow of solid fuel (coal) in the considered case was reduced by the equivalent of chemical energy supplied to the combustion chamber with the syngas. The decrease in fuel mass in the system was complemented by the same increase in the mass of inert material.

The following three variants of gas supply to the combustion chamber were analyzed (Figure 18):

• The use of nozzle no. 1, a total 4.6 kg/s of syngas was supplied through two side start-up burners (variant “K1”),
• Simultaneous use of nozzles No. 1 and No. 2—a total of 9.2 kg/s of syngas was supplied to the combustion chamber through four side start-up burners (variant “K2”),
• Simultaneous use of nozzles No. 1, No. 2, and No. 3—a total of 13.8 kg/s of syngas was supplied to the combustion chamber through four side start-up burner and two front burners.

The reference variant, corresponding to the conventional, monocombustion of bituminous coal, is marked with the symbol “K0”.

The syngas mas flow rate in each case was equal to 2.3 kg/s, which corresponds to 14.01 MJ of energy from gas, which corresponds to 0.839 kg/s of the considered coal.

6. Validation

Based on reports from the boiler thermal measurement, which have been shared with authors by the boiler operator, it was possible to compare data obtained based on the numerical simulation with measurement results. Figure 19 shows the comparison between the measured and calculated pressures and temperatures.

The developed CFD model was successfully validated against experimental data [54]. The good performance of the developed CFD model was achieved. The maximum relative error was lower than 10% for pressures and 5% for temperatures. Such low relative errors form a solid basis for the possibility of using the developed model in practice.

The simulation results are described in Part 2 of the paper.

7. Conclusions

An interesting idea of multi-fuel combustion in a circulating fluidized bed was considered in the paper. A comprehensive 3D CFD model of bituminous coal and syngas co-combustion in an existing large-scale OFz-425 CFB boiler was proposed. The model considers complex processes that occur in the furnace of the boiler, including dynamics of the fluidized bed, reactions, and heat transfer inside the boiler.

Four different operating scenarios were considered, including the reference variant, corresponding to the conventional, monocombustion of bituminous coal, and three tests involving the replacement of secondary air and part of the coal stream with syngas fed by start-up burners.

The model was successfully validated on the experimental results. The maximum relative error between measured and calculated data was lower than ±10%.

The detailed results of the simulations were described in part II of this paper.