CFD-DEM Simulation of Biomass Pyrolysis in Fluidized-Bed Reactor with a Multistep Kinetic Scheme

Abstract

The pyrolysis of biomass in a fluidized-bed reactor is studied by a combination of a CFD-DEM algorithm and a multistep kinetic scheme, where fluid dynamics, heat and mass transfer, particle collisions, and the detailed thermochemical conversion of biomass are all resolved. The integrated method is validated by experimental results available in literature and a considerable improvement in predicting the pyrolysis product yields is obtained as compared to previous works using a two-fluid model, especially the relative error in the predicted tar yield is reduced by more than 50%. Furthermore, the evolution of light gas, char and tar, as well as the particle conversion, which cannot easily be measured in experiments, are also revealed. Based on the proposed model, the influences of pyrolysis temperature and biomass particle size on the pyrolysis behavior in a fluidized-bed reactor are comprehensively studied. Numerical results show that the new algorithm clearly captures the dependence of char yield on pyrolysis temperature and the influence of heating rate on light gas and tar yields, which is not possible in simulations based on a simplified global pyrolysis model. It is found that, as the temperature rises from 500 to 700 °C, the light gas yield increases from 17% to 25% and char yield decreases from 22% to 14%. In addition, within the tested range of particle sizes (<1 mm), the impact on pyrolysis products from particle size is relatively small compared with that of the operating temperature. The simulations demonstrate the ability of a combined Lagrangian description of biomass particles and a multistep kinetic scheme to improve the prediction accuracy of fluidized-bed pyrolysis.

1. Introduction

The thermochemical conversion of biomass has become more and more important in the renewable energy industry in recent years, an industry which contributes above 10% of the total power generation nowadays in European countries [1]. The major advantage of biomass over coal lies in its CO₂-neutral, renewable characteristics and large reservation around the world. Among various routes of biomass utilization, fluidized-bed pyrolysis is one of the most promising technologies to produce bio-oil and other high-quality fuels [2,3,4,5]. These fuels can be used for power generation or serve as an alternative to fossil fuels in the transportation sector. Recently, the conversion of biomass to valuable energetic products such as hydrogen is a good example [6,7]. Contemporary biomass conversion technology in fluidized-bed reactors also include oxy-combustion, chemical looping combustion, etc., [8,9,10], which could remarkably improve the conversion efficiency and reduce CO₂ emissions that otherwise constitute a major cause of global warming in traditional coal combustion [11,12,13]. In addition, according to recent studies, by either torrefication or blending with coal, the energy conversion efficiency can reach a much higher value than in traditional methods because of the diverse components of biomass materials [14,15,16]. However, due to the inherent complexity of biomass pyrolysis, the application to large-scale fluidized-bed reactors still faces many challenges [17,18]. For example, tar removal and pollutant formation are some of the most serious issues, which both reduce the biomass conversion efficiency and also cause air pollution that affects human health and further influences global warming [19,20].

In the last few decades, many kinetic models have been developed to explore the mechanism of biomass pyrolysis [21,22]. The earliest global type pyrolysis models are still widely used until today due to their high computational efficiency [23], especially in the simulation of large-scale reactors [24,25,26]. However, these models are too rough to capture the variation of pyrolysis yields with the change of operating conditions. To solve this problem, Miller and Bellan proposed the so-called competitive mechanism [27], where the formation of light gas and char competes with that of tar by using different kinetic rates. This model is applicable for each component of biomass (cellulose, hemicellulose, and lignin) and allows for the variation of product with pyrolysis conditions. However, the model is only able to predict the general trends of pyrolysis yields under different temperatures, as detailed information for a specific species cannot be obtained. To get a more accurate prediction of pyrolysis yields, Niksa established bio-FLASHCHAIN pyrolysis model [28], where the thermal degradation of biomass is modeled as four basic reactions: bridge scission, spontaneous charring, bimolecular recombination, and char decomposition. The coefficients of species yield in this model can also be determined from experimental measurements.

During the thermochemical conversion of biomass, tar generation is one of the most challenging issues [29,30,31], which could cause serious technical problems in real applications such as degradation and blocking of reactors [32] and erosion of equipment [33]. Although a lot of efforts have been made, the mechanism of tar generation is still an open problem. Experimental works in literature have shown that tar contains a wide range of chemical species (C₃H₆O₂, C₆H₅OH, C₇H₈…) and that these species vary a lot depending on the operating conditions and biomass types [34,35]. However, in theoretical studies, simplifications must be made to get a general understanding of the mechanism of tar generation. The multistep kinetic scheme developed by Ranzi’s group in recent years is one such mechanism [36,37]. The main idea of the multistep scheme is to model the pyrolysis of biomass with a series of parallel sub-mechanisms applied to cellulose, hemicellulose, and lignin, where light gas and tar are lumped to twenty representative species. The medium complexity and the detailed species yield information of the model make it very popular in the analysis of biomass pyrolysis. For example, Gentile et al. [38] studied the pyrolysis behavior of a single biomass particle by using the multistep kinetic model. Good agreement between simulation and experimental data was obtained for both light gas and tar species. In addition, the work in our group also demonstrates the universality of the multistep kinetic model in predicting the unsynchronized evolution for different light gas and tar species [39]. Due to the good performance of the multistep kinetic model, researchers have recently begun to use it in the study of reactor scale biomass pyrolysis. Sommariva et al. [40] revealed the important influence of inlet temperature on the start-up performance of a fixed-bed reactor from numerical simulations with the multistep kinetic model. Peng et al. [41] used the multistep model to investigate the influence of operating temperature and the mass fraction of basic components on fast pyrolysis in fluidized-bed reactors. They found that the pyrolysis reaction rate increased remarkably when the operating temperature exceeds 400 °C. In addition, the tar yield reached a maximum of 70% at 500 °C. Shi et al. [42] also used the multistep kinetic model to investigate the pyrolysis yields in a screw reactor. It was observed that both reactor temperature and char yield showed a reasonable agreement with the experimental measurement. Ranganathan and Gu [43] simulated fluidized-bed pyrolysis with different kinetic schemes and found that the multistep kinetic model performs much better than other kinetic models.

In addition to the accuracy of the adopted pyrolysis model, the multiphase flow algorithm is another important issue in reactor scale pyrolysis simulations. As has been pointed out by Ranzi et al. [44], biomass pyrolysis is a multi-scale problem, where the interphase transport phenomena both at particle scale and reactor scale should be carefully considered. In general, there are two types of multiphase flow algorithms to deal with interphase exchange: the Eulerian–Eulerian method and the Eulerian–Lagrangian method [45]. The Eulerian–Eulerian method is especially suitable for large scale simulations with huge numbers of particles due to a simplified modeling of the particle collisions and the lack of detailed information about the environment surrounding the individual particles [46]. Instead of tracking each particle in the reactor, a pseudo-stress model is generally used to simulate particle interactions which remarkably reduces the computational cost. However, due to the averaging over many particles, the heat/mass transfer between the gas phase and the particle phase can only be approximately accounted for. The Eulerian–Lagrangian method differs from the Eulerian–Eulerian approach in that it models particle–particle interaction by using Newton’s second law of motion, thus allowing for a more accurate consideration of interphase exchange [47,48]. For example, it is straightforward to take the dynamic effects of particle shrinkage and varying particle sizes into account, which is much more difficult in Eulerian–Eulerian algorithms [49,50,51]. Existing studies also illustrate that the Eulerian–Lagrangian method predicts more realistic fluidization phenomena such as segregation rate and bubble structures [52,53]. However, most of the above-mentioned simulations of fluidized-bed pyrolysis are based on the Eulerian–Eulerian method. Therefore, the influence of the accuracy of the interphase exchange on the prediction of the chemical kinetics in these works is still not well known.

From the above-mentioned works, one can see that the existing numerical studies for biomass pyrolysis in fluidized-bed reactors are mainly based on Eulerian–Eulerian and Eulerian–Lagrangian multiphase flow models. The Eulerian–Lagrangian model is usually more accurate than a Eulerian–Eulerian model, and it is also more time-consuming. Therefore, one-step pyrolysis models such as global reaction mechanisms are usually used for simplicity and to save computational time. However, the simplified one-step pyrolysis model is incapable of describing many of the intricate details of the biomass conversion process, and is therefore far less accurate in the study of some of the key issues related to tar evolution and pollutant formation (e.g., NOx and SOx) [14,20]. On the other hand, the studies with a detailed pyrolysis model are mainly based on Eulerian–Eulerian models. The lack of an assessment of a combined Eulerian–Lagrangian multiphase flow model and a detail pyrolysis model contributes to the research gap and thus motivates the current work.

To obtain a comprehensive understanding of biomass pyrolysis in fluidized-bed reactors, this paper is aimed at developing a detailed algorithm based on the combination of the Eulerian–Lagrangian approach denoted CFD-DEM (discrete element method) and the multistep pyrolysis kinetics, where both the evolution of syngas and tar species as well as the interphase heat/mass transfer can be accurately considered. The performance of the new algorithm will be validated and compared with existing experimental results and numerical predictions from a Eulerian–Eulerian multiphase flow model with a multistep pyrolysis scheme. Fluidization characteristics, such as the particle temperature distribution, which are not easily obtained in experiments will also be discussed. The current work could improve the mechanism exploration of biomass pyrolysis in fluidized-bed reactors and also lay a foundation for future studies of more complex issues such as tar removal and pollutant formation. The paper is organized as follows: Section 2 presents the mathematic model. The implementation of the algorithm is presented in Section 3. After that, model validation and discussions are carried out in Section 4. Finally, conclusions are drawn in the last section.

2. Computational Model

Biomass pyrolysis in a fluidized-bed reactor includes several simultaneous physical and chemical conversion processes. The final pyrolysis products depend to a large extent on the heating process of each biomass particle. To capture the detailed conversion behavior of the biomass particles, the current pyrolysis algorithm is established based on a CFD-DEM model developed in our previous works [54,55,56]. In this section, the governing equations of gas and particle phases, the interphase heat/mass exchange and the multistep pyrolysis kinetics are briefly presented. The gas phase governing equations which include the mass, momentum, energy, and species transport equations as well as the k-ε turbulence model, are listed in

Table 1. Gas phase governing equations.

$\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}\right)+\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}\right)=S_{\mathrm{p},\mathrm{m}}$

$\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}\right)+\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}\right)=-\nabla p+\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\mathbf{\tau }}_{\mathrm{eff}}\right)+{\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}\mathit{g}+S_{\mathrm{p},\mathrm{mom}}$

$\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}E\right)+\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}({\rho }_{\mathrm{g}}E+p)\right)=\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\alpha }_{\mathrm{eff}}\nabla h_{\mathrm{sen}}\right)+S_h+S_{\mathrm{p},h}+S_{\mathrm{rad}}$

$\frac{\partial }{\partial t}\left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}{Y}_i\right)+\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}{\mathit{u}}_{\mathrm{g}}{Y}_i\right)=\nabla \cdot \left({\epsilon }_{\mathrm{g}}{\rho }_{\mathrm{g}}{D}_{\mathrm{eff}}\nabla Y_i\right)+S_{\mathrm{p},Y_i}+S_{Y_i}$

$\begin{array}{l}\frac{\partial }{\partial t}({\epsilon }_g{\rho }_{g}k)+\nabla \cdot ({\epsilon }_g{\rho }_g{\mathbf{u}}_{g}k)=\nabla \cdot ({\epsilon }_g({\mu }_g+\frac{{\mu }_t}{{\sigma }_k})\nabla k)+{\epsilon }_g{G}_k-{\epsilon }_g{\rho }_g\epsilon \\ \frac{\partial }{\partial t}({\epsilon }_g{\rho }_g\epsilon )+\nabla \cdot ({\epsilon }_g{\rho }_g{\mathbf{u}}_g\epsilon )=\nabla \cdot ({\epsilon }_g({\mu }_g+\frac{{\mu }_t}{{\sigma }_{\epsilon }})\nabla \epsilon )+{\epsilon }_g\frac{\epsilon }{k}(C_{\epsilon 1}{G}_k-C_{\epsilon 2}{\rho }_g\epsilon )\\ {\mu }_t={\rho }_g{C}_{\mu }\frac{k^2}{\epsilon }\end{array}$

[55]. Note that in the energy equation, P-1 model is applied to compute the radiation term [57]. For the k-ε turbulence model, the values of the constant parameters are set as follows: Cε1 = 1.44, Cε2 = 1.92, Cμ = 0.09, σk = 1.0 and σε = 1.3. Further detailed information about the gas phase evolving equations can be found in our earlier publication [55].

For both particle phases (biomass and bed material), each particle is treated as a discrete entity.

Table 2. Particle phase governing equations.

Mass:

$\frac{d{m}_p}{dt}=\frac{d{m}_{vapor}}{dt}+\frac{d{m}_{devol}}{dt}$

Momentum:

$\begin{array}{l}{m}_p\frac{d{\mathit{v}}_p}{dt}={\mathit{f}}_g+{\mathit{f}}_c+m_p\mathit{g}\\ I_p\frac{d{\mathbf{\omega }}_p}{dt}={\mathbf{T}}_p\end{array}$

Energy:

$m_p{c}_p\frac{d{T}_p}{dt}=h{A}_p(T_g-T_p)+\frac{e_p{A}_p}{4}(G-4\sigma T_p^4)+Q_p$

Drag correlation:

$\begin{array}{l}{\mathit{f}}_g=\frac{V_p\beta }{{\epsilon }_p}({\mathbf{u}}_g-{\mathit{v}}_p)\\ \beta =\{\begin{array}{ll}150\frac{{\epsilon }_p^2{\mu }_g}{{\epsilon }_g^2{d}_p^2}+1.75\frac{{\epsilon }_p{\rho }_g}{{\epsilon }_g{d}_p}\left|{\mathbf{u}}_g-{\mathit{v}}_p\right|& {\epsilon }_g<0.8\\ 0.75C_d\frac{{\epsilon }_p{\rho }_g}{d_p}\left|{\mathbf{u}}_g-{\mathit{v}}_p\right|{\epsilon }_g^{-2.65}& {\epsilon }_g\ge 0.8\end{array}\\ C_d=\{\begin{array}{ll}\frac{24}{{\mathrm{Re}}_p}(1+0.15R{e}_p^{0.687})& {\mathrm{Re}}_p<1000\\ 0.44& {\mathrm{Re}}_p\ge 1000\end{array}\\ {\mathrm{Re}}_p={\epsilon }_g{\rho }_g{d}_p\left|{\mathbf{u}}_g-{\mathit{v}}_p\right|/{\mu }_g\end{array}$

presents the particle phase governing equations [55], which also include mass, momentum and energy evolving equations as well as the Gidaspow drag correlation used to calculate interphase momentum exchange [58]. In the energy equation, it is assumed that the temperature throughout the particle radius is uniform during the heating and conversion process because of the small particle size. Moreover, every biomass particle undergoes a continuous shrinkage with a constant particle density. During fluidization, both the collisions between particles (e.g., biomass and sand particles) and the collisions between particles and reactor wall are simulated by a soft-sphere model [59].

For modeling biomass pyrolysis, a multistep kinetic scheme is implemented in the CFD-DEM model. The applicability of the multistep kinetic model has been demonstrated in our previous work [39], in which the pyrolysis of a single biomass particle was investigated and a good agreement between the simulated yields of char, tar, and light gas with the experimental measurements was achieved. However, it was also found that there is still a small amount of intermediate species such as G{CO} left in the particle after pyrolysis. Note that these intermediate virtual species only influence the release rate of the corresponding gas species without changing the total pyrolysis yield and without taking part in the subsequent char reactions. Therefore, in the present paper, the activation energies for the thermal decomposition of the intermediate species G{H2}, G{CO}, G{CH4}, G{COH2}, G{C2H4} and G{CH3OH} (please see ref [37]) are slightly changed to 45, 40, 40, 40, 50 and 40 kcal/mol, respectively.

3. Numerical Algorithm

The integrated algorithm has been implemented into the open source CFD toolbox OpenFOAM [60], and the final flow chart is shown in Figure 1. The particle motion is solved by a first-order Euler time integration scheme. For the gas phase governing equations, the standard finite-volume method in OpenFOAM is utilized for spatial discretization and an Euler scheme is used for time discretization. For all the simulations, a small time step of 1.0 × 10⁻⁵ s is adopted to ensure temporal convergence. During each time step, the particle mass loss is calculated by the moisture evaporation and the multistep pyrolysis submodels. Meanwhile, the particle properties such as temperature, diameter, heat capacity, and components are correspondingly updated. The interphase coupling is resolved by the various interphase source terms (i.e., Sp,m, Sp,mom, Sp,h, and Sp,Yi). In addition, a parallelization technique with Massage Passing Interface (MPI) is also used to accelerate the simulation.

4. Results and Discussions

In this section, the established pyrolysis model is first validated with experimental data reported in literature. The product distributions of both light gas and tar species that are not easily obtained in the experiment are also analyzed. Then, the influence of pyrolysis temperature and biomass particle size on the pyrolysis behavior is systematically explored.

4.1. Validation

In our previous works, the CFD-DEM algorithm has already been systematically validated [54,55,56]. Good performance was achieved in the parameter sensitivity analysis for a series of operation conditions such as temperature, gasification agent, and biomass particle size. In addition, in a different work, the multistep pyrolysis model has also been demonstrated to have a high accuracy in predicting pyrolysis yields of tar, char, and light gases [39]. Moreover, there are many other works that have shown the validity of the multistep kinetic scheme in fluidized bed pyrolysis simulations [41,43]. For simplicity, the model validation attended here only concerns the integrated CFD-DEM model with the multistep pyrolysis model. To this end, a biomass fluidized-bed pyrolysis experiment carried out by Kantarelis et al. [61] is adopted for validation. The numerical simulation for this experiment by the same research group is used for additional comparisons with the current work. Figure 2 shows a schematic diagram of the reactor. The fluidizing gas is nitrogen with a temperature of 463 °C. The pyrolysis products including light gases (H₂, CH₄, CO, CO₂, and C₂H₄) and tar leave the top outlet boundary and are collected for analysis. Initially, the bed material consists of sand particles with a total number of 12,000. Note that in order to reduce the computational cost, the thickness of the bed is set equal to the particle diameter which means that the bed has only one layer of particles. Moreover, due to the complexity of the screw feeder, the biomass feeding position is adjusted from the sidewall to the reactor centerline with the same feeding height (i.e., 57.5 mm) as that in the experiment. Before introducing the biomass, the simulation is running for one second with only sand particles. Afterwards, biomass particles with a diameter of 850 μm are continuously injected and the simulation is run for 10 s. More details of the simulation setup are presented in

Table 3. Simulation setups of the fluidized-bed reactor.

	Parameter	Value
Reactor	geometry (m)	0.072 × 0.65 × 0.00085
	wall temperature (°C)	500
	inlet gas temperature (°C)	463
	inlet gas flow rate (g/s)	0.05
Sand particle	shape	sphere
	diameter (μm)	850
	density (kg/m3)	2650
Biomass particle	shape	sphere
	diameter (μm)	450, 650, 850
	material	mixture of pine and spruce
	density (kg/m3)	600
	feeding rate (g/s)	0.00835
Resolution	grid	20 × 130 × 1
	time step (s)	1.0 × 10−5

[61,62] and the biomass components are provided in

Table 4. Component mass fraction of the mixed pine and spruce (%).

Biomass		Lignocellulose
moisture	9.8	cellulose	44.0297
ash	0.3067	hemicellulose	21.9245
lignocellulose	89.8933	LIG-H	16.2854
		LIG-O	6.4921
		LIG-C	11.2683

[39,61].

Figure 3 presents the comparison of the predicted pyrolysis products with the experimental data from Kantarelis et al. [61] as well as the Eulerian–Eulerian simulation results of Mellin et al. [62]. Note that the char yield contains carbon, ash, and unconverted biomass. As shown in Figure 3a, the current pyrolysis model gives a better prediction of steam, tar, and char mass fractions than that of Mellin et al., while the light gas yield is also underestimated due to the overpredicted char yield, in which there is still a little amount of intermediate gas species left. Compared with the Eulerian–Eulerian method, the predicted steam yield of the Eulerian–Lagrangian method increases from 11% to 14%, while the difference of tar yield with that of experiment is decreased by more than 50%. Figure 3b reveals that our model performs better in the concentrations of four light gas species (i.e., H2, CO, CO2 and C2H4), although the hydrogen concentration is still overpredicted. In particular, the predicted CO and C₂H₄ yields in this work are improved significantly, with a reduction in the relative error to the experimental data to approximately 30% compared to the Eulerian–Eulerian method. The above comparison indicates that the multistep kinetic model together with the Lagrangian description of the biomass particles improves the prediction of pyrolysis yields. In addition to the comparison, Figure 3c also provides the mole fractions of tar species obtained from the simulations. It is seen that formaldehyde (CH₂O), methanol (CH₃OH), HAA (C₂H₄O₂), and acetone (C₂H₅CHO) are the most abundant products among the light tar species, while xylose monomer (C₅H₈O₄), HMFU (C₆H₆O₃), and LVG (C₆H₁₀O₅) account for the largest portion among the heavy tar species.

4.2. Pyrolysis Phenomena

The previous validation demonstrates that the established algorithm improves the overall prediction of pyrolysis products in biomass fluidized-bed pyrolysis. In this subsection, pyrolysis phenomena of the validation case are further analyzed to depict the instantaneous behavior. Figure 4 shows the particle distributions at different time instants near the end of the simulation. The formation of bubble structures during fluidization is clearly observed. As the bubbles rise, bed materials experience a strong up-and-down movement resulting in an efficient mixing of biomass and sand particles. Above the bed, some biomass particles are more easily thrown into the freeboard due to their lower particle density as compared with the sand particles. After the breakdown of bubble structures, the particles fall back into the bed. When the biomass particles are heated up, the moisture and volatile contents are released into the reactor, which further promotes the local gas phase mixing and contributes to the tendency of the biomass particles to “float” to the surface.

Figure 5 presents the corresponding temperature distributions. At the bottom of the reactor, the temperature is high because the introduced fluidizing gas has a high temperature. In the bed section, the temperature is comparably lower because of the cooling effect of the injected biomass particles. Downstream in the freeboard section, the temperature gradually reaches the operating temperature (i.e., 500 °C) due to the dilution of particles. Generally, the bed temperature decreases from about 470 °C at the bottom to 410 °C towards the bed surface and then increases gradually to 500 °C in the freeboard section. Consequently, the heating process of the biomass particles is actually rather complex, as the particle may move repeatedly through the different regions during fluidization. Moreover, Figure 6 and Figure 7 depict CO and LVG (C₆H₁₀O₅) mass fraction distributions in the reactor, respectively. In the bed section, it is observed that both light gas (CO) and tar species (C₆H₁₀O₅) have a higher concentration near the wall and a lower one in the middle area of the reactor. This effect arises due to the longer retention time of the biomass particles on their way downwards in the near-wall regions as compared to their faster rise velocities in the core where they are lifted by the bubbles. Downstream in the freeboard section, the concentrations of both species become relatively uniform.

4.3. Influence of Pyrolysis Temperature

Pyrolysis temperature is one of the most important operating parameters during biomass fluidized-bed pyrolysis [63]. In this subsection, three temperatures, i.e., 500, 600 and 700 °C, which cover the typical temperature range encountered in fluidized-bed pyrolysis [64], are tested while the other operating parameters remain the same as those of the validation case except for the inlet gas temperature which is set to the pyrolysis temperature. Moreover, the case of 500 °C is considered as the base case. Figure 8 compares the pyrolysis products for the three different temperatures. As shown in Figure 8a, increasing the pyrolysis temperature raises the light gas production from 17% to 25% and decreases char mass fraction from 22% to 14%, while it slightly affects the steam and primary tar concentrations, in which the steam yield shows an increasing trend and the primary tar yield shows a decreasing trend. This demonstrates that the multistep pyrolysis model captures the dependence of the char yield on the pyrolysis temperature very well, unlike a one-step global reaction mechanism where the char yield has to be pre-assumed before the simulation [65]. From Figure 8b, it is found that the concentrations of H₂, CH₄, CO, and CO₂ are moderately reduced when the pyrolysis temperature increases from 500 to 600 °C, while the C₂H₄ volume fraction increases from 0.06 to 0.1. When the temperature further increases to 700 °C, only a small variation is observed for all the light gas species. From Figure 8c, when the pyrolysis temperature rises from 500 to 700 °C, the concentrations of most low-molecular-weight tar species (C1–C3) show an increasing trend while most high-molecular-weight tar species (e.g., C₅H₈O₄, C₆H₁₀O₅, C₇H₈O, and C₁₁H₁₂O₄) exhibit a decreasing trend.

Figure 9 presents the transient production rates of H₂, CH₄, CO, and CO₂ for different pyrolysis temperatures. Increasing the pyrolysis temperature advances the starting point of light gas releases from about 2 s to 1 s. Moreover, for all light gas species, the average steady state production rate shows an apparent increase when the pyrolysis temperature enhances from 500 to 600 °C. Taking CH₄ for example, the average production rate is promoted from 0.12 mg/s at 500 °C to 0.17 mg/s at 600 °C. However, when the pyrolysis temperature further increases to 700 °C, the average production rate only exhibits a slight increase although the amplitude of the fluctuations becomes larger. In addition to the gas products, Figure 10 compares the instantaneous biomass particle distribution at t = 10 s together with the temperature contour for the three temperatures. It is shown that most of the biomass particles undergo a fast heating process and reach the operation temperature very quickly. In addition, the particles lose mass faster at higher temperature and some of them are blown to the top area of the reactor at 700 °C (Figure 10c). Figure 11 shows the average thermochemical conversion processes of all biomass particles for the three different pyrolysis temperatures. The temperature curves in Figure 11a illustrates that the average heating process inside the biomass particles is very sensitive to the change in pyrolysis temperature. At 500 °C, a fast temperature rise occurs before 0.8 s, to be followed by a significant slowdown in the temperature increase. At 600 °C, after the first stage of rapid temperature increase, a second heating stage can be observed from 0.7 s to 1.8 s. Finally, a third stage characterized by a slower temperature rise is seen. At 700 °C, a small plateau period is observed in the heating process before 0.7 s, after which the temperature rise is accelerated again and lasts until 1.1 s.

4.4. Influence of Particle Size

The CFD-DEM model is a good tool to study the influence of particle size on pyrolysis behavior, which is not as easily accounted for by using a two-fluid model [66,67]. In this subsection, pyrolysis simulations with biomass particle size ranging from 450 to 850 μm are compared, while the size of sand particles is kept at 850 μm. The other operating conditions such as pyrolysis temperature and biomass material are the same with that of the base case. Figure 12 shows the pyrolysis product distributions. It is found that the char yield is slightly promoted with increasing biomass particle size, which results in a reduction in light gas and tar yields. This effect is likely caused by a lower heating rate inside the biomass particles, which influences the final pyrolysis products. Therefore, it is demonstrated that the multistep model can also capture the dependence of pyrolysis yields on the heating rate, unlike a traditional one-step global reaction mechanism [68]. In addition, the relative productions of the light gases are almost the same for the different particle sizes. In addition, the tar yield in Figure 12c illustrates that the relative production of light tar species (C1–C3) is similar for the three cases, except that CH2O and CH3OH productions slightly decrease with increasing particle size. Among heavy tar species, the C5H8O4 and C6H10O5 production exhibit a small increase for larger biomass particles. These changes are most likely caused by the variation in pyrolysis reactions in the multistep kinetic model due to the different heating processes. Note that there are basically two parallel reactions with different activation energies and pre-exponential factors for each components of the biomass [37]. As a result, the relative reaction rate depends strongly on the heating rate.

Figure 13 presents the production rate of four light gas species. It is clearly shown that, the volatile release for smaller biomass particles reaches steady state much earlier due to a faster heating process inside the particles. At the steady state, the release rates for different sized biomass particles are close to each other due to the same fuel feeding rate. Finally, Figure 14 presents the biomass particle distributions for the three cases. With reducing particle size, the number of biomass particles increases and the particles become lighter. As a result, some of the initially 450 μm sized particles are blown into the freeboard region after reaching the operation temperature. For a particle diameter larger than 650 μm, the biomass particles undergo a slower heating process and most of them stay inside the sand bed until the end of pyrolysis. Figure 15 shows a detailed comparison of the heating process and the mass loss history for these three cases. It is observed that reducing the particle size remarkably accelerates the heating and mass loss rates of the particle. However, the final conversion stays at the same level of about 0.84 for all particle sizes.

5. Conclusions

A comprehensive simulation algorithm for biomass pyrolysis in a fluidized-bed reactor is established based on a previously well-validated CFD-DEM model developed in our group. With the new algorithm, the influence of variations in the operating conditions, such as pyrolysis temperature and biomass particle size, is studied systematically. The main conclusions of the paper are as follows:

(1)

The novelty of the integrated algorithm lies in the utilization of a multistep kinetic scheme for the detailed simulation of particle pyrolysis, thus allowing for a practical consideration of product variation with the change of operating conditions. The current work is the first step for future studies of more complex issues such as tar removal and pollutant reduction.

(2)

The integrated CFD-DEM and multistep pyrolysis model is validated with experimental data and is also compared with the predictions of a two-fluid model. An overall improvement is found in light gas, tar, and char yields, among which, the relative error of the predicted tar yield is reduced by more than 50%, while that of the CO and C₂H₄ yields is reduced by approximately 30%. In addition, fluidization behaviors, such as the biomass particle distribution and tar evolution, which cannot easily observed in experiments, are also obtained.

(3)

Under the specific pyrolysis conditions at hand, light gas yield increases considerably from 17% to 25% when the pyrolysis temperature rises from 500 to 700 °C. In addition, the primary tar yield stays around 46% and only shows a slight influence by the pyrolysis temperature. The variation in biomass particle size is found to have a significant impact on particle flow pattern, while pyrolysis products only show small variations within the tested range of particle diameter (<1 mm).