A fast-solving particle model for thermochemical conversion of biomass

A fast-solving particle model for thermochemical conversion of biomass

a,∗ b

Tian Li , Henrik Thunman , Henrik Strömc

aDepartment of Energy and Process Engineering, Faculty of Engineering, NTNU — Norwegian University of Science and Technology, Trondheim, Norway

bDivision of Energy Technology, Department of Space, Earth and Environment, Chalmers University of Technology, G¨oteborg, Sweden

cDivision of Fluid Dynamics, Department of Mechanics and Maritime Sciences, Chalmers University of Technology, G¨oteborg, Sweden

Abstract

Computational fluid dynamics (CFD) simulations of large-scale furnaces or reactors for thermal conversion of solid fuels remains challenging partially due to the high computational cost related to the particle sub-models. Ow- ing to the thermally thick nature, it is particularly expensive to simulate the conversion of large fuel particles such as biomass particles. To address this issue, a fast-solving particle model was developed in this work with special attention to the computational efficiency. The model spatially discretizes a fuel particle in one homogenized dimension. The conversion process of the fuel particle is treated as a reactive variable-volume one-dimensional tran- sient heat conduction problem. The model also utilizes several features that are typically found in sharp interphase-based models to reduce the compu- tational cost. Validation of the model was carried out by comparing with experimental results under both pyrolysis and combustion conditions. The

∗Corresponding author

Email address: [email protected] (Tian Li)

Author’s Post-print version. Article published in Combustion and Flame, Volume 213, March 2020, Pages 117-131.

DOI: 10.1016/j.combustflame.2019.11.018

Published version: https://www.sciencedirect.com/science/article/pii/S0010218019305231

accuracy and computational efficiency of the model was thoroughly exam- ined by varying the degrees of temporal and spatial discretization. It was found that the model well predicted pyrolysis and combustion of a single biomass particle within a broad range of temporal and spatial discretization.

The time used to simulate the conversion of a biomass particle using the developed model can be more than one order of magnitude smaller than the conversion process itself. It was also revealed that a well-predicted conduc- tive heat transfer inside the particle is essential for a precise simulation of the drying and devolatilization process. The char conversion process, however, is less sensitive to the external heat transfer as it is mainly controlled by the mass diffusion process. Further studies showed that a time step of 1×10⁻³ s and a spatial discretization of 20 cells were sufficient for simulating the con- version of typical fuel particles in grate-fired and fluidized-bed furnaces. We also demonstrated that when the particle model was implemented in a CFD solver, only 2.2% of computational overhead was introduced by the model.

As the model can efficiently employ fixed time stepping, optimal load bal- ancing during parallel computing of many simultaneous conversion processes becomes trivial. This performance opens up new possibilities for treating fuel polydispersity in Eulerian CFD simulations of biomass conversion.

Keywords: biomass, combustion, mathematical modeling, CFD

Nomenclature Symbols

A pre-exponential factor (1/s)

AC ash content (-)

a matrix coefficient (J/K s) b matrix coefficient (J/s)

C molar concentration (mol/m³) c_p heat capacity (J/kg K)

d initial diameter (m) dt time step (s)

∆H heat of reaction (J/kg) E activation energy (kJ/mol)

h convective heat transfer coefficient (J/m² K s) k_{ef f} effective rate constant (m/s)

l initial length of a finite cylinder (m) l₁,l₂,l₃ initial dimensions of a parallelepiped (m) M molar mass (kg/mol)

m mass of a given species in a cell (kg) mp particle mass (kg)

M C moisture content (-) N p number of cells (-)

r homogenized radial position (m) rb particle radius (m)

RCT relative calculation time (-) RE relative error (%)

S mass source during a given conversion stage (kg/s) S_g specific gravity (-)

S_p implicit coefficient for the source term of energy equation (J/K s) S_u explicit contribution to the source term of energy equation (J/s) S_energy source term of energy equation (J/s)

SA surface area (m²)

SF shrinkage factor during a given conversion stage (-) t time (s)

T_g gas temperature (K) T_s solid temperature (K)

T_w reactor wall temperature (K)

tHT timescale of conductive heat transfer (s) t_total,cal total time used for executing the code (s) t_total,con total conversion time of a particle (s)

V volume of a given species in a cell (m³) V_c volume of a cell (m³)

Vp particle volume (m³) w mass fraction (-)

xb location of the outer cell boundary (m) xc location of the cell mass center (m) Greek

α fraction of oxygen consumed by volatile gases (-) β mass transfer coefficient (m/s)

Γ particle shape factor (m²) emissivity (-)

θ radiation temperature (K) κ thermal conductivity (J/m K s) ρ density (kg/m³)

σ Stefan Boltzmann constant (kg/s³ K⁴) φ volume fraction (-)

χ coefficient of the explicit discretization of the energy equation (J/K s) Ω stoichiometric coefficient (-)

Subscript

∞ far field ash ash

b cell outer boundary i.e. west boundary BL baseline

C carbon atom char char

charc char conversion process charf char front

db dry based

dev devolatilization process dry drying process

dry wood dry wood E east cell gas gas

H hydrogen atom H2O water molecule O oxygen atom

O2 oxygen molecule P cell center sur particle surface tar tar

vol volatiles W west cell wb wet based wet wood wet wood Superscript

∗ current time step 0 previous time step

Abbreviations

CFD computational fluid dynamics DEM discrete element method IBM interface-based model MBM mesh-based model

ODE ordinary differential equation PDE partial differential equation

SFOR single first-order

SIMPLE semi-implicit method for pressure linked equations

1. Introduction

Biomass is a promising renewable source of heat, power, liquid fuels, and chemicals. It has gained popularity due to the increased awareness of sustain- ability, especially for the lignocellulosic biomass, which has been widely used in the Nordic countries to replace fossil fuels [1]. Due to the fibrous nature of lignocellulosic biomass, it is energy-intensive to reduce the size of a biomass particle. Consequently, biomass fuel particles used in industrial-sized fur- naces are usually much bigger than coal particles. Moreover, the proportion and rate of release of volatiles differ significantly between coal and biomass.

This situation makes the often-used thermally thin assumption (i.e. negligi- ble internal temperature gradients) in the modelling of coal conversion less suitable when simulating conversion of biomass. In addition, a recent study has suggested that temperature gradients inside a biomass particle need to be considered even under suspension-firing conditions [2], where the major- ity of fuel particles are less than 500µm. The internal temperature field not only influences the rate of conversion, but also to a great extent the biomass decomposition products.

To more accurately simulate the conversion of biomass particles, non- isothermal models accounting for the thermal thickness have been devel- oped, ranging from simplified one-dimensional modelling approaches to com- prehensive multi-dimensional models. A comprehensive overview of these various modelling approaches has been presented recently [3], in which multi-

dimensional models are recommended only when there is a need to account for the anisotropy of the fuel particle or when the boundary conditions of the fuel particle are varying significantly in space [3]. Both these scenarios are, however, rare in large-scale biomass-fueled reactors. Despite the anisotropy of the raw lignocellulosic biomass, the fuel particles used in industrial-size fur- naces can generally be characterized as isotropic due to pretreatments such as size reduction and densification. Even where this is not the case, such as for wood chip fuels, the accuracy obtained from one-dimensional particle mod- els is usually satisfactory without accounting for particle anisotropy [3]. The assumption of homogeneous boundary conditions for the fuel particle is also typically acceptable [4] since the size of a fuel particle is much smaller than that of the fuel bed. Moreover, one-dimensional particle models are also typ- ically linked to the surface area per unit volume of the particle phase, rather than to the surface area of an individual particle, when used in Eulerian fixed-bed reactor models [5]. Consequently, one-dimensional models with their associated reduction in computational cost are in general well suited to simulations of fuel particle conversion in large-scale biomass-fueled reactors.

This fact is even more pronounced when the particle conversion model is to be coupled to a computational fluid dynamics (CFD) solver, as but one out of many sub-models in what necessarily becomes a rather computationally expensive configuration.

Generally, one-dimensional models can be categorized into two concep- tual groups: interface-based models (IBMs) [5, 6, 7, 8, 9] and mesh-based models (MBMs) [10, 11, 12, 13, 14, 15, 16, 17]. In the IBM, the fuel particle is divided into three (moist wood, dry wood, and char) or four (moist wood,

dry wood, char, and ash) distinct layers. The chemical reactions and phase changes are then described as occurring only in narrow regions between the different layers (i.e. at sharp interfaces). With such a spatial discretization, the governing equations for each layer become ordinary differential equations (ODEs) coupled via boundary conditions at the interfaces. Thus, very few equations and variables need to be solved and stored respectively, and the relatively low memory requirements were also a main driving force for the original development of these types of methods. In the MBM, a particle is discretized along the radius into an ensemble of mesh points over which the conservation equations are spatially discretized. Unlike the IBMs, each discretized unit/cell may contain several different species in the MBM. All governing partial differential equations (PDEs) are converted into an alge- braic equation system by discretization in space and time. Depending on the gradients of temperature and species mass fractions, various numbers of grid points are needed for a reasonable representation of the particle con- version process. In some combustion cases, the grid points for a single fuel particle can be as many as 1200 [10], which potentially provides informa- tion in much greater spatial detail compared to the more coarsely discretized interface-models.

The IBMs are intuitively believed to have a higher numerical efficiency than the MBMs due to the reduced number of equations. Comparative stud- ies between IBMs and MBMs are however rarely reported. It was briefly stated in a previous study [7] that the developed IBM could achieve a similar level of accuracy as a more comprehensive MBM [15] but at lower compu- tational cost. However, it is worth noting that IBMs have some inherent

characteristics, which may limit their computational performance. For in- stance, two different sets of energy equations need to be solved sequentially at the boundaries and layers, since the solutions of the algebraic equations at the boundaries provide boundary temperature data that need to be known to advance the layer temperatures in time. Typically, a higher-order explicit method [7] or an iterative implicit method [6] is used to solve the ordinary differential equation for the layer temperatures. However, the accuracy of both methods, as well as the robustness of the explicit method, is in prac- tice limited by the smallest thermal mass of any present layer, indicating that very small time steps are needed in the beginning and towards the end of each conversion process (drying, devolatilization and combustion). Some sort of adaptive time stepping is therefore typically required, as well as a numerical limiter on the smallest layer or particle size allowed [5, 18]. On the other hand, it is common that linearly implicit extrapolation methods are utilized in the MBMs, which can be very efficient particularly when the size of the matrix is small. One reason for the high computational cost of the traditional MBMs is the consideration of intraparticle flow, which is usually ignored in the IBMs. In traditional MBMs, five sets of transport equations need to be solved of which four are related to the gas phase (continuity, mo- mentum, energy, and species). The convective timescale of the off-gas may be smaller, if not comparable, to the timescale of conductive heat transfer in the solid phase, thus, resulting a lower maximum allowed time step. A simplified MBM without solving intraparticle flow may therefore achieve a good balance between accuracy and computational efficiency.

Even though the different types of particle models discussed here form the

basis for comprehensive CFD simulations of reactors for biomass conversion, the implications of the particle model formulation for the computational ef- ficiency of the overall simulation framework are in general not discussed in depth. In particular, biomass combustion systems are typically characterized by a strong flow-chemistry coupling and chemical time-scales that are sev- eral orders of magnitude smaller than those relevant for the global conversion process [19]. In such situations, it is customary to invoke operator-splitting and compute the source terms in the energy and species transport equations for the duration of a fluid-flow time step [20, 21]. This repeated stopping and re-initialization of the biomass conversion sub-model has profound con- sequences for the global efficiency of the numerical method. For example, higher-order time integration methods require storing of information from multiple previous time levels and are associated with a larger computational effort concentrated to the beginning of each CFD time interval [21]. Such methods, which may exhibit superior performance to lower-order alterna- tives in stand-alone codes, can therefore, perhaps somewhat unexpectedly, deteriorate the overall performance of the final CFD implementation. In the present work, we therefore choose to focus on the balance between computa- tional efficiency, numerical accuracy and predictive capability achievable for biomass particle conversion models intended for use in biomass conversion reactor assessment.

Our long-term ambition with the current work goes further than a mere reduction in the computational cost of particle sub-models in biomass con- version simulations. We target our particle model development at a CFD simulation scenario where the fluid flow time step is in the order of 1×10⁻³

s or less. If a particle model can be derived such that the computational overhead in a typical fixed-bed conversion simulation is lower than e.g. 10%, then new opportunities will present themselves. More specifically, if the cost of executing one instance of the particle model is small enough, more in- stances can be afforded with little penalty to the overall computational cost.

These additional instances could enable a more elaborate treatment of par- ticle polydispersity or finer spatial resolution in Eulerian simulations. The success of such a computational setup would also in practice depend on the possibilities to achieve efficient parallelization of many simultaneous particle conversion processes. We shall return to these possibilities in the discussion of the results.

The remainder of this paper is organized as follows. In Section 2, a sim- plified MBM is presented. Besides neglecting the demanding solution of the intraparticle gas flow, the proposed model combines the concept of sharp reaction fronts from the IBM and the advantage of the traditional MBM by using an efficient matrix solver (LU factorization). The modeling details of phase change, reaction, thermophysical parameters, and solution strategy are also described. Validation of the model is discussed in Section 3 by com- paring with the experimental data and results from an IBM. In Section 4, comprehensive analysis of the model is carried out to systematically scru- tinize effects of model formulation on the predictability and computational efficiency in various conditions. A specific aim is to establish a solid theoret- ical foundation for interpretations of the observed results. Moreover, general guidelines for CFD simulations using one-dimensional thermally thick parti- cle sub-models are given, and new possibilities offered by the herein presented

particle model are discussed.

2. Theory

2.1. Model framework

As the present model is developed based on an, in our opinion, well- established baseline model [5, 6], implementation details regarding the physic- ochemical modeling of the conversion processes are not discussed in detail here. We have strived to achieve as similar implementations as the overar- ching model frameworks allow. The same argument can be made for the choice of material and transport properties – it is not the main purpose here to assess the relative merits of different correlations and literature sources on these values, but to illustrate, analyze and contrast the model behaviors for an acceptable set of property values typical of thermochemical conversion of woody biomass. The main differences between the present MBM, com- prehensive MBMs, and IBMs are summarized in Table 1 (the characteristic features of comprehensive MBMs and IBMs are exemplified using previous works [6, 15]). It should be noted that, although there is no direct resolution of gas-phase species concentration fields in the present MBM or in classical IBMs, effects of chemical reactions on outflowing gases may still be included in these modelling frameworks by integration of the reaction rate expressions, given that the residence time through the exterior layers and an approximate temperature profile in these layers are always known.

2.2. Particle heat balance

In the discretized particle model proposed in the current work, the ther- mochemical conversion of a biomass particle is essentially treated as a re-

Table 1: Main differences between the present MBM and typical comprehensive MBMs and IBMs.

Present MBM Comprehensive MBM IBM Hydrodynamic con-

sideration of intra- particle flow

Not considered Darcy’s law or full Navier-Stokes equa- tions

Not considered

Direct resolution of gas-phase species transport inside par- ticle

Not considered PDEs for each gas- phase species includ- ing temporal and spa- tial gradients, convec- tion and source terms

Not considered

Chemical reaction of gas-phase species in- side particle

Not considered Arrhenius expressions Not considered

Energy conservation PDEs for the solid phase

Combined PDEs for both gas and solid phases

ODEs for the tem- peratures of the solid phase layers and alge- braic equations at the boundaries between layers

Solver Matrix solver (Gaus-

sian elimination by LU factorization)

Matrix solver (Gaus- sian elimination, e.g.

tridiagonal matrix algorithm; possibly with auxiliary algo- rithms for handling pressure-velocity cou- pling, e.g. SIMPLE)

ODE solver

active variable-volume one-dimensional transient heat conduction problem.

The conceptual basis for the model is a previously developed IBM [5, 6], which has here been fundamentally reformulated in adaptation to the MBM framework. As illustrated in Figure 1, the biomass particle is divided into a number of segments. In a conventional IBM, these segments would be few and correspond to moist wood, dry wood, char and ash. Here, they merely represent a spatial discretization of the particle and they may, at any given point in time, correspond to either of the four previously listed constituents, or a mixture thereof. The basis for the MBM formulation is the treatment of the local particle properties as mass-weighted averages of the four basic constituents, in effect turning the MBM into a “one-solid” model in the vein of “one-fluid” models for segregated multiphase flows.

Figure 1: Schematics of one-dimensional MBM.

The governing equation to be solved is thus:

m_pc_p∂T_s

∂t = V_p Γ

∂

∂r

Γκ∂T_s

∂r

+S_energy (1)

where m_p is the particle mass (kg), c_p is the particle heat capacity (J/kg K), T_s is the particle temperature (K), t is the time (s), V_p is the particle volume (m³), Γ is the particle shape factor (m²),r is the homogenized radial position (m), κ is the particle thermal conductivity (J/m K s) andS_energy is the source term (J/s) due to the conversion processes (drying, devolatiliza- tion and char combustion). Expressions for the particle shape factor are given in Table 2. Radiative heat transfer inside the particle may be included via a semi-empirical contribution to the particle thermal conductivity [15].

However, this contribution is neglected in the current work, as the process is mostly heat-transfer controlled at lower temperatures (during drying and devolatilization) and mass-transfer controlled during char combustion (when the temperatures are high).

The particle properties in Equation (1) are mass-weighted averages based on the local presence of moist wood, dry wood, char and ash. Equation (1) is discretized using a fully implicit finite-volume method, with a first- order backward differencing scheme for the transient term and a second- order central differencing scheme for the conduction term. The source term linearization (S_energy =S_u+S_pT_s) uses the temperature at the new time level to harmonize with the implicit treatment, and the discretized equation [22]

for the temperature becomes:

T_s,P = χ_WT_s,W +χ_ET_s,E+χ⁰_PT_s,P⁰ +S_u

χ⁰_P +χ_W +χ_E −S_p (2)

Table 2: Particle shape factors.

Γ Infinite plate¹ 1 Infinite cylinder¹ 2πr

Sphere 4πr²

Finite cylinder² 2π(3r²+r(l−d))

Parallelepiped³ 24r²+8r((l₂−l₁)+(l₃−l₁))+2(l₂−l₁)(l₃−l₁)

1 For an infinite cylinder, Γ is related to a unit length of the cylinder, and for infinite plates to a unit area.

2 Initial length land initial diameter d.

3 With initial dimensionsl₁×l₂×l₃, wherel₁ represents the shortest end.

where the coefficientsχ(J/K s) and temperaturesT_shave subscripts that re- fer to the east (E), west (W) and cell center (P) positions, and a superscript of 0 refers to the previous time level. The coefficients and all ancillary expres- sions used in the discretized model are presented in the Appendix A to this paper. The boundary conditions employed are zero gradient (symmetry) at the particle center and a known temperature at the outer particle boundary.

The particle surface temperature is obtained from a heat balance where the conductive heat flux into the particle is equated to the external (convective and radiative) heat flux, under the assumption of no heat accumulation at the surface. A uniform temperature field is provided as the initial condition.

The spatial discretization is performed so that the biomass particle is di- vided into a number of cells of equal initial mass in its radial direction. Once

the coefficients and source terms for all cells have been calculated, a system of algebraic equations can be assembled. By neglecting the temperature- induced variation of the coefficients over the time step, this equation system can thereafter be solved directly by LU factorization with a matrix solver (Meschach [23]). The combined method is stable, owing to its implicit char- acter, and the employment of a first-order temporal scheme minimizes over- head losses due to re-initialization after each fluid flow time step [21]. The limitations in accuracy of the method brought about by the implicit treat- ment [24] and the further approximation of slowly varying coefficients can be offset by employing a small enough time step. Here, we take advantage of the fact that a CFD simulation of biomass conversion anyhow must include outer loops where the fluid flow solution is advanced in time with the newly updated source terms from the particle conversion sub-model. This iteration between the particle and the fluid level sets an upper limit on the time step for the sub-model, implying that higher-order methods that would allow time steps longer than the time scales on which the fluid flow must be updated are anyway unnecessary.

2.3. Drying, devolatilization and char conversion

Two types of drying models are common in the literature: thermal drying models and kinetic drying models. In this study, a single first-order (SFOR) kinetic model with an Arrhenius type rate expression (A_dry = 5.13×10¹⁰ 1/s and E_dry = 88 kJ/mol [10, 15]) is chosen in order to guarantee that the tem- perature field is everywhere differentiable and to improve numerical stability.

Some previous studies have suggested that two types of water are present in woody biomass (free and bound water), the evaporation processes of which

should be modeled separately [15] since they may behave differently. How- ever, the determination of free water and bound water contents is nontrivial, requiring complicated analytical methods (e.g. nuclear magnetic resonance scanning). Considering the low availability of the required data and addi- tional uncertainties it may bring, we do not distinguish between free water and bound water in the current model. Previous experimental work has also indicated that the effect of more tightly bound water on the overall drying process is insignificant under conditions relevant to the current applications [25]. Furthermore, the condensation of water vapor inside the particle is also neglected in the current model, but there is nothing in the model formulation that precludes an extension to include this effect in the future.

Devolatilization of biomass is a very complex process involving decom- position of carbonaceous materials into permanent gases, tar, and char. An accurate prediction of this process is usually considered very difficult. The promising multispecies multistep models [26, 27] and structural models (Bio- FLASHCHAIN [28], Bio-CPD [29], and FG-Bio [30]) are simply too expensive to be used in the CFD simulation of large-scale reactors, despite their poten- tials to predict product distributions over a wide range of temperature. On the other hand, simple and fast-solving SFOR devolatilization models require information of char yield in advance, making their application less reliable.

Hence, a competing devolatilization scheme is selected in this study to bal- ance the computational cost, accuracy, and applicability. The formation of light gases, tar and char are described by three competing SFOR reactions (A_dev,gas = 1.11×10¹¹ 1/s, E_dev,gas = 177 kJ/mol,A_dev,tar = 9.28×10⁹ 1/s, E_dev,tar = 149 kJ/mol, A_dev,char = 3.05×10⁷ 1/s, E_dev,char = 125 kJ/mol),

respectively, which were derived from conditions close to this study [31, 32].

In addition, the selected kinetic parameters have been successfully applied to simulate conversion of large biomass particles [2, 15]. In the current work, this set of the devolatilization kinetics is used in all simulations to show col- lective behaviors, as it is typically difficult to find kinetics matching both wood species and devolatilization conditions.

This study aims to develop a model for combustion of biomass in fixed or fluidized beds where a steep temperature gradient exists inside the par- ticle. Therefore, char is assumed to be converted in an infinitely narrow region similar to the method used in a previously developed IBM [5]. The adoption of this assumption in the MBM avoids the requirement of knowing concentrations of gas-phase reactants at every cell inside the particle, thus, greatly reducing the computational cost. The mass source for char during its conversion (S_charc,char) is given by the following equation:

S_charc,char =−(1−α)CO2,∞Ωk_{ef f}SA_charM_c (3) where α is the fraction of oxygen that is consumed by the volatile gases (-), C_O2,∞ represents the oxygen concentration in the far field (mol/m³), Ω is the stoichiometric factor of the char combustion (-),k_{ef f} is the effective rate constant estimated by taking both chemical kinetics and mass diffusion into consideration [5, 19] (m/s), SA_char is the surface area of the char front (m²), and M C is the molar mass of carbon (kg/mol). As shown in a previous work [33], conversion of relatively large char particles at high temperature is usually dominated by the oxidation reaction. Therefore, the gasification re- actions are not included in this study. Various reaction schemes for modeling

of char conversion have been proposed previously. With the goal of achieving an acceptable balance between accuracy and computational cost, we choose a global scheme [34]. This is an efficient model that has been widely used to simulate combustion of biomass [5, 19, 35, 36, 37, 38, 39]. Moreover, it has been well documented that the external mass transfer is the rate-limiting mechanism for char oxidation in biomass-fueled grate-fired and fluidized-bed furnaces [33, 40, 41], so different kinetic parameters may only have limited effect on the overall conversion process of the char. Since the current single particle model is intended to be coupled with a CFD solver where gas phase reactions are handled, homogeneous gas phase reactions inside the particle are not explicitly calculated. The consumption of oxygen by volatile gases is instead estimated by a factor α:

α= min





−Sdev,dry wood(1−AC_{dry wood})

wC

MC +_4M^wH

H −_2M^wO

O

SA_surCO2,∞β

−S_dev,char(1−AC_char)_M¹

C

SA_surC_O2,∞β ,1

!

(4)

whereSdev,dry wood and S_dev,char are mass source terms for dry wood and char during the devolatilization process respectively (kg/s),w_C,w_H, w_O are mass fractions of carbon, hydrogen, and oxygen in the dry wood, respectively (-), AC_{dry wood} and AC_char are ash contents in dry wood and char, respectively (-), M_H and M_O are molar mass of hydrogen and oxygen (kg/mol), SA_sur is the particle surface area (m²), and β is the mass transfer coefficient (m/s).

To simplify the estimation, char is considered to contain only carbon and volatile gases are assumed to be converted fully into CO₂ and H₂O.

2.4. Additional model information

Shrinkage is a non-negligible phenomenon accompanying the entire con- version process of a biomass particle including drying, devolatilization, and char conversion. It is influenced by several factors, such as type of biomass, water content, biomass density, size of the particle, and external heating con- ditions [42]. Given that shrinkage depends on many parameters, it should ideally be measured experimentally at the conditions of interest [43]. Here, the volumetric shrinkage from dry wood to char (SF_dev) and from char to ash (SF_charc) is estimated by constants (28% and 95%, respectively), as de- scribed in a previous study [5]. Shrinkage during drying is calculated based on wood species and moisture content as explained in the literature [42]. It is important that the thermal properties are chosen in a consistent way with regard to the changes the particle undergoes during conversion (mass loss and shrinkage) [44]. This is, however, not a straightforward task particularly for the determination of the thermal conductivity due to the variations in wood species and moisture content. A general equation describing thermal conductivity perpendicular (tangential and radial directions) to the grain [43]

was used in this study, whereas the thermal conductivity parallel to the grain was estimated as 1.8 times the value in the perpendicular direction [43]. The overall thermal conductivity was calculated by averaging in all three direc- tions [44]. Regarding the thermal conductivity of char, a wide range of values have been reported with magnitudes (in J/m K s) ranging from 10⁻² to 10⁰ [2, 3, 15, 19, 43, 44, 45]. Thus, a median value of 0.1 was chosen in this study. Details about the thermal conductivity and other model parameters are listed in Table 3.

Table 3: Model parameters.

Parameter Expression Ref.

Specific heat capacity c_p (J/kg K)

wet wood:

(4.206T_s−37.7) (1−M C_wb) + 4309M C_wb +(23.55T_s−1320M C_db−6191)M C_db

[45]

dry wood: 4.206T_s−37.7 [45]

char: −119×10⁻¹²T_s⁴+ 1010×10⁻⁹T_s³

−3160×10⁻⁶T_s²+ 4410×10⁻³T_s−334.0

[45]

ash: 754 + 0.586 (T_s−273.15) [19]

Thermal conductivity κ (J/m K s)

wet wood:

3.8 (S_g(0.1941 + 0.4064M C_db) + 0.01864)/3 Sg is the specific gravity

[46]

dry wood: 3.8 (0.1941S_g+ 0.01864)/3 [46]

char: 0.1

ash: 754 + 0.586 (T_s−273.15) [19]

Initial densityρ (kg/m³)

beech: 718 (M C_wb = 20%) [42]

poplar: 844 (M C_wb = 40%) [42]

boxwood: 830 (M C_wb = 0%)¹ [47]

SFdry (-)

beech: 14% (M C_wb= 20%) [42]

poplar: 13% (M C_wb = 40%) [42]

boxwood: 0% (M C_wb = 0%)¹ Particle emissivity

(-)

0.85 [15]

1 M C_wb is assumed to be zero, as it was not directly stated in the original

All the simulations were calculated sequentially using one of the compu- tational cores from an Intel^R Xeon^R X5650 CPU (2.66 GHz, Intel^R Turbo Boost Technology disabled, Intel^R Hyper-Threading Technology disabled).

The overall solving procedures for the IBM and the MBM are outlined in Figure 2.

Initialize solution

Calculate source terms

Update particle surface temperature

Determine ma- trix coefficients

Solve matrix system to obtain particle internal tempera- ture field att+∆t Advance

solution tot+∆t

Update mass in every cell

Initialize solution

Update interface temperatures

Calculate derivative of the layer temperature with respect to time

Test update of layer temperature tot+∆t

Update interface temperatures

Advance solution tot+∆t

Convergence? Yes No

(a) (b)

Figure 2: Solving procedures for the MBM proposed in the current work (a) and the IBM on which it is based (b).

3. Model validation

The developed MBM has been validated against experimental data from three different studies [14, 15, 47]. The first experiment is pyrolysis of a

spherical beech particle (d = 0.02 m, M C_wb = 20%) in a nitrogen atmo- sphere at a temperature of 1098 K. The ambient gas temperature (T_g) and reactor wall temperature (T_w) in the simulation was considered to be the same as they were not distinguished in the original study [14]. The second experiment is combustion of a near-spherical poplar particle (d= 0.0095 m, M C_wb= 40%) in air. The reactor wall temperature and ambient gas temper- ature were configured as 1050 K and 1273 K, respectively in the simulation according to the experiment [15]. Due to insufficient information provided in the experimental studies, the convective heat transfer coefficient was esti- mated from a correlation [48]. In both simulations, the particle is discretized into 50 cells. A relatively small time step of 1.0×10⁻⁴ s was chosen to ensure stability and accuracy of the solution.

In addition to the experimental data, simulation results from other well- validated models were also added. These include three comprehensive MBMs (MBM 1 [14], MBM 2 [17], and MBM 3 [15]) and two IBMs (IBM 1 [5, 6]

and IBM 2 [18]) developed focusing on the numerical efficiency. We aim here to illustrate the general performance of the single-particle model with the assumption of constant boundary conditions. It should also be noted that the same property and kinetic parameters were used in IBM 1 and the proposed MBM, so that a direct comparison between the two different modeling philosophies can be made.

The evolutions of particle surface temperature, particle central tempera- ture, and mass loss for pyrolysis of a beech particle are shown in Figure 3.

As shown in Figure 3, a good agreement between the model predictions and the experimental data can be found for the evolutions of both temperatures

and mass loss ratios. In general, with thermochemical parameters as listed in Table 3, the proposed MBM is able to reasonably capture both surface temperature and mass loss during pyrolysis of a woody biomass particle.

In fact, all the included models behave similarly, despite being constructed in significantly different ways. Nearly identical surface temperature profiles were produced using different models. However, relatively large differences are shown in Figure 3(b) regarding the center temperature. It is however interesting to see that the IBM 1, using a thermal drying model, and the proposed MBM, using a kinetic drying model, predicted very similar drying times (the time until the center temperature rises above the boiling point of water). This further proves the applicability of both models. The difference between model predictions of the drying time mainly stem from the selec- tion of the thermochemical parameters, such as heat capacity and thermal conductivity. Examples of the influence of such parameters are illustrated in Appendix B by varying the thermal conductivity of char. These examples also show that better agreement with experimental data for both the IBM 1 and the proposed MBM can be obtained by choosing different parameter values. Due to the overlapping of drying and devolatilization, the mass loss of the particle is over 80% at the end of drying in all simulations. Therefore, although different center temperatures were obtained after the end of dry- ing as shown in Figure 3(b), mass loss profiles are rather similar, with the exception of the simplified IBM 2 as indicated in Figure 3(c).

500 750 1000

Temperature (K)

(a)

0 50 100 150 200 250

Time (s) 500

750 1000

Temperature (K)

(b)

0 50 100 150 200 250

Time (s) 0.0

0.2 0.4 0.6 0.8 1.0

Mass loss ratio (­)

(c)

This paper IBM 1 [5,6]

IBM 2 [18]

Comprehensive MBM 1 [14]

Exp. [14]

Figure 3: Surface temperature profile (a), center temperature profile (b), and mass loss profile (c) during pyrolysis of a beech particle using the parameters in Table 3. Better agreement for the current MBM and IBM 1 can be obtained with other parameters, as shown in Appendix B. (d= 20 mm,Tw=Tg= 1098 K,M Cwb= 20%,φN2= 100%)

Similarly, comparisons between model predictions and experimental data of combustion of a poplar particle are shown in Figure 4. It can be seen in Figure 4(a) that relatively big differences exist between the simulated surface temperatures and the experimental data. All models give over-prediction of the surface temperature at the early stage of the conversion process (less than around 20 s). This disagreement likely arises from the experimental uncer- tainties, as explained by the authors of the original work [15]. Since the sur- face temperature was measured by a thermocouple buried next to the surface, the reading from the thermocouple might be slightly lower than the actual value. All models underpredict surface temperatures after around 10–20 s.

The lack of heat flux from the volatile flame perhaps contributes most to the underprediction of the surface temperature. Unfortunately, it is difficult to

accurately describe the process of volatile combustion using a single-particle model only. This issue, however, can be resolved by the coupling with a CFD solver. Additionally, authors of the experimental work have stated that the bead of the thermocouple might lose contact to the particle surface and ex- pose to the volatile flame during the tests [15], which potentially resulted in higher temperature readings. It should be noted that the surface tempera- ture calculated from the proposed MBM decreases slightly at around 50 s as shown in Figure 4(a), which is a result of the completion of the drying pro- cess. If a particle contains moisture, the center temperature of the particle stays relatively low. The devolatilization front moves gradually from the hot surface to the cold drying front, thus gently releasing volatile gases. A small amount of volatile gases does not consume all the oxygen transported to the particle surface. Therefore, there is still oxygen left for the exothermic char oxidation, which elevates the particle surface temperature. However, after the depletion of moisture inside the particle, the temperature of the particle increases dramatically. A large amount of volatile gases is released accord- ingly, and this causes a complete consumption of the oxygen that terminates further char oxidation temporarily. Without the heat released from char ox- idation, the surface temperature decreases slightly. Again, despite the large differences between measured and calculated temperatures, the mass loss ra- tios predicted by all models (apart from the IBM 2) agree well with the experiment as shown in Figure 4(c). This phenomenon is further discussed in Section 4.

500 1000 1500

Temperature (K)

(a)

0 20 40 60 80 100

Time (s) 500

1000 1500

Temperature (K)

(b)

0 20 40 60 80 100

Time (s) 0.0

0.2 0.4 0.6 0.8 1.0

Mass loss ratio (­)

(c)

This paper IBM 1 [5,6]

IBM 2 [18]

Comprehensive MBM 2 [17]

Comprehensive MBM 3 [15]

Exp. 1 [15]

Exp. 2 [15]

Figure 4: Surface temperature profile (a), center temperature profile (b), and mass loss profile (c) during combustion of a poplar particle. (d= 9.5 mm, Tw= 1273 K, Tg= 1050 K,M Cwb= 40%,φO2= 21%, φN2= 79%)

It can be seen from Figure 4(b) that all models produced a distinct tem- perature plateau until around 40–50 s owing to the drying process. In fact, it is promising to see that all models predict a similar drying time. On the other hand, the measured temperature rises gradually after only 10 s. Such behav- ior may be caused by heat conduction along the wire as detailed in a recent study [47]. This influence should be more severe in the combustion case [15]

than in the pyrolysis case [14] due to the smaller particle and higher gas tem- perature. In order to further test the proposed MBM under such conditions (small particle and high gas temperature), another experiment [47] was used for validation, which contains measured temperature profiles for pyrolysis of boxwood particles with three different diameters (d= 0.003, 0.005, 0.008 m). In the simulation, the ambient gas temperature was set to 1355 K, which is same as the measured gas temperature at the end of the experiment (50

s after the start of the experiment), whereas the reactor wall temperature was set to 723 K, as given in the experimental work [47]. Figure 5 shows the comparison between measured center temperature and calculated values using the proposed MBM. Note that all the shown experimental data have been corrected to minimize the effect of heat conducted along the thermo- couple wires [47]. It is clear that a better agreement is achieved compared to Figure 4(b).

0 10 20 30 40 50

Time (s) 400

600 800 1000 1200 1400

Particle center temperature (K)

3 mm 5 mm 8 mm This paper

Exp. [47]  

Figure 5: Center temperature profile during pyrolysis of boxwood particles using the parameters in Table 3. (d = 3, 5, 8 mm, Tw = 723 K, Tg = 1355 K, M Cwb = 0%, φN2= 100%)

4. Results and discussion

In order to better understand the effects of spatial and temporal resolution on the performance of the MBM, the two validation cases described in Section 3 (the pyrolysis of a 20 mm beech particle and the combustion of a 9.5

mm poplar particle) were investigated further. The influence of the spatial resolution of the MBM was investigated by varying the numbers of cells, and the influence of the temporal resolution was studied by changing the time step. In addition, two non-dimensional factors, relative error (RE, %) and relative calculation time (RCT, -), are defined accordingly as:

RE = 50%

v u u u t 1 n

n=b^ttotal,conc X

t=1,2,...,n

m_p(t)−m_p,BL(t) m_p,BL(t)

!2

+50%

v u u u t 1 n

n=b^ttotal,conc X

t=1,2,...,n

T_s,sur(t)−T_s,sur,BL(t) T_s,sur,BL(t)

2

(5)

RCT = ttotal,cal

t_total,con (6)

where m_p(t) and m_p,BL(t) are the particle masses (kg) calculated from a test case and a baseline case at time t respectively, T_s,sur(t) and T_s,sur,BL(t) are the particle surface temperatures (K) calculated from a test case and a baseline case at time t, respectively, t_total,con is the total conversion time of a particle (s), and t_total,cal is the elapsed real time (wall-clock time) for the calculation (s). When coupled with a CFD solver, the MBM exchanges energy and mass with the surrounding environment. Accurate estimations of surface temperature and mass loss of the particle are crucial to calculate fluxes of energy and mass coming out of the MBM. Therefore, a combined error function defined by both surface temperature and mass of the particle is constructed in Equation (5). The baseline case referred to is a case configured with very fine spatial and temporal resolutions. Details are given in the following sections.

4.1. Effect of the spatial resolution

Figure 6 shows predicted temperatures and mass loss ratios during com- bustion of a 9.5 mm poplar particle using the MBM with the same time step of 1×10⁻⁴ s but different degrees of spatial discretization (i.e. different numbers of cells). Note that the case denoted as 10^† was configured slightly different from the others, in that combustion of char was not allowed until all the dry wood in the particle had become char. The surface temperature profiles look overall similar in Figure 6(a). However, some oscillations can be observed for the cases with coarse spatial resolution. These oscillations are caused by intermittent char conversion in the outermost cell. Since the volatile flame is not accounted for in the particle model when de-coupled from a CFD solver, the available oxygen for char conversion is instead estimated from the amount of hydrogen being released (cf. Section 2.3). As indicated from the surface temperature of the 10^† case before 50 s, the outermost cell was heated to a high temperature even without char conversion. This sit- uation makes the char conversion in this cell very sensitive to the available oxygen concentration. Given the small thermal mass of this cell, even a small amount of oxygen can cause a big increase in the local temperature, and consequently also in the surface temperature. Similar oscillating behav- ior has also been reported previously using a different particle model [49].

Although the oscillation may appear dramatic, its influence on the overall particle conversion is nearly negligible as shown in Figure 6(c). Moreover, this oscillating phenomenon would be greatly reduced in a CFD simulation of a large furnace, in which the oxygen concentration is influenced by a collec- tive of particles and the homogenous gas phase combustion is simultaneously

accounted for.

500 1000 1500

Temperature (K)

(a)

0 20 40 60 80 100

Time (s) 500

1000 1500

Temperature (K)

(b)

0 20 40 60 80 100

Time (s) 0.0

0.2 0.4 0.6 0.8 1.0

Mass loss ratio (­)

(c)

10 10 30 50100 200Exp. 1 [15]

Exp. 2 [15]

Figure 6: Effect of the spatial resolution on surface temperature profile (a), center temper- ature profile (b), and particle mass loss ratio (c) during combustion of a poplar particle.

(d= 9.5 mm,Tw= 1273 K, Tg= 1050 K,M Cwb= 40%,φO2= 21%, φN2= 79%)

With an increased degree of spatial resolution, calculated surface temper- atures seem to converge to the same values. On the other hand, the particle center temperature varies more with the number of cells. This is due to that the center temperature shown in Figure 6(b) is in fact the temperature of the innermost cell. Since the particle is heated externally, the larger the inner cell is (i.e. the lower the degree of spatial discretization is), the sooner the center temperature appears to rise inside the particle. As shown in Figure 6(c), the effects of spatial resolution on the mass loss ratio are very small:

the predicted mass loss evolution is well characterized already at a resolution of 10 cells.

Both temperature and heating rate may affect product distribution dur- ing the devolatilization process. As shown in Figure 6, certain variations can

be observed in the temperature prediction with varying degrees of spatial dis- cretization. Since competing reactions have been implemented for resolving the devolatilization process, the spatial discretization may also have an influ- ence on the distribution of the devolatilization products. Figure 7 shows the normalized cumulated mass for gas, tar, and char under the same condition as used in Figure 6 as obtained with particle resolutions of 10, 50, and 200 cells. It seems that nearly identical product distributions are obtained. The heating rate of the particle is relatively low (estimated to be below 100 K/s during the devolatilization process) under the tested condition, and the slight differences in the temperature prediction caused by varying spatial resolution are too small to influence the product distribution.

0 20 40 60 80 100

Time (s) 0.0

0.2 0.4 0.6 0.8 1.0

Normalized mass (­)

1050 200

Gas   Tar Char   

Figure 7: Effect of the spatial resolution on distribution of the devolatilization products during combustion of a poplar particle. The product mass is normalized by the initial particle dry mass. (d= 9.5 mm,Tw= 1273 K, Tg = 1050 K,M Cwb = 40%,φO2= 21%, φN2= 79%)

By using Equations (5) and (6), quantitative differences of computational

accuracy and cost between calculations with various degrees of discretization were obtained, as shown in Figure 8. Besides the combustion test case, the pyrolysis of a 20 mm particle described in Section 3 was also simulated. The relative errors shown in Figure 8 were calculated by comparing with the base- line cases, respectively for the combustion (200-cell spatial discretization) and pyrolysis (300-cell spatial discretization) conditions. The low relative errors again demonstrate that the predictability of MBM model is not sensitive to the spatial resolution at the tested conditions. On the contrary, the calcu- lation time is heavily influenced by the number of cells (note a logarithmic scale is used). It can be estimated from Figure 8 that when the number of cells is lower than 25, the RCT is lower than 1 or, in other words, the time used for executing the code is shorter than the actual physical conversion process simulated.

The computational efficiency of the MBM is further analyzed by timing different operations, i.e. solving the matrix system and all other operations.

The effect of the spatial resolution on the time used for the matrix solver and the other operations are shown in Figure 9. As described in Section 2.2, the system of linear equations is directly solved by LU factorization to obtain the temperature of each discretized cell. The computational cost of solving such a system scales with the size of the system to the power of three. Consequently, as shown in Figure 9, the time used for solving matrix is almost proportional to cube of the number of cells. The rest of the operations, such as calculations of the reaction rates and updating the mass, need to be performed the same number of times for each individual cell at any given time. Therefore, the time used for those operations is linearly correlated to the number of cells.

0 100 200 300 Cell number (­)

0 1 2 3 4 5 6

Relative error (

RE

10¹

10⁰

10¹

10²

10³

Relative calculation time (

RC T

RE Combustion RE Pyrolysis RCT Combustion RCT ^Pyrolysis

Figure 8: Effect of the spatial resolution on relative error and relative calculation time.

(combustion: poplar,d= 9.5 mm,Tw= 1273 K,Tg= 1050 K,M Cwb= 40%,φO2= 21%, φN2= 79%; pyrolysis: beech,d= 20 mm,Tw=Tg= 1098 K,M Cwb= 20%,φN2= 100%)

0 50 100 150 200

Cell number (­) 0

500 1000 1500 2000 2500 3000 3500

Time used for matrix solver (s) 100

200 300 400 500 600 700

Time used for other operations (s)

Matrix solver

y = 0.0004x³

Others

Figure 9: Effect of the spatial resolution on computational efficiency. (poplar, d = 9.5 mm,Tw= 1273 K,Tg= 1050 K, M Cwb= 40%, φO2= 21%,φN2= 79%)

4.2. Effect of the temporal resolution

Figure 10 shows the effect of the temporal resolution (time step) on the prediction of temperature and mass loss using the MBM. The same condition of the combustion of a 9.5 mm poplar particle was used. All the calculations shown in Figure 10 were obtained using the same spatial discretization of 50 cells, which is sufficient to obtain reasonable results under the tested condi- tions as illustrated in Section 4.1. As shown in Figure 10, most of the profiles look nearly identical, except the ones calculated with time steps larger than 1×10⁻² s. Since a fully implicit finite-volume method is used in the MBM, an unconditionally stable solution should be expected. However, as mentioned before, solutions calculated with time steps larger than 1×10⁻² s differ sig- nificantly to other results. This suggests that some physical processes limit the largest tolerable time step with respect to solution accuracy. Detailed illustrations of this phenomenon are given in the next section.

The effect of the temporal resolution on computational accuracy and cost for both combustion and pyrolysis conditions are shown in Figure 11. It is not surprising that the computational cost is inversely proportional to the time step. When the time step is larger than 5×10⁻⁴ s, the computational time is shorter than the physical conversion time at the tested conditions and current spatial discretization. The relative errors were estimated by comparing to the base cases calculated with a time step of 1×10⁻⁵ s. For both tested combustion and pyrolysis conditions, a very small relative error can be obtained by using a time step smaller than 1×10⁻² s.

500 1000 1500

Temperature (K)

(a)

0 20 40 60 80 100

Time (s) 500

1000 1500

Temperature (K)

(b)

0 20 40 60 80 100

Time (s) 0.0

0.2 0.4 0.6 0.8 1.0

Mass loss ratio (­)

(c)

1 × 10¹s 5 × 10²s 1 × 10²s 1 × 10³s 1 × 10⁴s 1 × 10⁵s Exp. 1 [15]

Exp. 2 [15]

Figure 10: Effect of the temporal resolution on surface temperature profile (a), center temperature profile (b), and particle mass loss ratio (c) during combustion of a poplar particle. (d= 9.5 mm,Tw= 1273 K,Tg= 1050 K,M Cwb= 40%,φO2= 21%,φN2= 79%)

10⁵ 10⁴ 10 ³ 10 ² 10 ¹ 10⁰ Time step (s)

0 10 20 30 40

Relative error (

RE

10³ 10² 10¹ 10⁰ 10¹

Relative calculation time (

RC T

RE Combustion

RE Pyrolysis

RCT Combustion

RCT Pyrolysis

Figure 11: Effect of the temporal resolution on relative error and relative calculation time.

(combustion: poplar,d= 9.5 mm,Tw= 1273 K,Tg= 1050 K,M Cwb= 40%,φO2= 21%, φN2= 79%; pyrolysis: beech,d= 20 mm,Tw=Tg= 1098 K,M Cwb= 20%,φN2= 100%)

4.3. Optimal spatial and temporal resolutions

Despite the unconditionally stable nature of the implicit scheme used to solve Equation (1), small time steps are required to ensure accuracy. In order to find the optimal spatial and temporal resolutions, three scenarios (inert heating, pyrolysis, and combustion) were tested based on the condition of the combustion of a 9.5 mm poplar particle described in the previous sections.

The drying, pyrolysis, and char conversion were disabled in the inert heating cases, whereas the char conversion was turned off in the pyrolysis cases. The surface temperature profiles for the inert heating scenario with various level of spatial and temporal resolutions are shown in Figure 12.

0 100 200 300

Time (s) 0

200 400 600 800 1000 1200

Surface temperature (K)

1 × 10⁰s 1 × 10⁴s 1020

50100 200

10  50  200

Figure 12: Effect of spatial (cells: 10 – 200) and temporal resolutions (time step: 1×10⁰ and 1×10⁻⁴ s) on prediction of surface temperature during heat transport process.

As shown in Figure 12, when the time step is small (1×10⁻⁴ s), the obtained surface temperatures are identical for the cases with varying spatial resolutions. Due to the similarity, only three cases were plotted. However,

large discrepancies can be noted from the cases calculated using the time step of 1×10⁰ s. Counterintuitively, the results deviate more with an increased number of grid points. This observed phenomenon can be attributed to the comparison between the time step and the time scale of conductive heat transfer (t_HT) inside the particle, which is defined by the equation below:

t_HT = (length scale)²

thermal diffusivity (7)

The length scale of each discrete volume can be estimated by the distance between two adjacent grid points. Therefore, when the particle size remains constant, this length scale reduces as the number of grid points increases, and so does the associated time scale t_HT for cell-to-cell heat conduction. If t_HT is smaller than the time step used in the solution procedure, the heat transfer process cannot be resolved accurately. Figure 13 shows the heat transfer timescale of each discretized cell for all three scenarios (time step = 1×10⁻⁴ s and cell count = 50). The heat transfer time scales for all cells increase initially due to the decreasing thermal diffusivity when temperature is elevated. Since the spatial discretization is performed by dividing the particle into equal mass cells (i.e. equal volume cells) in its radial direction, the length scale of the cell decreases from particle center to particle surface.

For the inert heating case shown in Figure 13(a), it is apparent that majority of the timescales are smaller than 1×10⁰ s. Consequently, a large deviation can be found for the results calculated using the time step of 1×10⁰ s and the spatial discretization of 50. As the cell number increases, even smaller time steps are required, and thus more severe underprediction of the surface temperatures were obtained as shown in Figure 12.