Determination of Zero Dimensional, Apparent Devolatilization Kinetics for Biomass Particles at Suspension Firing Conditions

: As part of the strive for a carbon neutral energy production, biomass combustion has been widely implemented in retroﬁtted coal burners. Modeling aids substantially in prediction of biomass ﬂame behavior and thus in boiler chamber conditions. In this work, a simple model for devolatilization of biomass at conditions relevant for suspension ﬁring is presented. It employs Arrhenius parameters in a single ﬁrst order (SFOR) devolatilization reaction, where the effects of kinetics and heat transfer limitations are lumped together. In this way, a biomass particle can be modeled as a zero dimensional, isothermal particle, facilitating computational ﬂuid dynamic calculations of boiler chambers. The zero dimensional model includes the effects of particle aspect ratio, particle density, maximum gas temperature, and particle radius. It is developed using the multivariate data analysis method, partial least squares regression, and is validated against a more rigorous semi-2D devolatilization model. The model has the capability to predict devolatilization time for conditions in the parameter ranges; radius (39–1569 µ m), density (700–1300 kg/m 3 ), gas temperature (1300–1900 K), aspect ratio (1.01–8). Results show that the particle radius and gas phase temperature have a large inﬂuence on the devolatilization rate, and the aspect ratio has a comparatively smaller effect, which, however, cannot be neglected. The impact of aspect ratio levels off as it increases. The model is suitable for use as stand alone or as a submodel for biomass particle devolatilization in CFD models.


Introduction
An increased awareness of the need for a carbon neutral energy production has led to a corresponding increase in biomass combustion. The world's primary energy consumption in 2019 amounted to 584 EJ of which coal constituted 27% and renewables (including biofuel) to 5%. In the European Union, the numbers are 69 EJ, 11%, and 11%, respectively [1]. Combustion of biomass can be regarded as a carbon neutral source for energy [2]. Coal releases approximately 93 × 10 6 tons CO 2 /EJ, depending on coal type and origin [3,4]. Byutilizing suspension firing technology in existing coal power plants rebuilt to combust biomass particles, a decrease in carbon emissions can be obtained, which enables the European Union's Green deal of no net emission of green house gasses by 2050. In combined heat and power plants, up to 90% of the energy stored in woody biomass can be recovered as heat and electricity [2]. Thus, with this aim, the interest in retrofitting suspension firing units to combust biomass instead of coal has increased over the past decades. Suspension firing of biomass is typically done with small particle sizes (d p = 0.1-3 mm), [5] at high temperatures (T > 1000 K), [6], and at high heating rates (>10 3 -10 5 K/s) [7].
In order to predict biomass flame behavior and boiler chamber conditions, modeling of biomass suspension firing is used [8,9]. A reasonably accurate representation of the devolatilization process is crucial in obtaining correct modeling results. Flame characterization and, consequently, flame attachment to the burner quarl in full sized industrial burners is dependent on the time for onset of devolatilization and the amount of volatiles released [8]. A better flame characterization will allow for higher fuel flexibility and plant efficiency. These challenges can be investigated through CFD combustion modeling.
Modeling a suspension firing unit often involves CFD simulations [10]. To avoid too high computational costs, subprocesses in the particle combustion modeling in CFD require simplifications. Coal particles have historically been modeled as isothermal in CFD combustion simulations, and this approach has sometimes been extended to also include modeling of biomass pyrolysis [11][12][13]. However, model work validated against experimental data [14,15] shows that biomass particles at sizes and temperatures relevant for suspension firing cannot be regarded as thermally thin, i.e., isothermal. Some thermally thick particle conversion models exists in CFD focusing on efficient use of computational power [16,17], which with some limitations adequately predicts biomass devolatilization at fluidized and fixed bed conditions. More complicated particle models are not feasible in CFD for industrial modeling purposes, due to high computational costs [9].
Simplified devolatilization models [9,15,18], which account for the complicated process of biomass devolatilization in a simple lumped fashion, have previously been presented. In these isothermal models, apparent devolatilization kinetic parameters have included heat transfer limitations, which were not otherwise included. None of the models, however, take the effect of biomass density, gas temperature, or particle morphology into account. These are important properties, influencing flame stability and onset of devolatilization.
To compromise between the need for a simple devolatilization model and the need for describing the complicated phenomenon of biomass particle pyrolysis, this paper introduces lumped Arrhenius kinetic parameters for a single first order global pyrolysis model. The parameters of a zero dimensional, isothermal model are found by fitting to predictions of a semi-two dimensional model [19] that previously has been validated against experimental data. The comparison is done for different gas temperatures, particle sizes, particle aspect ratios, and densities. The combined effect of these experimental properties on the Arrhenius parameters is quantified using the multivariate data analysis method, partial least squares regression (PLS). To the knowledge of the authors, no experimental data valid for suspension firing of woody biomass exist, which has not been used in the development of the semi-two dimensional model. Thus, the present model is not directly evaluated further against experimental data.

Method and Model Description
The model development process is done in two parallel tracks, which are subsequently compared and used for final model determination. An overview can be seen in Figure 1. In a previously published paper, a rigorous semi-2D model was developed [19]. In this paper, a 0D model is set up. Subsequently the lumped 0D Arrhenius parameters, A and E, which provide the best agreement between the 0D and the semi-2D model for a range of selected cases are determined. Then, a PLS model is established, which calculates the lumped Arrhenius parameters from particle radius, aspect ratio, particle density, and maximum gas temperature. The model development process is described in more detail in the subsequent sections.

The Semi-Two Dimensional Model
The semi-two dimensional model is described in detail elsewhere [19][20][21], and only a short introduction will be given here. The model was first presented for devolatilization at medium heating rates by Thunman et al. [20], and subsequently developed further by Ström and Thunman [21]. It was recently adapted to account for biomass devolatilized under suspension firing conditions by Leth-Espensen et al. [19]. The model is a shell model capable of describing both cylindrical and spherical geometries. The cylinder is the preferred simple geometry to model biomass particles, [22] whereas coal particles tend to be almost spherical in nature. The cylindrical model takes changes during devolatilization in both radius and length into account using one variable only. The decrease in particle length is assumed comparable to the decrease in radius, thus particle shrinking is described only by the radius. Overall, the cylindrical model accounts for the internal temperature gradient, the volume, and surface area in two dimensions, whereas the heat transfer is solved in a one dimensional manner. The spherical model is one dimensional. Regardless of the shape, in the model, the particle is divided into three concentric shells. The innermost layer consists of moist biomass, the middle layer is dry biomass, and the outer layer is char. The shells move inwards during the devolatilization, transforming the model particle from consisting of practically only a moist shell in the beginning to be all char after full devolatilization. The model accounts for intraparticle heat and mass transfer.

The Zero Dimensional Model
The 0D model is based on the following assumptions: • The particle is isothermal. • The particle has constant density (shrinking particle). • The particle is dry. • The particle is spherical. • The initial particle diameter, R, is defined as the initial cylinder diameter, according to recommendations in [19]. • The devolatilization enthalpy is assumed to be 0. • Both the influences of kinetics and heat transfer are described by a single first order reaction model. • The wall radiation temperature is defined as T w = T g − 200 K. The devolatilization process in the isothermal particle is modeled as a global reaction, using the single first order reaction (SFOR) model, described in Equation (1) [23].
Here, t is the time, Υ is the fraction of volatiles released, Υ * is the fraction of volatiles present in the particle at t = 0, and k is an Arrhenius reaction rate constant given in Equation (2).
R g is the gas constant, T is the particle temperature, and A and E are the Arrhenius preexponential factor and activation energy, respectively. The temperature in the particle is uniform and determined by the radiation and convective heat transfer and is given in Here, ρ is the density, C p is the specific heat capacity of the fuel, R is the particle radius, σ is the Stefan-Boltzmann constant, is the emissivity, T w is the radiation temperature (reactor wall temperature), and h is the heat transfer coefficient. Expressions for the model input parameters are given in Table 1. The two coupled differential equations are solved using the ode45 solver in Matlab ® . ) 2 + ... [27] ...1.369784 × T g 1000

Fitting the Arrhenius Parameters
In order to approximate the semi-2D model with a 0D model, the Arrhenius equation must account both for the rate of the kinetics and for any heat transfer limitations in a given particle. The apparent kinetic parameters are obtained by fitting the Arrhenius equation. The pre-exponential factor, A, and the activation energy, E, in the Arrhenius equation are coupled, though, so the procedure suggested by Rawlings and Ekerdt [30] is used here. In order to minimize the correlation between the fitted parameters, a modified Arrhenius equation, as seen in Equation (4), has been used.
k(T re f ) is the rate constant at a reference temperature, here 1600 K, the midpoint in the temperature interval. k(T re f ) and E are fitted to the result from the semi-2D devolatilization model using the lsqcurvefit command in Matlab ® , which works by minimizing the residual sum of squares between the model results from the semi-2D model and the 0D model. The residual sum of squares is given in Equation (5). The value of A can be calculated from k(T re f ) and E.
An example of a fitted curve and the semi-2D cylindrical model output can be seen in Figure 2.

Chemometrics
Chemometrics is a statistical approach to extract data from chemical or biological data sets. A common method within chemometrics is partial least squares regression (PLS) [31,32]. In PLS, a correlation between a matrix of input variables (X) and an output variable (Y) is determined. PLS is a linear regression method. An in depth description of PLS is beyond the scope of this manuscript, but can be found elsewhere [31][32][33][34]. The PLS models presented here are calculated in PLS Toolbox version 8.1.1 and Matlab version 9.3.0 (R2017b).

Parameter Definition
The input parameters tested for the calculation of A and E in the 0D model are particle radius, particle density, gas temperature, and particle aspect ratio, as they have previously [10,24,35] been shown to influence the devolatilization process. The applied parameter spans for each variable can be seen in Table 2. They cover the parameter variations relevant for biomass suspension firing conditions. A total of 35 simulations with the semi-2D model have been made to span the parameter space. The applied parameter values for each simulation can be seen in the Appendix A. To show the correlation (negative or positive) between input parameters, the correlation coefficient chart is made. It is shown in Figure 3. The higher the absolute value in the correlation coefficient matrix, the higher degree of correlation between two input parameters. When the degree of correlation is high, the effect of each input parameter cannot be separated due to confounding. Here the degree of correlation is numerically low (≤0.23), which means the influences on the model results can largely be ascribed to each individual input parameter.

Preprocessing
Preprocessing is an important step in PLS, and the individual preprocessing methods tested can be seen in Tables 3 and 4. When studying, e.g., the effect of the particle radius on the Arrhenius parameters, it is clear that the correlation is nonlinear as seen in Figure 4. Thus, both R, log(R), R 1/2 , and R 1/3 have been tested, and the quality of the preprocessing is then determined based on the explained variance for the input matrix (ExpVar(X)), the output parameter (ExpVar(Y)), and the root mean squared error of cross validation (RMSECV) for both A and E, respectively. Similar considerations have been made for determining the preprocessing of ρ, T g , and AR. In these cases, the nonlinearity of the correlations between input and output variables are so weakly pronounced that no significant improvement is added to the model by choosing various more complex preprocessing methods and thereby increase model complexity beyond what is necessary.  Furthermore, the input variables here are not on the same scale, so to account for this, all input parameters have also been scaled to account for unit variance.

Cross Validation
The cross validation is made to ensure that the presented model is robust. It can be done in many ways; here, the random subset method [36] is used with 6 splits and 20 iterations, thus on average 17% of the data set is removed in each iteration. All models have been made with two PLS components. When doing the cross validation, preferably the obtained RMSECV values should be low, and ExpVar values high. The values in Tables 3 and 4 are the basis for choosing the relevant preprocessings for a model. Due to the lower RMSECV(A) and RMSECV(E) and higher ExpVar(X) (compared to the number of input data in X) and ExpVar(Y) values, model 2 has been chosen both for A and E. If models appear equally good or almost equally good based on RMSECV and ExpVar data, the simpler model is usually preferred within PLS.

Arrhenius Plots
The Arrhenius parameters for four different particle sizes (d p = 78.8 µm, d p = 800 µm, d p = 1.51 mm, and d p = 3.12 mm) found by fitting the 0D model to the semi-2D model for different aspect ratios have been used to generate the Arrhenius plots in Figure 5. They have been determined for T = 1600 K, and plotted with these A and E values for a short temperature interval to illustrate the effects of the different parameters. It can be seen that the aspect ratio only has a minor influence on the Arrhenius parameters, and that the effect of aspect ratio levels off for higher values of AR. The effect of aspect ratio on the apparent rate constant is more pronounced for the smaller particles.
Furthermore, the change in apparent rate constant as a function of temperature is smaller for the large particles. It is due to the increased heat transfer limitations in the large particles. This is in contrast to the smaller particles, which are kinetically controlled, and where the temperature, consequently, has more influence over the apparent rate constant. (d) d p = 3.12 mm. Figure 5. Arrhenius plots for the fitted values for A and E a for different aspect ratios for four particle sizes, ρ = 700 kg/m 3 , T g = 1600 K. The Arrhenius parameters have been used to generate Arrhenius plots for a smaller temperature interval to illustrate the effects of the input variables.

PLS Model
The PLS model presented here is developed to be able to predict the Arrhenius parameters in a lumped SFOR kinetic scheme accounting both for the reaction rate and the heat transfer limitations. A and E values have been found for multiple values spanning the parameter ranges (see Appendix A for all parameters combinations), by fitting the 0D model to the semi-2D model. A PLS model has then been generated using these new A and E values. Two graphs showing the 0D model predictions of log(A) and E as a function of the log(A) and E values given by fitting the 0D model to the semi-2D model (cf. Section 2.3) can be seen in Figure 6. The figure shows that there is good agreement between the values for log(A) and E predicted by the simple 0D model and by fitting to the semi-2D model. The Arrhenius parameters for the 0D model can be found through the regression vector. The regression vectors both for the preprocessed input parameters and for the raw input parameters can be seen in Table 5. The regression vectors for the preprocessed input parameters show that both for log(A) and E an increase in AR, T g , or log(R) would result in a decrease of log(A) and E. In other words, the PLS model analysis shows that if the gas temperature increases the apparent activation energy decreases. This is as expected, because the model covers both heat transfer limitations and kinetics. An increase in gas temperature means a higher absolute heat transfer as well as a higher reaction rate, resulting in a smaller apparent activation energy. The raw regression vector can be used for predictive purposes. In Equations (6) and (7), the values for log(A) and E can be estimated for pyrolysis situations similar to the ones used for model development.  Here, A is the preexponential factor in s −1 , E is the activation energy in J/mol, R is the particle radius in µm, ρ is the density in kg/m 3 , AR is the aspect ratio, and T g is the maximum gas temperature in Kelvin.

Evaluation of Particle Conversion Predicted by the 0D Model and Semi-2D Model
Mass loss profile obtained from experiments is desirable for validating the proposed kinetic model. However, time-resolved mass loss data during devolatilization of biomass at high heating rates and peak temperatures is scarce in open literature. Although experiments using a flat flame burner can reach the desired conditions, such as the studies conducted by Johansen et al. [15] as well as Lewis and Fletcher [37], mass loss data at the initial conversion stage are difficult to obtain and have not been reported to the best knowledge of the authors. Therefore, validation of the 0D model (Equations (1)-(3) and Equations (6) and (7)) is here done by comparison of predictions to the results from the semi-2D model presented by Leth-Espensen et al. [19]. The semi-2D model was validated using a range of different studies. Agreement between the 0D and the semi-2D model has been checked for all input parameter combinations given in the Appendix A. Examples of conversion comparisons obtained by the semi-2D and the 0D model for the different particle sizes can be seen in Figure 7. The 0D model shows good agreement with the more complicated semi-2D model and predicts devolatilization times well. The model predictions are best for smaller particles as they are heated up rapidly, and thus closer to being isothermal, but it should be noted that the conversion time is relatively accurate for all particle sizes.    (6) and (7) compared to a semi-2D devolatilization model [19]. ρ = 700 kg/m 3 , T g = 1600 K.

Comparison of Devolatilization Times Predicted by the 0D and Semi-2D Model
An overall comparison of the devolatilization time predicted by the semi-2D model and the 0D model can be seen in Figure 8. The devolatilization time is defined as the time, where 95% of full devolatilization is reached. The 0D model predicts most devolatilization times adequately, showing limitations only for particles with a large radius at a high gas temperature. In these cases the 0D model may predict negative activation energies, as can be deduced from Equation (7), see examples in the table in Appendix A. As negative E values entail that the reaction rate decreases for increasing temperatures, the model cannot be trusted in these instances and such values have been disregarded in Figure 8. From Figure 8, it can be seen that the 0D model tends to overpredict the devolatilization time for small particles and underpredict the devolatilization time for large particles, however, the prediction capability is overall good, and the trend is captured well by the model.

Comparison of the Arrhenius Parameters to Literature
The lumped Arrhenius parameters, which can be calculated from Equations (6) and (7) covers most of the area for Arrhenius parameters previously covered in literature [10,38,39] for simple first order reactions relevant for suspension firing. In Figure 9 a comparison can be seen. Previously reported SFOR Arrhenius parameters are often only determined short temperature span, but extrapolated outside this region, the model presented here has gas phase temperature as a parameter, thus allowing for a broader temperature interval, where the model is accurate. Consequently, it allows for better modeling of industrial sized suspension firing boilers. Furthermore, as the model also takes particle size into account one equation can describe all sizes relevant for suspension firing. The model's tendency to predict slightly too low reaction rates for small particles at high temperatures, also described in Section 4.2, can be seen in Figure 9 by comparing to the Arrhenius parameters for small particles from Johansen et al. [10] as well. The model's slightly higher predicted devolatilization times for small particles is the result of providing a model, which account for broad parameter ranges within particle size, density, gas temperature, and aspect ratio. Especially the particle size is an important factor in determining lumped kinetics for biomass devolatilization and covering almost two orders of magnitude results in less accuracy in the interval ends.  [10]. This model accounting for both kinetics and heat transfer is based on experimental data from two single particle combustors (HR = 10 2 -10 3 , T g = 1480-1831 K, d p = 1790-5800 µm), (HR = 30, T g = 1050-1267 K, d p = 10 900 µm) and an entrained flow reactor (HR = 10 5 , T g = 1405-1667 K, d p = 63-90 µm), and kinetic data from Wagenaar et al. [38] obtained in an entrained flow reactor (HR = 2-12 K/s, T g = 773-873 K, d p = 100-212 µm), and parameters from Ontyd et al. [39] obtained in a drop tube reactor (HR 3.7 × 10 3 -7.8 × 10 3 K/s, T g = 950-1300 K, d p = 100 µm). The reaction rate predicted by the model depends on particle size, density, and aspect ratio. Here, ln(k) is shown for AR = 1 for comparative reasons, as the literature data are reported for spherical particles. The model's ability to account for density and particle size changes means it covers a wider area, here shown in dark gray. Extrapolations of both the model presented here and literature models are shown in light gray.

Conclusions
A 0D model has been developed to describe devolatilization of cylindrical biomass particles using a lumped single first order Arrhenius equation to account for both kinetics and heat transfer limitations in the particle. The model accounts for aspect ratio, particle density, gas temperature, and particle radius, because of their influence on devolatilization rate and on flame shape and stability. Including these parameters in the devolatilization model allows for better CFD simulations of industrial sized units. Calculations show that of these parameters, especially particle radius and temperature, are important for the prediction of the degree of particle devolatilization. Model predictions are compared to those of a more complicated semi-2D model, showing good agreement, especially for the smaller particles, where the heat transfer limitation is less severe. The presented model also shows that it is important to include differences in aspect ratio, but that the influence of aspect ratio on devolatilization rate levels off as aspect ratio increases and consequently end effects become less important.
The devolatilization model presented is simple and can be implemented into CFD without adding substantially to the computational costs. By implementation into CFD simulations, as, e.g., the ones described by Johansen et al. [8], it is likely that an even greater match between numerical simulations and the physical phenomena during flame development in a suspension firing unit can be obtained. Future work should include making the CFD simulation described above, and expanding the model to include different kinds of biomass and municipal waste.

Data Availability Statement:
The data presented in this study are available in the Appendix A.

Acknowledgments:
The authors gratefully acknowledge the financial and advisory support received from Ørsted A/S, Burmeister, and Wain Scandinavian Contractors A/S, and Rambøll A/S. We also thank the Nordic Five Tech (N5T) alliance for financial support.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: