CFD Simulation of Biomass Combustion in an Industrial Circulating Fluidized Bed Furnace

ABSTRACT

In this study, a three-dimensional computational fluid dynamics (CFD) model is employed to investigate the hydrodynamic and combustion characteristics of biomass particles in an industrial-scale circulating fluidized bed (CFB) furnace. The CFD model considered here is based on the Eulerian-Lagrangian framework, the multi-phase particle-in-cell (MP-PIC) collision model, the coarse grain method (CGM), and a recently developed distribution kernel method (DKM). The challenge of simulating industrial-scale CFB furnaces using CFD lies in the large number of particles in the system. MP-PIC and CGM showed that local particle overloading could occur, causing the numerical simulation to diverge. The combination of MP-PIC with CGM and DKM was shown to overcome this problem. The CFD predictions werecompared with onsite temperature experiments in the furnace, and the predicted furnace temperature agreed fairly well with the measured data. Using the CFD results, the study analyzed the transient solids mixing and fluidization characteristics, as well as the thermochemical process in biomass combustion. The simulated individual particle provided insight into the physical and chemical processes of the granular flow in the dilute/dense regions of the CFB furnace. The simulated results revealed the CO and NOx emission processes in the furnace.

KEYWORDS:

• Hydrodynamic
• combustion
• CFD simulation
• industrial-scale CFB
• formation and emission of pollutants

Introduction

Biomass, which refers to all organic materials such as wood, agricultural residues, forestry residues, and energy crops, is a renewable energy source with the potential to substitute fossil fuels in the future (Chen et al. 2021; Di Blasi 2009). Direct burning of biomass is an important energy conversion technology for generating heat and power. Biomass combustion in circulating fluidized bed (CFB) furnaces is gaining attention due to its stable low-temperature combustion, high combustion efficiency, high fuel flexibility, and low environmental impact (Khan et al. 2009; Deng et al. 2021; Cai et al. 2018).

The physical and chemical processes that occur in industrial-scale CFB furnaces are complex and include granular motions, particle–particle and particle–wall collisions, heat and mass transfer, combustion and gasification of fuel particles, formation and emission of pollutants, etc. The CFB furnace is characterized by high fluidizing velocity, a large number of particles, complex flow structures, complex particle collisions, and a turbulent combustion process. Therefore, understanding the mechanisms of the transient hydrodynamic and combustion processes in CFB furnaces and developing efficient CFD simulation approaches are paramount to designing prototypes, scaling up furnaces, controlling the operating conditions, and optimizing the operating parameters.

Two main approaches, i.e., experimental measurement and numerical modeling, have been widely adopted to investigate granule-fluid flow and gasification/combustion processes in CFB furnaces. Numerous experimental studies have been conducted to investigate physical and chemical processes in industrial-scale CFB furnaces. Svensson, Johnsson, and Leckner (1996), Leckner (2017), and Johnsson et al. (2000) investigated the characterization of fluidization regimes using time-series analysis of pressure fluctuations. The method was demonstrated to be powerful in revealing the structure of cold granule-fluid flow at the macro level. Larsson et al. (2021) investigated the thermochemical conversion of solid fuels by steam gasification in different dual fluidized beds (DFBs), showing a strong correlation between the availability of active components in the reaction environment and the quality of the product gas. Kolbitsch et al. (2010) investigated H${ }_2$, CO, and CH${ }_4$ conversion processes in two different oxygen carriers in a 120 kW dual circulating fluidized bed (DCFB) reactor. A natural oxygen carrier, i.e., ilmenite, was shown to improve the conversion efficiency compared to a fabricated Ni-based oxygen carrier. Vainio et al. conducted an experiment on the fate of fuel nitrogen in the furnace of a full-scale bubbling fluidized bed boiler (Vainio et al. 2012), measuring and analyzing the main components of nitrogen species at various height levels in detail. Based on the experimental research on the large-scale fluidized bed, the experimental research focuses more on the phenomena and mechanism of the granule-fluid flow at the furnace level. Compared to the experimental approach, the simulation approach is considered to be a more efficient, economical, and robust method to investigate hydrodynamic and combustion processes at multiple space-time scales (Alobaid et al. 2021).

Table 1 summarizes the recent work on the numerical modeling of dense granule-fluid systems, predicting physical and chemical processes, e.g., gas–solid hydrodynamics, heat and transfer, and gasification/combustion. It can be seen that one-dimensional (1-D) and two-dimensional (2-D) simulations have been frequently conducted (Blaszczuk, Zylka, and Leszczynski 2016; Collado 2018; Deng et al. 2021; Dinh et al. 2021; Lu et al. 2018), taking into account relatively detailed physical and chemical processes of the gas and particle phases. The 1-D and 2-D models have the advantages of high computational efficiency, easy implementation, and flexible application. However, they only consider the variation of the physical parameters in the vertical direction of furnaces, while assuming a uniform distribution of the physical parameters along one horizontal direction (Deng et al. 2021). Three-dimensional (3-D) simulations, as shown in Table 1, have been used to study lab-scale fluidized bed furnaces, non-reactive two-phase flow, or reactive flows without taking into account fully granular collisions. However, 3-D CFD simulations of industrial-scale fluidized bed furnaces that take into account the full details of granular motions and thermochemical processes of the particle and gas phases are desirable but rarely performed. There are several challenges in 3-D modeling the solid fuel combustion process in industrial-scale CFB furnaces, such as a huge number of particles resulting in expensive computational cost, large particles in small-size grids resulting in numerical instability, and complex chemical kinetic mechanisms involved in the devolatilization, heterogeneous char reactions, and homogeneous gas-phase reactions.

Table 1. CFD simulations of FBs. D represent dimensional, H- Hydrodynamic, T- thermol, R-gasification/combustion.

Refs	Scale	D	Model	subModel	H	T	R
Gu et al. (2020)	12 MWth	3-D	E-L	MP-PIC	✓	✓	✓
Kong et al. (2020)	lab-scale	3-D	E-L	MP-PIC	✓	✓	✓
Lin et al. (2022)	1 MWth	3-D	E-L	MP-PIC	✓	✓	✓
Yang et al. (2020)	0.3 MWth	3-D	E-L	MP-PIC	✓	✓	✓
Li and Shen (2021)	lab-scale	3-D	E-E	TFM	✓	✓	✓
Cai et al. (2022)	lab-scale	3-D	E-E	TFM	✓	✓	✓
Ghadirian, et al. (2019)	lab-scale	3-D	E-E	TFM	✓	✓	✓
Wang et al. (2018)	lab-scale	3-D	E-L	MP-PIC	✓	✓	✓
Wang et al. (2017)	600 MWth	3-D	E-L	DEM	✓	✕	✕
Luo et al. (2015)	lab-scale	3-D	E-L	DEM	✓	✕	✕
Lee et al. (2022)	550 MWth	3-D	E-L	MP-PIC	✓	✕	✕
Tu and Wang (2018)	lab-scale	3-D	E-L	MP-PIC	✓	✕	✕
Kadyrov, et al. (2019)	lab-scale	3-D	E-L	MP-PIC	✓	✕	✕
Yang and Wang (2020)	lab-scale	3-D	E-L	MP-PIC	✓	✕	✕
Muhammad, Zhang, and Wang (2019)	lab-scale	3-D	E-L	MP-PIC	✓	✕	✕
Ma et al. (2017)	lab-scale	3-D	E-L	MP-PIC	✓	✕	✕
Li et al. (2021)	lab-scale	3-D	E-L	MP-PIC	✓	✕	✕
Liu et al. (2002)	lab-scale	3-D	E-E	TFM	✓	✕	✕
Lu et al. (2018)	lab-scale	2-D	E-E	TFM	✓	✕	✕
Deng et al. (2021)	350 Mw	1-D	✕	✕	✓	✓	✓

Current 3-D CFD approaches for modeling particle flow under CFB furnace conditions can be classified into two categories: the Euler–Euler approach, e.g., the two-fluid model (TFM) (Anderson and Jackson 1967), and the Eulerian-Lagrangian approach, e.g., the discrete element method (DEM) (Tsuji, Kawaguchi, and Tanaka 1993) and multi-phase particle-in-cell (MP-PIC) approaches (Gidaspow 1994). The TFM heavily relies on the constitutive or closure relations for the source terms that describe the exchange between the gas and solid phases. The solid phase is treated as a different phase when different types of particles are present (Zhou et al. 2010), making it challenging to model complex solid-phase motion and collision processes with acceptable accuracy. The DEM accurately models granular collisions and tracks individual particles; however, modeling the collisions of quadrillions of particles on the individual particle level in an industrial-scale CFB furnace is impractical. The MP-PIC approach, which uses the kinetic theory of granular flow (KTGF) (Chapman, George, and Cowling 1990) to eliminate difficulties in calculating the inter-particle interaction, is considered an efficient and suitable method for simulating industrial-scale CFBs (Snider 2001). Additionally, the computational cost of simulating large-scale CFB furnaces can be further reduced with the coarse grain method (CGM), which clusters a number of particles with similar physical locations and properties into a virtual parcel tracked in the Lagrangian framework (Hilton and Cleary 2014). Previous work on biomass gasification in a lab-scale fluidized bed employed the CGM coupling with the Eulerian-Lagrangian method (Qi et al. 2019). However, due to the complex structure of industrial CFB furnaces and the use of unstructured grids, grid refinement is necessary. The CGM approach can cause a local overloading issue where the volume of the solid exceeds its physically possible limit in some Eulerian cells, i.e., the volume of the solid is larger than that the local cell could accommodate. A recently developed distribution kernel method (DKM) (Yang et al. 2022) remedies this issue by spreading the solid volume and source terms of solid particles in the parcel to the Eulerian domain, improving the accuracy and robustness of the model.

To the best of the authors’ knowledge, few 3-D CFD simulations that take into account granular collisions and gas-solid phase coupling have been reported for investigating the hydrodynamic and combustion processes of biomass in an industrial-scale full-loop CFB furnace, cf. Table 1. In this study, a 3-D CFD model consisting of the MP-PIC collision model, CGM, and DKM, which considers gas/solid interactions, granular collisions, heat and mass transfer, radiation, and homogeneous and heterogeneous chemical reactions, was applied to investigate a biomass-fired industrial-scale CFB furnace. The objective is to improve the understanding of gas/solid two-phase flow and the thermochemical process of biomass in industrial-scale CFB furnaces.

Mathematical formulation

In the present model, the governing equations of the continuous and discrete phases involved in fluidized bed furnaces are described in the Eulerian and Lagrangian frameworks, respectively. The governing equations of the continuous and discrete phases and main sub-models involved in the 3-D model are described in the present study.

Governing equations for the continuous phase

Reynolds-averaged Navier–Stokes (RANS) approach is used to describe the mean gas flow in the FB reactors. The gas-phase governing equations consist of the Reynolds-averaged continuity, momentum, energy, and species transport equations (Zhou et al. 2010). The Reynolds averaged continuity equation is

(1) $\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\right)={\bar{S}}_{m,g},$(1)

where overbar and tilde denote Reynolds averaged and Favre averaged, respectively. ${\alpha }_g$, ${\rho }_g$, and ${\mathbf{u}}_g$ are the gas volume fraction, the gas density, and the velocity vector of the gas phase, respectively. $S_{m,g}$ represents the gas formation rate due to the thermochemical conversion of the fuel particles.

The Reynolds averaged momentum equations are

(2) $\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g{\tilde{\mathbf{u}}}_g\right)=-{\alpha }_g\mathrm{\nabla }{\bar{p}}_g+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\tau }}_g\right)+{\bar{\mathbf{S}}}_{u,g},$(2)

where $p_g$ is the gas pressure, ${\tau }_g$ is the sum of viscous stress and Reynolds stress, and ${\mathbf{S}}_{u,g}$ is the source term of momentum exchange from the solid phase.

The Reynolds averaged energy equation is

(3) $\begin{array}{c}\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g\left(\tilde{h}+\tilde{K}\right)\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\left(\tilde{h}+\tilde{K}\right)\right)={\alpha }_g\frac{\mathrm{\partial }{\bar{p}}_g}{\mathrm{\partial }t}\\ +\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\mathrm{\Gamma }}_g\mathrm{\nabla }\tilde{h}\right)+{\bar{\dot{Q}}}_r+{\bar{\dot{Q}}}_{com}+{\bar{S}}_{h,g},\end{array}$(3)

where $h$ denotes the specific enthalpy of the gas, and $K$ denotes the kinetic energy of the gas flow. ${\dot{Q}}_r$ denotes the source term due to radiative heat transfer, ${\dot{Q}}_{com}$ denotes the source term due to volatile chemical reactions, and $S_{h,g}$ denotes the source term due to thermochemical conversion of the solid fuel. Heat diffusion coefficient ${\mathrm{\Gamma }}_g$ is the sum of the molecular and turbulent heat diffusion coefficients given by

(4) ${\mathrm{\Gamma }}_g={\mathrm{\Gamma }}_l+\frac{{\mu }_t}{{\bar{\rho }}_{g}P{r}_t},$(4)

where ${\mathrm{\Gamma }}_l$ is the molecular heat diffusion coefficient, $P{r}_t$ is the turbulent Prandtl number and ${\mu }_t$ is the turbulent eddy viscosity.

The Reynolds averaged species transport equation is

(5) $\begin{array}{c}\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g{\tilde{Y}}_{g,k}\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g{\tilde{Y}}_{g,k}\right)=\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{D}_g\mathrm{\nabla }{\tilde{Y}}_{g,k}\right)+{\bar{\dot{\omega }}}_{g,k}+{\bar{S}}_{Y_{g,k}},\end{array}$(5)

in which $Y_{g,k}$ is the mass fraction of species $k$ in the gas mixture, and ${\dot{\omega }}_{g,k}$ denotes the chemical reaction rate of species $k$. $S_{Y_{g,k}}$ denotes the formation rate of species $k$ due to thermochemical conversion of solid fuel particles. The mass diffusion coefficient $D_g$ for species $k$ taking both the molecular and turbulent contributions into account and is given by

(6) $D_g=D_l+\frac{{\mu }_t}{{\bar{\rho }}_{g}S{c}_t},$(6)

where $D_l$ is the mass diffusion coefficient for species $k$ due to molecular contribution and $S{c}_t$ is the turbulent Schmidt number.

A Partially Stirred Reactor (PaSR) model is used to account for turbulence–chemistry interaction when computing the mean source terms due to gas-phase chemical reactions (${\bar{\dot{\omega }}}_{g,k}$, ${\bar{\dot{Q}}}_{com}$) (Chomiak and Karlsson 1996). In the PaSR model, the mean reaction rates are modeled as follows:

(7) ${\bar{\dot{\omega }}}_{g,k}=\kappa {\dot{\omega }}_{g,k}(\tilde{Y},\tilde{T},\bar{p}),$(7)

in which $T$ is the gas temperature, and $\kappa$ is the volume fraction of the reactive mixture and given by,

(8) $\kappa =\frac{{\tau }_c}{{\tau }_c+{\tau }_m},$(8)

where ${\tau }_c$ and ${\tau }_m$ denote the local chemical reaction time and the local mixing time, respectively. The chemical reaction time, ${\tau }_c$, is determined from the mean reaction rates of the fuel ${\dot{\omega }}_f(\tilde{Y},\tilde{T},\bar{p})$ and the oxidizer or the gasification agents ${\dot{\omega }}_o(\tilde{Y},\tilde{T},\bar{p})$,

(9) $\frac{1}{{\tau }_c}=max\left\{\frac{-{\dot{\omega }}_f}{Y_f},\frac{-{\dot{\omega }}_o}{Y_o}\right\},$(9)

where subscripts $f$ and $o$ denote the fuel and oxidizer or the gasification agents, respectively. The mixing time ${\tau }_m$ is modeled as

(10) ${\tau }_m=C_{mix}\sqrt{\frac{\nu }{\epsilon }},$(10)

where $C_{mix}$ is a model constant ($C_{mix}=1.0$ in this study). $\nu$ and $\epsilon$ denote the kinematic viscosity and the dissipation rate of turbulent kinetic energy, respectively.

The stress tensor ${\tau }_g$ in Eq. (2) is the sum of the viscous and Reynolds stresses and can be written as

(11) ${\bar{\tau }}_g={\bar{\tau }}_l+{\bar{\tau }}_t.$(11)

The stress tensor for a Newtonian fluid ${\tau }_l$ is expressed as

(12) ${\bar{\tau }}_l={\mu }_g((\mathrm{\nabla }{\tilde{\mathbf{u}}}_g)+(\mathrm{\nabla }{\tilde{\mathbf{u}}}_g{)}^T-\frac{2}{3}(\mathrm{\nabla }\cdot {\tilde{\mathbf{u}}}_g)\mathbf{I}),$(12)

and the Reynolds stress ${\tau }_t$ is modeled according to

(13) ${\bar{\tau }}_t={\mu }_t((\mathrm{\nabla }{\tilde{\mathbf{u}}}_g)+(\mathrm{\nabla }{\tilde{\mathbf{u}}}_g{)}^T-\frac{2}{3}(\mathrm{\nabla }\cdot {\tilde{\mathbf{u}}}_g)\mathbf{I})-\frac{2}{3}{\bar{\rho }}_{g}k\mathbf{I},$(13)

Standard $k-\epsilon$ model is employed to determine the eddy viscosity,

(14) ${\mu }_t={\bar{\rho }}_g\frac{C_{\mu }{k}^2}{\epsilon },$(14)

where $k$ is the turbulent kinetic energy. $k$ and $\epsilon$ are modeled using the following transport equations:

(15) $\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_{g}k\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_{g}k\right)=\mathrm{\nabla }\cdot \left({\alpha }_g({\mu }_g+\frac{{\mu }_t}{{\sigma }_k})\mathrm{\nabla }k\right)+{\alpha }_g{P}_k-{\alpha }_g{\bar{\rho }}_g\epsilon ,$(15)

(16) $\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g\epsilon \right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\epsilon \right)=\mathrm{\nabla }\cdot \left({\alpha }_g({\mu }_g+\frac{{\mu }_t}{{\sigma }_{\epsilon }})\mathrm{\nabla }\epsilon \right)+{\alpha }_g\frac{\epsilon }{k}(C_{{\epsilon }_1}{P}_k-C_{{\epsilon }_2}{\bar{\rho }}_g\epsilon ),$(16)

where $P_k={\bar{\tau }}_t:\mathrm{\nabla }{\tilde{\mathbf{u}}}_g$ is the production rate of turbulent kinetic energy. Standard values of model constants are used, $C_{\mu }$ = 0.09, $C_{\epsilon 1}$ = 1.44, $C_{\epsilon 2}$ = 1.92, $C_{\sigma k}$ = 1.0 and $C_{\sigma \epsilon }$ = 1.3 (Ku, Li, and Løvås 2015; Yan et al. 2016). The mean source terms are due to the particle/gas interaction in Equations 1(1) $\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\right)={\bar{S}}_{m,g},$(1) , 2(2) $\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g{\tilde{\mathbf{u}}}_g\right)=-{\alpha }_g\mathrm{\nabla }{\bar{p}}_g+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\tau }}_g\right)+{\bar{\mathbf{S}}}_{u,g},$(2) , 3(3) $\begin{array}{c}\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g\left(\tilde{h}+\tilde{K}\right)\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g\left(\tilde{h}+\tilde{K}\right)\right)={\alpha }_g\frac{\mathrm{\partial }{\bar{p}}_g}{\mathrm{\partial }t}\\ +\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\mathrm{\Gamma }}_g\mathrm{\nabla }\tilde{h}\right)+{\bar{\dot{Q}}}_r+{\bar{\dot{Q}}}_{com}+{\bar{S}}_{h,g},\end{array}$(3) and 5(5) $\begin{array}{c}\frac{\mathrm{\partial }\left({\alpha }_g{\bar{\rho }}_g{\tilde{Y}}_{g,k}\right)}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{\tilde{\mathbf{u}}}_g{\tilde{Y}}_{g,k}\right)=\mathrm{\nabla }\cdot \left({\alpha }_g{\bar{\rho }}_g{D}_g\mathrm{\nabla }{\tilde{Y}}_{g,k}\right)+{\bar{\dot{\omega }}}_{g,k}+{\bar{S}}_{Y_{g,k}},\end{array}$(5) ), i.e., ${\bar{S}}_{m,g}$, ${\bar{\mathbf{S}}}_{u,g}$, ${\bar{S}}_{h,g}$ and ${\bar{S}}_{Y_{g,k}}$, require the modeling of particle phase as discussed below.

Solid phase governing equations

In the Eulerian-Lagrangian approach, biomass and sand particles are tracked using the Lagrangian approach. The interactions between the particles and the surrounding gas are through mass and momentum exchange and heat transfer. The mass, momentum, and energy conservation equations for the solid phase in the Lagrangian framework are presented in the following. For simplicity, the Reynolds/Favre averaged gas properties are indicated without using over-bars or tildes.

Mass conversion of solid phase

Biomass particles undergo thermochemical conversion reactions, i.e., drying, pyrolysis, and the heterogeneous reaction of char, while sand particles are assumed to be chemically inert. The mass conservation equation for the $i$-th biomass particle is written as

(17) $\frac{d{m}_i}{dt}={\dot{m}}_{vapor,i}+{\dot{m}}_{devol,i}+{\dot{m}}_{char,i},$(17)

where $m_i$, ${\dot{m}}_{vapor,i}$, ${\dot{m}}_{devol,i}$ and ${\dot{m}}_{char,i}$ denote the mass of $i$-th biomass particle, the evaporation rate, the devolatilization rate, and the char conversion rate, respectively.

Drying

The moisture evaporation rate (Ku et al. 2014; Yan et al. 2016) is modeled as

(18) ${\dot{m}}_{vapor,i}=-{\varphi }_{vapor,i}A{s}_i{M}_v,$(18)

where ${\varphi }_{vapor,i}$, $A{s}_i$, and $M_v$ represent the molar flux of vapor, the surface area of the particle, and the molar weight of the vapor, respectively. ${\varphi }_{vapor,i}$ is given by

(19) ${\varphi }_{vapor,i}=k_c(C_{vapor,i}-C_{vapor,g}),$(19)

where $k_c$, $C_{vapor,i}$ and $C_{vapor,g}$ denote the mass transfer coefficient, vapor concentration at the particle surface, and the vapor concentration in the bulk gas, respectively. $k_c$, $C_{vapor,i}$ and $C_{vapor,g}$ can be described as

(20) $k_c=\frac{Sh{D}_{diff,va}}{d_i},$(20)

(21) $C_{vapor,i}=\frac{P_{sat,T_i}}{R_u{T}_i},$(21)

and

(22) $C_{vapor,g}=X_v\frac{p_g}{R_u{T}_g},$(22)

where $Sh$ is the Sherwood number modeled using Ranz–Marshall correlation (Ranz and Marshall 1952)

(23) $Sh=(2+0.6R{e}_i^{1/2}S{c}^{1/3}),$(23)

where $Sc$ is the Schmidt number of the surrounding gas and $R{e}_i$ is the Reynolds number of $i$th particle. $D_{diff,va}$, $P_{sat,T_i}$, $T_g$, and $X_v$ represent the vapor diffusion coefficient, the saturation pressure, the gas temperature, and the molar fraction of vapor in the surrounding gas, respectively. $R_u$ is the universal gas constant, $T_i$ is the particle temperature, and $d_i$ is an equivalent spherical particle diameter computed based on the particle real-time mass $m_i$ and a constant particle density ${\rho }_i$ and is given by

(24) $d_i=(6m_i/\pi {\rho }_i{)}^{1/3},$(24)

where ${\rho }_i$ is the particle density.

Pyrolysis

There are four different types of pyrolysis models reviewed by Hameed et al. (2019), Vikram, Rosha, and Kumar (2021), and Fatehi et al. (2021), i.e., (a) single-step model, (b) three-parallel-step model with secondary tar cracking reactions, (c) FG-DVC model, and (d) multicomponent pyrolysis model. In fluidized bed furnaces, especially large-scale industrial furnaces, the number of biomass particles is enormous. It would require tremendously long computational time to carry out numerical simulations of a 3-D fluidized bed furnace if multicomponent pyrolysis models were used. Three-parallel step model and FG-DVC model have some drawbacks, as they are developed based on specific experimental conditions. Therefore, using them to predict other experimental conditions may result in substantial experimental errors. Thus, single-step models have been employed in fluidized bed simulations (Gómez et al. 2014; Karim and Naser 2018; Ku, Li, and Løvås 2015; Luo et al. 2022, 2020; Qi et al. 2019; Wang et al. 2018; Yang et al. 2022; Yang et al. 2019; Zhou et al. 2006).

The rate of devolatilization is computed based on the pyrolysis reaction model,

(25) ${\dot{m}}_{devol,i}=-A_{d}exp\left(-\frac{E_d}{R_u{T}_i}\right)m_{volat,i},$(25)

where $A_d$ and $E_d$ are rate constants (Ku, Li, and Løvås 2015), $m_{volat,i}$ is the mass of the volatile remaining in the particle.

A one-step pyrolysis model involving nitrogen conversion is written as

(26) $\begin{array}{c}Biomass\to x_1\mathrm{C}\mathrm{O}+x_2\mathrm{C}\mathrm{O}{}_2+x_3\mathrm{H}{}_2\mathrm{O}+x_4\mathrm{H}{}_2+x_5\mathrm{C}\mathrm{H}{}_4\\ +x_6\mathrm{H}\mathrm{N}\mathrm{C}\mathrm{O}+x_7\mathrm{H}\mathrm{C}\mathrm{N}+x_8\mathrm{N}\mathrm{H}{}_3+x_9\mathrm{N}\mathrm{O}+x_{10}\mathrm{C}{}_{(s)},\end{array}$(26)

where C${ }_{(s)}$ denotes char in the solid phase. $x_j$ are the stoichiometric constants, i.e., $x_1$ = 0.5014, $x_2$ = 0.0954, $x_3$ = 0.0864, $x_4$ = 0.0512, $x_5$ = 0.1060, $x_6$ = 0.0021, $x_7$ = 0.0043, $x_8$ = 0.0067, $x_9$ = 0.0005, and $x_{10}$ = 0.1458. In this model, volatile nitrogen-containing species released during the pyrolysis process include NH${ }_3$, NO, HCN, and HNCO, whose release rates are proportional to biomass pyrolysis.

Biomass NO${ }_{\mathrm{x}}$ formation mechanism has been investigated for several decades. Winter, Wartha, and Hofbauer (1999) investigated the NO${ }_{\mathrm{x}}$ formation of different biomass fuels in a fluidized bed combustor and a grate furnace. NO, N${ }_2$O, HCN, and NH${ }_3$ were measured in the flue gas shortly after biomass combustion, while N${ }_2$O was rapidly converted to N${ }_2$. HCN was formed in quantities similar to NH${ }_3$ during woody biomass combustion and the HCN/NH${ }_3$ ratios depend on the H/N ratio in biomass fuels. According to the measurements of Bassilakis, Carangelo, and Wojtowicz (2001) and Hansson et al. (2004), HNCO is a significant intermediate product for NO${ }_{\mathrm{x}}$ formation during biomass combustion. In the study of Bassilakis, Carangelo, and Wojtowicz (2001), the mass ratios (dry basis) of NH${ }_3$/HCN/HNCO at a heating rate of 30 K/min are 37/43/20 for wheat straw and 35/26/39 for tobacco, respectively. Hansson et al. (2004) reported that the mass ratios (dry basis) of NH${ }_3$/HCN/HNCO are 57/28/15 at a pyrolysis temperature 973 K and 31/60/9 at 1273 K, respectively. According to studies by Leppalahti and Koljonen (1995) and Weissinger, Fleckl, and Obernberger (2004), NH${ }_3$ is the main nitrogen-containing intermediate product during biomass pyrolysis. Zhou et al. (2006) showed that up to 1–4% of nitrogen is directly converted to NO during biomass pyrolysis. Despite numerous studies on NO${ }_{\mathrm{x}}$ formation in biomass combustion, there is no general consensus on the ratio of NH${ }_3$/HCN/HNCO/NO in the published literature. Based on the above studies, the components containing nitrogen in the pyrolysis products are NH${ }_3$, HCN, HNCO, and NO in descending order. NH${ }_3$ is the main component of pyrolysis nitrogenous products, while the ratio of HCN/HNCO is approximately 2. The mass ratio of NH${ }_3$/HCN/HNCO/NO during biomass pyrolysis is estimated to be 51/31/15/3 in the present work. This ratio is used to determine the model constants $x_j$ in Equation (26)(26) $\begin{array}{c}Biomass\to x_1\mathrm{C}\mathrm{O}+x_2\mathrm{C}\mathrm{O}{}_2+x_3\mathrm{H}{}_2\mathrm{O}+x_4\mathrm{H}{}_2+x_5\mathrm{C}\mathrm{H}{}_4\\ +x_6\mathrm{H}\mathrm{N}\mathrm{C}\mathrm{O}+x_7\mathrm{H}\mathrm{C}\mathrm{N}+x_8\mathrm{N}\mathrm{H}{}_3+x_9\mathrm{N}\mathrm{O}+x_{10}\mathrm{C}{}_{(s)},\end{array}$(26) .

Char conversion

Char conversion is a complex process in which chemical reactions occur at the surface of the porous medium structure with complex interior and microstructures. The heterogeneous rates of char conversion are affected by the fundamental components, e.g., surface area, surface accessibility, carbon active sites, added inorganic matter, and the gaseous reactant concentration (Di Blasi 2009). The rate of char conversion is computed based on all heterogeneous reactions,

(27) ${\dot{m}}_{char,i}=-\sum _{j=1}^3{\dot{m}}_{char,ij},$(27)

where ${\dot{m}}_{char,ij}$ represent the char consumption rates by reactions with O${ }_2$, H${ }_2$O, and CO${ }_2$, respectively.

(28) ${\dot{m}}_{char,ij}=-A{s}_i{p}_j{R}_{d,j}Da{{ }^{\mathrm{\prime }}}_j,$(28)

where $A{s}_i$ denotes the particle surface, and $p_j$ represents the partial pressure of the gasifying agents or oxidizers in the gas surrounding the particle. A normalized Damköhler number $Da{{ }^{\mathrm{\prime }}}_j$, which is the ratio of the kinetic reaction rate to the mass transport rate (Hazenberg and van Oijen 2021), is defined to take into account the contribution of the kinetic and the diffusion rates,

(29) $Da{{ }^{\mathrm{\prime }}}_j=-\frac{R_{kin,j}}{R_{d,j}+R_{kin,j}},$(29)

where $R_{d,j}$ and $R_{kin,j}$ represent, respectively, the diffusion rate coefficient and kinetic rate coefficient. $R_{d,j}$ and $R_{kin,j}$ are defined as follows:

(30) $R_{d,j}=C_j\frac{{[0.5(T_g+T_i)]}^{0.75}}{d_i},$(30)

and

(31) $R_{kin,j}=A_{j}exp\left(-\frac{E_j}{R_u{T}_i}\right),$(31)

where $A_j$ and $E_j$ represent the pre-exponential factor and activation energy for the char gasification reactions, respectively. $C_j$ is the mass diffusion rate constant and $C_j=5\times {10}^{-12}$ (s/K${ }^{0.75}$) (Ku, Li, and Løvås 2015).

Simplified homogeneous reactions of volatile gas and heterogeneous reactions of char used in this study are listed in Table 2. In this model, thermal NO formation is neglected because the maximum temperature in the furnace is lower than 1600 K, i.e., fuel-NO${ }_{\mathrm{x}}$ from nitrogen in the biomass is the main source of NO${ }_{\mathrm{x}}$ formation. This chemical kinetic model is selected mainly due to its high computational efficiency. The NO chemistry (R9–R16) has been used by H. Zhou et al. (2006) to predict NO formation in straw combustion in a fixed bed furnace showing good accuracy.

Table 2. Homogeneous and heterogeneous reactions considered in biomass combustion and gasification. Note: C${ }_{(s)}$ is solid phase char. $C_k$ represents the molar concentration of gas species $k$. References for each reaction: R1 (Yan et al. 2016, 2018), R2 (Yan et al. 2016, 2018), R3 (Yan et al. 2016, 2018), R4 (Yan et al. 2016, 2018), R5 (Yan et al. 2016, 2018), R6 (Yan et al. 2016, 2018), R7 (Brink, Kilpinen, and Hupa 2001), R8 (Brink, Kilpinen, and Hupa 2001), R9 (Zhou et al. 2006), R10 (Zhou et al. 2006; Ma et al. 2021), R11 (Zhou et al. 2006; Ma et al. 2021), R12 (Ma et al. 2021; Zhou et al. 2006), R13 (Ma et al. 2021; Zhou et al. 2006), R14 (Ma et al. 2021; Zhou et al. 2006), R15 (Ma et al. 2021; Zhou et al. 2006), R16 (Ma et al. 2021; Zhou et al. 2006), R17 (Nikoo and Mahinpey 2008; Yang et al. 2019), R18 (Nikoo and Mahinpey 2008; Yang et al. 2019), R19 (Nikoo and Mahinpey 2008; Yang et al. 2019.).

Reference	Reactions	Reaction rate [Kmol/m${ }^3$/s]
R1	CH${ }_4$ + H${ }_2$O $\to$ CO + 3H${ }_2$	$0.312exp(\frac{-15,098}{T_g})C_{C{H}_4}{C}_{H_{2}O}$
R2	CO + H${ }_2$O $\to$ CO${ }_2$ + H${ }_2$	$2.5\times {10}^{8}exp(\frac{-16,597}{T_g})C_{CO}{C}_{H_{2}O}$
R3	CO${ }_2$ + H${ }_2$ $\to$ CO + H${ }_2$O	$9.43\times {10}^{9}exp(\frac{-20,563}{T_g})C_{C{O}_2}{C}_{H_2}$
R4	CH${ }_4$ + 2O${ }_2$ $\to$ CO${ }_2$ + 2H${ }_2$O	$2.119\times {10}^{11}exp(\frac{-24,379}{T_g})C_{C{H}_4}^{0.2}{C}_{O_2}^{1.3}$
R5	CO + 0.5O${ }_2$ $\to$ CO${ }_2$	$1.0\times {10}^{10}exp(\frac{-15,154}{T_g})C_{CO}{C}_{O_2}^{0.5}$
R6	H${ }_2$ + 0.5O${ }_2$ $\to$ H${ }_2$O	$2.2\times {10}^{9}exp(\frac{-13,109}{T_g})C_{H_2}{C}_{O_2}$
R7	NH${ }_3$ + O${ }_2$ $\to$ NO + H${ }_2$O + 0.5H${ }_2$	$1.21\times {10}^{9}exp(\frac{-8,000}{T_g})C_{O_2}{C}_{N{H}_3}$
R8	NH${ }_3$ + NO $\to$ N${ }_2$ + H${ }_2$O + 0.5H${ }_2$	$8.73\times {10}^{17}exp(\frac{-8,000}{T_g})C_{NO}{C}_{N{H}_3}$
R9	CH${ }_4$ + NO $\to$ HCN + H${ }_2$O + 0.5H${ }_2$	$2.7\times {10}^{6}exp(\frac{-9,466}{T_g})C_{NO}{C}_{C{H}_4}$
R10	CH${ }_4$ + NO $\to$ CO + 2H${ }_2$ + 0.5N${ }_2$	$2.7\times {10}^{6}exp(\frac{-9,466}{T_g})C_{NO}{C}_{C{H}_4}$
R11	HCN + 0.5O${ }_2$ $\to$ CNO + 0.5H${ }_2$	$2.14\times {10}^{5}exp(\frac{-10,000}{T_g})C_{O_2}{C}_{HCN}$
R12	CNO + 0.5O${ }_2$ $\to$ CO + NO	$2.14\times {10}^{5}exp(\frac{-10,000}{T_g})C_{O_2}{C}_{HCN}$
R13	CNO + NO $\to$ N${ }_2$ + 0.5O${ }_2$ + CO	$1.02\times {10}^{9}exp(\frac{-25,460}{T_g})C_{O_2}{C}_{HCN}$
R14	HNCO $\to$ NCO + 0.5H${ }_2$	$2.14\times {10}^{5}exp(\frac{-10,000}{T_g})C_{O_2}{C}_{HNCO}$
R15	NCO + 0.5O${ }_2$ $\to$ NO + CO	$2.14\times {10}^{5}exp(\frac{-10,000}{T_g})C_{O_2}{C}_{HNCO}$
R16	NCO + NO $\to$ N${ }_2$ + 0.5O${ }_2$ + CO	$1.02\times {10}^{9}exp(\frac{-25,646}{T_g})C_{O_2}{C}_{HNCO}$
R17	C${ }_{(s)}$ + 0.5O${ }_2$ $\to$ CO	$0.046\times {10}^{7}exp(\frac{-13,523}{R_u{T}_i})$
R18	C${ }_{(s)}$ + H${ }_2$O $\to$ CO + H${ }_2$	$1.71\times {10}^{7}exp(\frac{-211,000}{R_u{T}_i})$
R19	C${ }_{(s)}$ + CO${ }_2$ $\to$ 2CO	$9.1\times {10}^{6}exp(\frac{-166,00}{R_u{T}_i})$

Momentum equation of solid phase

The velocity of the $i$-th particle ${\mathbf{u}}_i$ is governed by Newton’s second law,

(32) $m_i\frac{d{\mathbf{u}}_i}{dt}={\mathbf{f}}_{d,i}+{\mathbf{f}}_{{\mathrm{\nabla }}_p,i}+m_i\mathbf{g}+{\mathbf{f}}_{\tau ,i}.$(32)

The right-hand side terms represent the sum of all forces acting on the $i$-th particle by the surrounding gas and particles. The forces considered include, from left to right, the drag ${\mathbf{f}}_{d,i}$, pressure gradient ${\mathbf{f}}_{{\mathrm{\nabla }}_p,i}$, gravity $m_i\mathbf{g}$, and interparticle stress ${\mathbf{f}}_{\tau ,i}$. With a given ${\mathbf{u}}_i$, the position vector of the particle ${\mathbf{x}}_i$ is computed by integration of the equation

(33) $d{\mathbf{x}}_i/dt={\mathbf{u}}_i.$(33)

Drag model

The drag force model widely used for the $i$-th individual particle ${\mathbf{f}}_{d,i}$ is given by (Gidaspow 1994; Ku, Li, and Løvås 2015; Yang et al. 2019)

(34) ${\mathbf{f}}_{d,i}=V_{\mathrm{\Omega }}\beta ({\mathbf{u}}_g-{\mathbf{u}}_i),$(34)

where $V_{\mathrm{\Omega }}$ is the volume of the computational cell and $\beta$ is the drag force parameter, which is modeled using the Wen & Yu drag correlation (Gidaspow 1994; Wen 1966) and is given as

(35) $\beta =\left(\begin{array}{c}150\frac{{(1-{\alpha }_g)}^2{\mu }_g}{{\alpha }_g^2{d}_i^2}+1.75\frac{(1-{\alpha }_g){\rho }_g}{{\alpha }_g{d}_i}|{\mathbf{u}}_g-{\mathbf{u}}_i|{\alpha }_g<0.8\\ \frac{3}{4}{C}_d\frac{(1-{\alpha }_g){\rho }_g}{d_i}|{\mathbf{u}}_g-{\mathbf{u}}_i|{\alpha }_g^{-2.65}{\alpha }_g\ge 0.8\end{array}\right.,$(35)

where ${\alpha }_g$ is the gas-phase volume fraction, and the drag coefficient $C_d$ is modeled as (Gidaspow 1994)

(36) $C_d=\left(\begin{array}{c}\frac{24}{R{e}_i}(1+0.15R{e}_i^{0.687})R{e}_i<1000\\ 0.44R{e}_i\ge 1000\end{array}\right.,$(36)

where the particle Reynolds number $R{e}_i$ is defined as

(37) $R{e}_i={\alpha }_g{\rho }_g{d}_i|{\mathbf{u}}_g-{\mathbf{u}}_i|/{\mu }_g.$(37)

Interparticle stress

The particle stress ${\mathbf{f}}_{\tau ,i}$ is given by

(38) ${\mathbf{f}}_{\tau ,i}=V_{\mathrm{\Omega }}\mathrm{\nabla }\tau ,$(38)

where the contact normal stress $\tau$ can be given by the model of Lun et al. (1984),

(39) $\tau =[{\theta }_s{\bar{\rho }}_s+{\theta }_s^2{\bar{\rho }}_s(1+e)g_0]{\mathrm{\Theta }}_s,$(39)

where $g_0$, ${\bar{\rho }}_s$, and $e$ represent, respectively, the radial distribution function, the mean density of particles in a local cell, and the coefficient of restitution. ${\mathrm{\Theta }}_s$ is the granular temperature, and $g_0$ is the radial distribution function.

The solid volume fraction of particles ${\theta }_s$ satisfies ${\theta }_s+{\alpha }_g=1$ and is modeled based on the particle distribution function $f(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)$ and can be given by

(40) ${\theta }_s=\int \int f\frac{m_s}{{\rho }_s}d{m}_{s}d{\mathbf{u}}_s.$(40)

where ${\mathbf{u}}_s$ and $m_s$ denote the particle velocity and particle mass, respectively. The particle velocity ${\mathbf{u}}_s$ is the particle velocity in the Eulerian frame, which is different from ${\mathbf{u}}_i$ in Equation (32)(32) $m_i\frac{d{\mathbf{u}}_i}{dt}={\mathbf{f}}_{d,i}+{\mathbf{f}}_{{\mathrm{\nabla }}_p,i}+m_i\mathbf{g}+{\mathbf{f}}_{\tau ,i}.$(32) that represents the velocity of the $i$-th particle in the Lagrangian framework. This also applies to the distinctions between solid-phase parameters with subscripts of $s$ and $i$.

In the MP-PIC model, $f$ is obtained from the Liouville equation, which is a mathematical expression of the conservation of particle numbers per volume moving along dynamic trajectories in the particle-phase space (Andrews and O’Rourke 1996),

(41) $\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(f{\mathbf{u}}_s\right)+{\mathrm{\nabla }}_u\cdot \left(f{\mathbf{A}}_s\right)=\frac{f_G-f}{{\tau }_G}+\frac{f_D-f}{{\tau }_D}.$(41)

The first term in the RHS of Eq. (41) denotes the collision return-to-isotropy effect and the second term denotes the collision damping effect. Physically, particle collision tends to dampen out the velocity fluctuations. The collision model assumed that within a damping relaxation time, the particle velocity approaches a mean value and the distribution function $f(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)$ approaches $f_D(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)$. The collision-damping particle distribution function $f_D(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)$ is given by (O’Rourke and Snider 2010, 2012)

(42) $f_D(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)=\delta ({\mathbf{u}}_s-{\bar{\mathbf{u}}}_s)\int fd{\mathbf{u}}_s,$(42)

where $\delta$ is the Dirac function. The mean value of particle velocity can be given by

(43) ${\bar{\mathbf{u}}}_s=\frac{\int \int f{m}_s{\mathbf{u}}_{s}d{m}_{s}d{\mathbf{u}}_s}{\int \int f{m}_{s}d{m}_{s}d{\mathbf{u}}_s}.$(43)

The particle collision could result in a Gaussian distribution of particle velocity occurring within a relaxation time ${\tau }_G$. The Gaussian distribution is described by the equilibrium-isotropic particle distribution function $f_G(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)$,

(44) $f_G(m_s,{\mathbf{u}}_s,{\mathbf{x}}_s,t)=G({\mathbf{u}}_s;{\bar{\mathbf{u}}}_s,{\sigma }^2)\int fd{\mathbf{u}}_s,$(44)

where $G$ is a Gaussian velocity distribution with the mean ${\bar{u}}_s$ and variance ${\sigma }^2$, which can be obtained by enforcing that the variance of $f_G$ is equal to that of $f$.

Energy equation of solid phase

The particle temperature is obtained from the energy conservation equation for the $i$-th particle,

(45) $\begin{array}{c}{m}_i{C}_{p,i}\frac{d{T}_i}{dt}=q_{c,i}+q_{r,i}-q_{vapor,i}+q_{devol,i}+q_{char,ij}.\end{array}$(45)

where $C_{p,i}$, $q_{c,i}$ and $q_{r,i}$ denote, respectively, the particle heat capacity and convective and radiative heat transfer. $q_{vapor,i}$, $q_{devol,i}$ and $q_{char,ij}$ represent the heat transfer of latent, pyrolysis, and char reactions.

The convection heat $q_{c,i}$ and radiation $q_{r,i}$ are given by

(46) $q_{c,i}=h_{i}A{s}_i(T_g-T_i),$(46)

and

(47) $q_{r,i}=\frac{{\epsilon }_{i}A{s}_i}{4}(G-4\sigma T_i^4),$(47)

where $h_i$, ${\epsilon }_i$, $\sigma$, and $G$ represent interphase thermal transfer coefficient, emissivity, Stefan-Boltzmann constant, and incident radiation, respectively. The interphase thermal transfer coefficient $h_i$ can be given

(48) $h_i=\frac{Nu{\lambda }_g}{d_i},$(48)

where ${\lambda }_g$ is the thermal conductivity of the surrounding gas. $Nu$ is the Nusselt number computed using the Ranz-Marshall correlation given by

(49) $Nu=2+\frac{3}{5}R{e}_i^{1/2}P{r}^{1/3},$(49)

where $Pr$ is the Prandtl number of the surrounding gas. The incident radiation $G$ is obtained from the P-1 radiation model.

The heat fluxes due to evaporation, pyrolysis, and char reactions are, respectively, $q_{vapor,i}$, $q_{devol,i}$ and $q_{char,ij}$, and can be given by

(50) $q_{vapor,i}=h_{vapor,i}{\dot{m}}_{vapor,i},$(50)

(51) $q_{devol,i}=h_{devol,i}{\dot{m}}_{devol,i},$(51)

and

(52) $q_{char,ij}=\sum _{j=1}^3{h}_{i,j}{\dot{m}}_{char,ij},$(52)

where $h_{vapor,i}$, $h_{devol,i}$, and $h_{i,j}$ represent the latent heat, the heat of pyrolysis, and the heat of char reactions, respectively.

Numerical methods and computation cases

Solution procedure for solid-phase governing equations

In order to achieve a converged statistical solution in the MP-PIC method, a sufficiently large number of particles need to be simulated. This is not an issue in fluidized bed reactors since there is a sufficiently large number of particles. In practical simulations, not all these particles could be simulated. A coarse grain method (CGM) is employed in this study to reduce the computational cost. In the CGM approach, a finite number of virtual particles (hereafter referred to as parcels) are simulated. Assume that the number of parcels is $N_p$. The $i$-th parcel contains multiple real particles; however, all particles have the same properties, i.e., each real particle in the $i$-th parcel has the same mass $m_i$, velocity ${\mathbf{u}}_i$, temperature $T_i$ and diameter $d_i$.

The governing equations for the individual real particles in the $i$-th parcel are presented in subsection 2.2. These equations are integrated to compute the particle quantities, i.e., $m_i$, ${\mathbf{u}}_i$, $T_i$, and $d_i$. Implicit backward Euler scheme is used in the temporal integration of these equations.

As an example, the velocity of the $i$-th particles is obtained by integrating Eq. (32). The discrete form of the velocity equation for the $i$-th particles can be written as (O’Rourke and Snider 2012).

(53) $\frac{{\mathbf{u}}_i^{n+1}-{\mathbf{u}}_i^n}{\mathrm{\Delta }t}=D_i^n({\mathbf{u}}_{g,i}^{n+1}-{\mathbf{u}}_i^{n+1})+S_i^{n+1}+\frac{{\bar{\mathbf{u}}}_i^{n+1}-{\mathbf{u}}_i^n}{2{\tau }_D}$(53)

where superscript $n$ denotes the quantities at time $t^n$ and $\mathrm{\Delta }t$ is the time step. The first term on the RHS is the drag term, and the second term $S_i^{n+1}$ is the sum of source terms due to the pressure gradient force, gravity, and interparticle stress. The last term is explicitly added to model the effect of the collision damping term in the Liouville equation (41)(41) $\frac{\mathrm{\partial }f}{\mathrm{\partial }t}+\mathrm{\nabla }\cdot \left(f{\mathbf{u}}_s\right)+{\mathrm{\nabla }}_u\cdot \left(f{\mathbf{A}}_s\right)=\frac{f_G-f}{{\tau }_G}+\frac{f_D-f}{{\tau }_D}.$(41) , where ${\mathbf{u}}_i$ is the mass-weighted average of $u_i$ (O’Rourke and Snider 2012).

In this discrete form of particle velocity equation, ${\mathbf{u}}_{g,i}^{n+1}$ is the gas velocity at (${\mathbf{x}}_i,t^{n+1}$), which is computed from the gas velocity in the Euler grid around the particle position ${\mathbf{x}}_i$ at time $t^{n+1}$. A trilinear interpolation scheme is used to interpolate the Eulerian field quantities defined in the Euler cells to the discrete Lagrangian particle location ${\mathbf{x}}_i$. Figure 1 illustrates the interpolation procedure.

Figure 1. Diagram of solution procedure MP-PIC.

Once the mass, velocity, temperature, and position are computed for all particles, the solid volume fraction (which is an Euler field quantity) can be computed from the ensemble average of the particles. Assume that the number of real particles per unit volume that pertains to the $i$-th parcel is $n_i$. The solid volume fraction at ($\mathbf{x},t$) is

(54) $\theta (\mathbf{x},t)=\frac{1}{V_{\mathrm{\Omega }}}\sum _{i=1}^{N_p}{n}_i{V}_i({\mathbf{x}}_i)S(\mathbf{x},{\mathbf{x}}_i),$(54)

where $V_i$ is the volume of the $i$-th particle. ${\mathbf{x}}_i$ is the location of the $i$-th particle, whereas $S(\mathbf{x},{\mathbf{x}}_i)$ is the trilinear interpolation function that computes the Euler field properties at $\mathbf{x}$ from the Lagrangian quantities at ${\mathbf{x}}_i$.

The source terms due to the gas/solid interaction for the continuity equation, momentum equations, enthalpy equation, and species transport equations are computed similarly,

(55) $\begin{array}{c}{S}_{m,g}=-\frac{1}{V_{\mathrm{\Omega }}}\sum _{i=1}^{N_p}{n}_i{\dot{m}}_{i}S(\mathbf{x},{\mathbf{x}}_i),{\mathbf{S}}_{u,g}=-\frac{1}{V_{\mathrm{\Omega }}}\sum _{i=1}^{N_p}{n}_i{\mathbf{f}}_{i}S(\mathbf{x},{\mathbf{x}}_i),\\ S_{h,g}=-\frac{1}{V_{\mathrm{\Omega }}}\sum _{i=1}^{N_p}{n}_i{q}_{i}S(x,x_i),S_{Y_{g,k}}=-\frac{1}{V_{\mathrm{\Omega }}}\sum _{i=1}^{N_p}{n}_i{\dot{m}}_{i,k}S(\mathbf{x},{\mathbf{x}}_i),\end{array}$(55)

where $q_i$ is the heat exchange rate from solid particle, ${\mathbf{f}}_i={\mathbf{f}}_{d,i}+{\mathbf{f}}_{{\mathrm{\nabla }}_p,i}$ is the sum of the drag force and the pressure gradient force, and ${\dot{m}}_{i,k}$ is due to pyrolysis and char reactions.

Figure 1 shows the solution procedure of the MP-PIC model. In this figure, $C_1$ and $C_2$ are the Euler cell centers. Black and yellow particles denote the biomass and sand particles, respectively. The solution procedure involves the following four steps (the order of steps is not the order of execution in the CFD code):

• (1) As shown in Figure 1(a), the Euler field quantities in cell $C_1$ and $C_2$, e.g., the solid volume fraction and the source terms of the gas phase equations due to the particles, are computed from the Lagrangian particles in the cells.

• (2) As shown in Figure 1(b), the gas phase transport equations are numerically solved using the finite-volume method described in subsection 2.1.

• (3) As shown in Figure 1(c), the Euler field quantities in cell $C_1$ and $C_2$, e.g., ${\mathbf{u}}_g^{n+1}$, is interpolated to the position of the $i$-th Lagrangian particles, i.e., to compute ${\mathbf{u}}_{g,i}^{n+1}$ in Equation (53)(53) $\frac{{\mathbf{u}}_i^{n+1}-{\mathbf{u}}_i^n}{\mathrm{\Delta }t}=D_i^n({\mathbf{u}}_{g,i}^{n+1}-{\mathbf{u}}_i^{n+1})+S_i^{n+1}+\frac{{\bar{\mathbf{u}}}_i^{n+1}-{\mathbf{u}}_i^n}{2{\tau }_D}$(53) .

• (4) As shown in Figure 1(d), the properties of the Lagrangian particles are computed by temporal integration of the particle phase equations, using the method described earlier.

In the CGM-PCM approach, a large number of parcels in a small cell contribute to large source terms and local overloading of solids, i.e., the solid volume fraction (${\theta }_s$) is larger than the physically allowable value, e.g., ${\theta }_s>0.62$ (Sun and Xiao 2015; Yang et al. 2022). Since ${\alpha }_g=1-{\theta }_s$, a large ${\theta }_s$ leads to a small ${\alpha }_g$. Too large source terms and too small ${\alpha }_g$ can result in numerical instability of the gas phase governing equations. Thus, a threshold in the CFD solver is often employed, e.g., when ${\theta }_s>0.62$, ${\theta }_s$ is set to the value of 0.62. The use of such a threshold can result in mass loss of the solid phase in the gas-phase governing equations in PCM (due to the increased ${\alpha }_g$). Hence, a distribution kernel method (DKM) is developed to address this issue, as well as the cell searching strategy and parallel computation method for the DKM.

Spatial redistribution of parcels and source terms

The distribution kernel method (DKM) model was proposed and validated in our previous work (Yang et al. 2022, 2023), providing a detailed description of the method. A filtering kernel function $g(x,t)$ was introduced to redistribute the source terms to surrounding cells containing the parcel based on,

(56) ${\varphi }_r(\mathbf{x},t)=g(\mathbf{x},t){\int }_{\mathrm{\Omega }}{\varphi }_o(\mathbf{x},t)dV,$(56)

where ${\varphi }_r(\mathbf{x},t)$ represents redistributed source terms, $S_{m,g}$, ${\mathbf{S}}_{u,g}$, $S_{h,g}$, $S_{Y_{g,k}}$, and the solid phase volume fraction ${\theta }_s$. A distance of distribution ($d_{max}$) is prescribed, within which the solid volume fraction and source terms will be redistributed from the position of the centroid of the local cell (${\mathbf{x}}_0$). A simple filtering kernel function $g^{\prime }(\mathbf{x},t)$ is employed,

(57) $g^{\prime }(\mathbf{x},t)=(1-\frac{|\mathbf{x}-{\mathbf{x}}_0|}{d_{max}}{)}^2,$(57)

The filtering kernel function $g(\mathbf{x},t)$ is then obtained from the normalization of $g^{\prime }(\mathbf{x},t)$,

(58) $g(\mathbf{x},t)=g^{\prime }(\mathbf{x},t)/{\int }_{\mathrm{\Omega }}{g}^{\prime }(\mathbf{x},t)dV,$(58)

satisfying,

(59) ${\int }_{\mathrm{\Omega }}g(\mathbf{x},t)dV=1.$(59)

Numerical scheme

The governing equations are numerically solved using an open-source CFD code, Open-FOAM v6 (Weller et al. 1998). The MP-PIC collision model was implemented for the discrete phase, and the DKM model was implemented for the coupled source terms. To facilitate source term redistribution, an efficient cell search algorithm was implemented, and a message-passing interface (MPI) strategy was developed for the parallel computation using DKM. The finite volume method (FVM) was used to solve the governing equations of the continuous phase, with spatial derivatives calculated using second-order “Gauss-Limited linear” schemes and temporal integration using first-order Euler schemes. Velocity-pressure coupling in the continuity and momentum equations was performed using the PIMPLE algorithm, which combines the advantages of the PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) algorithms. A semi-implicit algorithm was employed to handle source terms between the solid phase and gas phase.

Computational cases

The CFD model was applied to simulate the granular flow and combustion process in an industrial-scale biomass-fired CFB furnace, which has a thermal power output of 110 MW. A schematic illustration of the furnace is shown in Figure 2, with only half of the computational domain displayed due to the symmetry of the geometry.

Figure 2. Computational domain of the 110 MWth industrial-scale CFB furnace. Only half of the computational domain is shown due to the symmetry of the furnace geometry.

The furnace has a height of 30 m ($z-$direction), a width of 8.7 m ($x-$direction), and a depth of 5.4 m ($y-$direction). The primary air inlet, located at the bottom of the furnace, was rectangular with a cross-section of 8.7 m $\times$ 2.4 m, while the eight secondary air inlets were circular with a diameter of 0.4 m for the lower row and 0.475 m for the upper row at a height of $z=1.0$ m and $z=5.5$ m, respectively, as shown in Figure 2. The combustion chamber occupied the upper portion of the furnace, with a height of up to 28 m. The fluidized bed region containing sand particles was located in the lower part of the combustion chamber and had a height of 6.5 m. Above that was the region for volatile gas and char combustion, which had a height of 20 m and a cross-section of 5.4 m $\times$ 8.7 m. The cylindrical cyclones, which had a height of 9.25 m and a diameter of 4.0 m, were connected to the furnace via two 0.5 m diameter pipes for circulating solid particles. Finally, flue gas was directed to the top-box above the furnace and flowed out through an outlet region as indicated in Figure 2.

The proximate analysis data and physical properties of the initial biomass and sand particles are presented in Table 3. The biomass is a mixture of waste wood, wood chips, sawdust, and bark with a mass ratio of 6 to 3 to 2 to 1. At the start of the simulation, a total mass of 60,000 kg of sand particles was fed to the furnace, while the biomass was supplied from the secondary air inlets at the lower row nozzles at a mass injection rate of 12.7 kg/s. The properties of the sand and biomass particles, including size distribution described by the Rosin-Rammler distribution function, are also shown in Table 3. Although the sand particles are chemically inert, their temperature varies in space and time. The biomass particles are modeled as having constant density, but their size changes during the thermochemical conversion process, and they are removed from the computational domain once burned out.

Table 3. Initial biomass and sand particles used in the CFB furnace.

Item -	$d_i$ [mm]	$\overline{d_i}$ [kg/m${ }^3$]	${\rho }_i$ [J/kg$\cdot$K]	$C{p}_i$ [kg/s]	feed rate [t]	total mass [${ }^{\circ }$C]	$T_i$
Biomass	0.50 $\sim$ 5	3.15	750	1500	12.70	–	25
Sand	0.125 $\sim$ 0.5	0.22	2300	860	–	60	800
Ultimate [wt%]	C	H	O	N	Ash
Biomass	50.5	0.57	38.6	0.97	4.1

Note: $\overline{d_i}$ is the mean diameter.

Regarding boundary conditions, the computational domain’s outlet boundary (located on the left surface of the top-box) is set as a fixed-pressure boundary condition, with a zero-gradient condition assumed for other variables. The air inflow boundary is prescribed with Dirichlet boundary conditions, with the inlet flow velocity computed from the mass-flow rate condition. The wall boundary is non-slip and maintained at a constant temperature of 1173 K. A CFL number of 0.2 is used in the gas phase flow calculations.

The CFD mesh was generated using the ANSYS Workbench V17.2 package. The fluidization gas from the primary and secondary air inlets had a temperature of 200 ${ }^{\mathrm{o}}\mathrm{C}$, with a mass flow rate of 15.21 m${ }^3$/s and 22.8 m${ }^3$/s, respectively. To evaluate mesh independence, three sets of unstructured grids were used: a fine mesh with 604,634 cells, 988,551 sand parcels, and 100,000 biomass parcels; a medium-mesh with 512,286 cells, 743,568 sand parcels, and 80,000 biomass parcels; and a coarse mesh with 413,541 cells, 497,781 sand parcels, and 50,000 biomass parcels. The simulation results obtained using these mesh resolutions did not show any systematic difference nor any significant difference in cost.

Figure 3 compares the mean gas temperature along the centerline of the furnace, i.e., along $z-$direction at $x=0$ and $y=0$ (Figure 2). It is shown that the results from the three meshes are rather similar, with the results from the fine mesh and the medium mesh showing closer agreement. The results from the fine mesh are discussed in the following.

Figure 3. Mean gas temperature along the centerline of the furnace ($x=0$ and $y=0$).

The gas temperature in the combustion furnace is measured using thermocouples at 20 monitoring locations. The coordinates of the monitoring locations are presented in Table 4. The origin of the coordinate is at the center of the primary air inlet as shown in Figure 2.

Table 4. Locations of temperature measurement points in the furnace. $x$ and $y$ represent the two horizontal coordinates and $z$ represents the vertical coordinate, see Figure 2.

Location	P1	P2	P3	P4	P5	P6	P7	P8	P9	P10
x [m]	3	0	−3	−3	0	3	−4	4	−4	−4
y [m]	−0.5	−0.5	−0.5	0.5	0.5	0.5	0	1	1	0
z [m]	0.4	0.4	0.4	0.4	0.4	0.4	8.5	9.3	9.3	9.3
Location	P11	P12	P13	P14	P15	P16	P17	P18	P19	P20
x [m]	−4	4	4	−4	−4	−4	−4	4	4	4
y [m]	−1	−1	0	0	−1	0	1	1	0	−1
z [m]	9.3	9.3	9.3	13.5	24.5	24.5	24.5	24.5	24.5	24.5

Results and discussion

The MP-PIC and CGM methods were initially used to simulate the granular flow and combustion process. This method is also known as the particle centroid method (PCM) since source terms in a mesh cell are calculated from all parcels in the cell. However, the numerical simulation was found to be unstable, and no converged solution could be obtained. In contrast, the MP-PIC with DKM resulted in a stable numerical solution. Therefore, the following discussion will first address the local overloading problem in PCM and DKM.

Performance of PCM and DKM

Particle local overloading can occur when the solid volume fraction, ${\theta }_s=V_s/V_c$, exceeds the physical packing limit in a mesh cell. The physical packing limit refers to the maximum solid volume fraction that a local cell can accommodate. For spherical particles, this limit is approximately 0.62 (Sun and Xiao 2015; Yang et al. 2022), meaning that the lowest gas volume fraction value in a cell is 0.38. Figure 4 presents the distribution of gas volume fraction ${\alpha }_g$ with and without the use of DKM, as well as the local particle load $1/{\alpha }_{c/p}$ without DKM.

Figure 4. Distribution of gas volume fraction ${\alpha }_g$ and particle load $1/{\alpha }_{c/p}$ at the different heights of the furnace ($z=$ 2, 3 and 4 m), predicted using PCM and DKM.

(60) ${\alpha }_g=1-\frac{V_s}{V_c},{\alpha }_{c/p}=(\frac{V_c}{V_s}{)}^{1/3},$(60)

where $V_c$ and $V_s$ denote respectively the cell volume and the solid phase volume in the local cell. For the furnace height between $z$ = 2 $\sim$ 4 m, where particles are densely concentrated, local overloading is clearly evident. Without DKM (termed as PCM), the local particle load ($1/{\alpha }_{c/p}$) can reach as high as 0.94 at $z=2$ m. In this case, the gas volume fraction is as low as 0.17, which is below the physical limit of 0.38. In order to maintain numerical stability and avoid nonphysical solutions in the OpenFOAM solver, a numerical limiter is applied. Specifically, when ${\alpha }_g<0.38$, it is set to the value of 0.38. In the simulation with PCM, a large portion of the domain has this limiter applied. In the simulation with DKM, the region requiring the limiter is significantly smaller, especially at higher furnace height. The source terms used in the PCM method are rather “noisy,” meaning they are similar to the distribution of ${\alpha }_g$ and $1/{\alpha }_{c/p}$. This is likely the reason why the numerical simulation could not converge. Additionally, the numerical limiter applied to ${\alpha }_g$ causes it to be numerically increased, which can lead to errors in the numerical solution (Yang et al. 2022).

Figure 5 compares the gas temperature at different locations obtained from numerical simulations using the MP-PIC and CGM/DKM models. Table 4 lists the locations P1–P6 at the fuel-supplying region with $z=0.4$ m, where sand and biomass particles exchange heat, and the biomass particles initiate the thermochemical conversion. The gas temperature in this region is relatively low, ranging from 1100 K to 1150 K, and is uniform in the horizontal plane. P7 – P13 are at the furnace height $z=9.3$ m, above the secondary air inlets. The higher gas temperature in this region indicates that the thermochemical conversion process has progressed, and exothermic volatile reactions are taking place. However, the gas temperature at this furnace height is non-uniform in the horizontal plane. P7, P10, and P13 are on the symmetric plane with $y=0$ and near the side wall of the furnace $|x|=4$ m, where the gas temperature is relatively low and similar. P8, P9, P11, and P12 are off the symmetric plane with $|y|=1$ m and $|x|=4$ m, where the gas temperature is higher than those on the symmetric plane. P14 is at furnace height $z=13.5$ m, where the gas temperature is the highest. At higher furnace height ($z=24.5$ m), the gas temperature is slightly lower than that at $z=13.5$ m, but it is uniform in the horizontal plane, as shown at P15–P20. The numerical simulations using the MP-PIC and CGM/DKM models replicate the experimentally observed trend of gas temperature well.

Figure 5. Comparison of gas temperature between the numerical simulation using MP-PIC and DKM and experiment at different monitoring locations. The spatial coordinates of the 20 locations are given in Table 4.

Granular flow and characteristics of fluidization

To understand the temperature distribution discussed above, it is essential to consider the granular flow and fluidization process of sand and biomass particles. Figure 6 illustrates the instantaneous distribution of sand and biomass particles in the furnace and cyclones at an arbitrary time after the numerical simulation reached a statistically steady state. The figure also shows the distribution of biomass particles colored with their size and temperature, and the gas velocity streamlines colored with the gas temperature and gas flow velocity. Initially, the sand particles are deposited in the lower part of the furnace ($z=0\sim 6$ m), and the biomass particles are then injected. The fluidization air flow, supplied from the primary air inlet at a velocity of about 1 m/s and the secondary air inlet at a higher velocity ($>$ 10 m/s), fluidizes the sand and biomass particles until quasi-steady-state fluidization is reached after 15 s of physical time.

Figure 6. Spatial distribution of sand and biomass particles (left panel), biomass particles colored with particle size and temperature (second and third panels), gas flow streamlines colored with gas flow speed and gas temperature (fourth and fifth panels), at an instant of time during the stationary operation stage. The results are obtained using MP-PIC and DKM models.

The furnace can be divided into two distinct regions: the dense particle region located in the lower part of the furnace, within 8 m ($z<8$ m) above the bottom plane, and the dilute particle region located further up in the furnace and in the cyclones, i.e., $z>8$ m. Most of the particles are concentrated in the dense particle region, where the size of the biomass particles is relatively larger. The temperature of the biomass particles is relatively low near the inlet and increases quickly in the dense particle region. The gas flow in the dense particle region is rather complex due to the gas/solid interaction resulting from the granular flow. In the dilute particle region, the gas flow is accelerated when it enters the cyclones, forming a swirling flow motion. The swirling gas flow exhibits a “vortex breakdown” structure upon entering the top-box, where an inner recirculating zone can be observed. The gas flow exits the furnace at the outlet located on the left surface of the top-box. The particles in the cyclones are observed to be separated from the gas flow and returned to the furnace via the two connecting pipes. The gas temperature exhibits a locally cold region and a hot region in the lower part of the furnace, indicating the non-uniform nature of the dense particle region. Further up in the furnace, the gas temperature becomes more spatially uniform, which explains the larger spatial variation of gas temperature observed in the experiments at the furnace height of $z=9.3$ m (thermocouple locations P8–P13) and the more uniform temperatures observed at the furnace height of $z=24.5$ m (thermocouple locations P15 – P20).

Figure 7 illustrates the instantaneous distributions of biomass particles in the dense particle region at four different times. The first instant, $t=0$ s, corresponds to an arbitrary time after the furnace has reached a statistically stationary operation state. The figure also displays the cross-sectional averaged Sauter mean diameters of biomass particles, biomass particle velocity, and temperature at different furnace heights. Gas bubbles can be identified, for instance, in the bottom row of the figure and the location and size of bubbles evolve over time. The gas bubbles in the upper part of the dense particle region are large in size, and they periodically break up (e.g., at $t=0$ s) and form (e.g., at $t=3$ s). It can be seen that larger particles are located near the bottom of the furnace, due to gravity. These particles undergo drying, pyrolysis (devolatilization), char oxidation, and gasification while moving around in the bottom of the furnace. The particle temperature is higher near the primary air inlet than that further up in the furnace, due to the exothermic reactions of the particles. As the particles become smaller and lighter, they are blown upward in the furnace; hence, the mean diameter of the particles tends to decrease along the furnace height.

Figure 7. Spatial distribution of biomass particles in the dense particle region, colored with the particle diameter (upper row) $d_{32}$ and temperature (bottom row), and velocity vector of the biomass particles at $t=$0, 3 s, 6 s, and 9 s. $t=0$ is an arbitrary time after the flow and combustion process reach statistically steady states. The right panel shows the cross-section averaged Sauter mean diameter ($d_{32}$) of the particles, the velocity of the particles, and the temperature of the particles along the furnace height at $t=$0, 3 s, 6 s, and 9 s. The results are from numerical simulations using MP-PIC and DKM.

At the lower row of the secondary air inlet (close to the furnace height indicated by the tube connecting the cyclones), cold fresh biomass particles are injected, resulting in a relatively low mean temperature in this furnace height. Small particles tend to be found in the center of the furnace, where the particle velocity is low. The larger particles tend to move at a higher velocity and are found in the near-wall region around the gas bubbles, where the gas velocity is higher. The mean diameter of particles in the region from the particle inlet to the upper surface of the dense particle region may be increasing along the furnace height due to the larger particles flowing upward around the boundaries of the gas bubbles in the furnace. Further up in the dilute particle region, the biomass particles are smaller and hotter (due to loss of mass during thermochemical conversion). It is worth noting that the bubble formation and breakup are highly unsteady, leading to temporally evolving particle properties (diameter, velocity, and temperature). However, the overall trend of the particle characteristics discussed above is similar at different times, as shown in the right panel of Figure 7.

Thermochemical conversion process of biomass particles

Figure 8 displays cross-sectional, averaged gas-phase properties from numerical simulations, including gas temperature $T_g$, gas pressure drop $\mathrm{\Delta }{P}_g$, gas velocity $U_g$, gas volume fraction ${\alpha }_g$, and mass fractions of H${ }_2$ and CO${ }_2$. The pressure drop increases along the furnace height, varying rapidly in the dense particle region and reaching a plateau in the dilute particle region ($z>8$ m). The gas volume fraction also varies significantly in the dense particle region, becoming nearly 1 in the dilute particle region. Gas temperature increases slowly along the furnace height in the dense particle region, where fuel particles undergo drying and pyrolysis, releasing volatile gases such as CO, H${ }_2$, and CH${ }_4$, as well as CO${ }_2$ and H${ }_2$O, along with char. Combustion of CO, H${ }_2$, CH${ }_4$, and char in the dilute particle region is responsible for the continued increase in gas temperature along the furnace height until $z=18$ m and the rapid decrease in H${ }_2$ mass fraction. The rapid decrease in CO${ }_2$ at the surface of the dense particle region is likely due to rapid mixing with air that erupted from the gas bubbles. Further downstream, the gas temperature decreases slightly along the furnace height due to heat loss to the walls. This result is consistent with the temperature measurement shown in Figure 5, where the highest gas temperature is around P14 ($z=13.5$ m).

Figure 8. Cross section averaged gas properties at different heights of the furnace including gas temperature ($T_g$), pressure drop ($\mathrm{\Delta }\mathrm{P}$), gas velocity ($U_g$), gas volume fraction (${\alpha }_g$), and mass fractions of H${ }_2$ and CO${ }_2$.

In the dense particle region, the mean gas velocity ($U_g$) increases due to the supply of secondary air. The gas velocity decreases when it erupts from the dense particle region, reaching a relatively constant flow speed before accelerating upon entering the cyclones, as shown in Figure 6. The gas velocity profile provides insight into the particle velocity as shown in Figure 7. In most parts of the furnace, the gas velocity exceeds the particle velocity, indicating that the particles are dragged by the gas flow and accelerate/decelerate along with the gas flow. In the upper part of the furnace, the particle velocity is similar to the gas velocity, likely due to the smaller size of the particles, which tend to follow the gas flow.

NO^x formation process

Figure 9 shows the spatial distribution of mass fractions of nitrogen-containing species, including NO, NH${ }_3$, HCN, HNCO, CNO, and NCO, while Figure 10 shows the cross-section averaged nitrogen-containing species. The NO concentration exhibits a local peak at $z=2$ m, followed by a decrease to a low level at $z$ = 4–14 m, resulting from reactions with NH${ }_3$ and CH${ }_4$ through reactions R8–R10. This local peak is caused by biomass pyrolysis, which releases NO along with other volatile species, as shown in Equation (26)(26) $\begin{array}{c}Biomass\to x_1\mathrm{C}\mathrm{O}+x_2\mathrm{C}\mathrm{O}{}_2+x_3\mathrm{H}{}_2\mathrm{O}+x_4\mathrm{H}{}_2+x_5\mathrm{C}\mathrm{H}{}_4\\ +x_6\mathrm{H}\mathrm{N}\mathrm{C}\mathrm{O}+x_7\mathrm{H}\mathrm{C}\mathrm{N}+x_8\mathrm{N}\mathrm{H}{}_3+x_9\mathrm{N}\mathrm{O}+x_{10}\mathrm{C}{}_{(s)},\end{array}$(26) and Figure 9. Further up the furnace, the NO concentration gradually increases, reaching a level of about 80 ppm at $z$ = 24–30 m. The formation of NO primarily occurs in the dilute particle region ($z>8$ m), due to the oxidation of NH${ }_3$ (R7), CNO (R12), and NCO (R15).

Figure 9. Spatial distribution of mass fractions of nitrogen-containing species including NO, NH${ }_3$, HCN, HNCO, CNO, and NCO at different heights of the furnace.

Figure 10. Cross-section averaged mass fractions of nitrogen-containing species at different heights of the furnace including NO, NH${ }_3$, HCN, HNCO, CNO, and NCO.

A similar tendency of NO was observed in the study by Vainio et al. (2012), indicating reasonable predictions in the current simulation. The NH${ }_3$ concentration also exhibits a local peak at $z=$ 2 m due to the release of volatiles during pyrolysis. NH${ }_3$ is rapidly consumed along the furnace height by reactions with O${ }_2$ (forming NO via reaction R7) and with NO (consuming NO through reaction R8). Above $z=$ 10 m, NH${ }_3$ concentration is almost negligible. HNCO and HCN initially increase (due to biomass pyrolysis) and then decrease (due to volatile combustion, e.g., R11, R14) along the furnace height. CNO and NCO gradually increase along the furnace height and peak near $z$ = 24 m, largely due to volatile combustion (R11, R14).

As a summary, biomass pyrolysis mainly occurs in the dense region of the furnace ($z<8$m), resulting in extreme values of NO, NH${ }_3$, HCN, and HNCO, as seen in Equation 26(26) $\begin{array}{c}Biomass\to x_1\mathrm{C}\mathrm{O}+x_2\mathrm{C}\mathrm{O}{}_2+x_3\mathrm{H}{}_2\mathrm{O}+x_4\mathrm{H}{}_2+x_5\mathrm{C}\mathrm{H}{}_4\\ +x_6\mathrm{H}\mathrm{N}\mathrm{C}\mathrm{O}+x_7\mathrm{H}\mathrm{C}\mathrm{N}+x_8\mathrm{N}\mathrm{H}{}_3+x_9\mathrm{N}\mathrm{O}+x_{10}\mathrm{C}{}_{(s)},\end{array}$(26) . After pyrolysis, NO is reduced via reactions R8, R9, R10, R13, and R16. Following the injection of a large amount of air into the secondary inlets, NH${ }_3$ is converted to NO by R7 in large quantities, leading to an increase in NO and a decrease in NH${ }_3$. Above the secondary air inlets, HCN and HNCO are oxidized and converted to N${ }_2$ by reactions R11–R16. It should be noted that NO concentrations, as well as CNO and NCO, fluctuate significantly at $z$ = 24–30 m due to interaction with the cyclones, which induces the formation of unsteady rotational flow (swirling flow) structures, as shown in Figure 9.

Discussion

The information provided about the granular flow and thermochemical conversion in the CFB furnace can enhance our understanding of the flow and combustion process in the furnace. The chemical kinetic model can be used to explain the NO formation process, while the bubble formation, breakup process, and division of dense and dilute particle regions can explain the gas temperature field and the interaction of the particles with secondary air. By analyzing CFD simulations, the operation of the furnace can be improved by optimizing the primary air and secondary air supply for different biomass fuels to achieve better fluidization in the dense particle region and combustion in the furnace. However, it is important to note that the present CFD results need to be thoroughly validated under industrial CFB boiler operating conditions. Due to the lack of experimental data, the current CFD results are only validated against mean gas temperature experimental data at a few sampling locations. Therefore, more experimental data on the gaseous species and particles in the furnace are desirable for further validation of the model.

Conclusions

A recently developed 3-D computational fluid dynamics (CFD) model has been employed for numerical simulations of biomass combustion in an industrial-scale circulating fluidized bed (CFB) furnace. The CFD model is based on the multi-phase particle-in-cell (MP-PIC) collision model, coarse grain method (CGM), and a distribution kernel method (DKM) that aims to resolve the local particle overloading issue typically found in the conventional particle centroid method (PCM). The hydrodynamic and combustion properties of the solid and gas phases are analyzed to provide insight into the physical and chemical processes in the furnace. The main conclusions drawn are as follows:

• The CFD model can well overcome the dense particle local overloading problem in the industrial CFB furnace. Without the use of DKM, the CFD simulation could not achieve any stable solution; with DKM, the model can simulate the dynamics of the granular flow and thermochemical conversion of the particles. The predicted temperature field agrees well with the thermocouple experiments in various locations in the furnace.

• The CFD results show that the CFB furnace can be divided into different regions based on the characteristics of the granular flow. In the lower part of the furnace, there is a dense particle region where most particles are located. In this region, gas bubbles are formed and evolve in space and time. The bubbles break up near the upper boundary of the dense region. Above the dense particle region is the dilute particle region, where the particles are smaller and lighter and tend to follow the gas flow. Further downstream, the tiny particles are separated in the cyclones and returned to the furnace through connecting pipes.

• When biomass particles are supplied to the furnace in the dense particle region, the larger particles tend to follow the high-speed gas flow in the boundary of the gas bubbles or fall to the bottom of the furnace due to gravity. The Sauter mean diameter of the particles is relatively low in the fuel injection region due to the fall of larger particles toward the bottom of the furnace.

• Drying and pyrolysis of the biomass particles occur mainly in the dense particle region. Oxidation of volatile gas and char particles continues in the dilute particle region. This explains why the highest temperature in the furnace is in the mid-height region, where most of the volatile gas is combusted. Further downstream, the gas temperature becomes more uniform in space, and the gas temperature is slightly lower than that in the mid-height region.

• Biomass pyrolysis in the dense particle region contributes to the release of NH${ }_3$, HCH, and HNCO. The combustion of volatile gas further up in the dilute region contributes to converting the nitrogen-containing species to CHO, NCO, and NO.

Acknowledgements

The authors would like to thank Kraftringen AB for supporting the detailed information on the CFB furnace. Miao Yang is sponsored by the China Scholarship Council (201808410350). The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at PDC (Beskow).

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Additional information

Funding

This work was supported by the Swedish Energy Agency (STEM) through KC-CECOST, project Nr [22538-4], and the Knut & Alice Wallenberg foundation (KAW COCALD project)

Notes on contributors

Miao Yang

Miao Yang Methodology, Investigation, Visualization, Writing-Original draft preparation.

Shenghui Zhong

shenghui Zhong Methodology, Investigation, Visualization

Shijie Xu

Shijie Xu & Leilei Xu Methodology, Visualization, Writing-Original draft preparation.

Leilei Xu

Peter Ottosson Investigation, Writing-Original draft preparation.

Peter Ottosson

Hesammedin Fatehi Methodology, Investigation, Writing-Original draft preparation.

Hesameddin Fatehi

Xue-Song Bai Funding acquisition, Writing-Reviewing and Editing, Resource, Supervision.