Developing an all-round combustion kinetics model for nonspherical waste-derived solid fuels

The utilization of challenging solid fuels in the energy industry (especially the ones derived from wastes) has a big priority nowadays, as it is a valid option to keep the recent EU directive related to the decrease of landfills. However, there are serious technical challenges, connecting to the lack of knowledge about the behavior of these fuels in the combustion chamber. This paper discusses the specific aspects of developing particle models concerning the combustion of these non-conventional fuels. A new modeling approach is presented, using which it is possible to develop an all-round particle model that includes every significant influencing process. Moreover, it does not have any restrictions regarding the shape, size and the origin of the particle. As an integral component of this model, the distinctive aspects of intrinsic reaction kinetics related to waste fuels are presented as well.


Introduction
Combustion of municipal solid wastes is a feasible process, which follows the most recent EU directive (Circular Economy Package 2018) that aims to lower the rate of landfill drastically. The main benefit of combustion is that it reduces the waste's volume (and the area needed for the landfill as well), while energy is produced during 1 3 the process. However, it is essential to highlight that only the non-recyclable components should be utilized in this way, as recycling is more reasonable environmentally and economically as well.
In most cases, waste combustion is performed in older boilers designed for conventional fuels such as coals or simple biomasses, commonly co-combusted with those (Gómez- Barea and Leckner 2010). This kind of complex operation requires careful design choices to prevent any possible malfunctions, which are impossible without indepth knowledge of the processes present in the combustion chamber. Designing a boiler often includes modeling these, which is a computation-heavy task, if every process is included. One of the most important aspects is the behavior of the particle. A detailed particle combustion model is affordable in case of a more straightforward (and computationally cheap) 1D boiler model. However, the complex 3D CFD models have so many submodels that simplifications are required (Niemi and Kallio 2018).
In a particle model, it is an important task to cover the combustion reactions, which is usually done based on the classical Arrhenius equation. In its original form, it assumes an ideal sample and environment, without any limiting effect. It is a reasonable assumption in the case of thermogravimetric (TG) measurements, but modeling an industrial combustion chamber needs a more sophisticated approach (Table 1). These expanded models are usually referred to as apparent. In contrast, the ones concerning only pure reaction kinetics are called intrinsic. Also, the latter provides the theoretically highest reaction rate, only achievable in TG conditions (Di Blasi 2009;Gómez-Barea and Leckner 2010;Khodaei et al. 2015;Jiang et al. 2017). The difference between these two approaches generally increases with the size of the sample [measured by Liu et al. (2016)].
The intrinsic kinetic submodels rely on model-free and model-fitting methods; multiple papers collect the commonly used ones (Çepelioğullar et al. 2016;Radojevic et al. 2018). The proper way should be selected based on the purpose of the results, as there are quite strict limits regarding the compatibility of the parameters of different reaction models. Because of this, utilizing kinetic data available in the literature is not always possible, considering the kinetic model used during the initial evaluation is essential (Várhegyi et al. 2018).
Several classical models are available concerning the effect of the particle size. Field (1969) developed a method by measuring a k s surface reaction rate for various coal chars, and used a modified Arrhenius equation to describe the shrinking of the spherical particle. His results [and the supplementary measurements published by Smith (1971)] are still widely used for coal chars. The limit of this method is the actual surface area S(m(t)) , because formulating its change during the conversion is possible only for spheres and some other simple shapes.
More complex particle models were developed throughout the years; however, most are only capable of describing spherical samples and still failing in the case of geometrically challenging ones. The most common classical sphere models are the uniform conversion model (UCM) and the shrinking unreacted core model (SUPM) (Gómez-Barea and Leckner 2010). The first assumes constant diameter, homogenous reaction rate, and decreasing density during the conversion, which are the conditions of TG measurements. The SUPM, however, assumes constant density and decreasing diameter, and the reactions happen only on the sample surface. This approach is close to the conditions in a real boiler (Gómez-Barea and Leckner 2010). Modeling coal chars that have a porous structure is usually performed using Heterogenous, because the heat conductivity of the particle slows down the heat transfer Reaction rate in the particle Homogenous, the conversion is the same in the whole sample Changing with the local temperature Oxygen availability The optimal amount is artificially maintained; available in the whole sample, because of the small particle size Often sub-optimal amount; available only on the surface Devolatilization controlling process Only the reaction kinetics, the temperature of the sample changes almost instantly The particle heat-up is slow and heterogenous, which limits the reaction rate Char combustion controlling process Only the reaction kinetics, because of the precisely maintained atmosphere and the small particle size The oxygen transfer is much slower then the reaction; the reaction front is narrowed to the particle surface the random pore model (RPM) method developed by Bhatia and Perlmutter (1980). Several research groups have used this method successfully (Kreitzberg et al. 2020;Bibrzycki et al. 2016;Beckmann et al. 2017). Several particle models cover the combustion of nonconventional chars. It is a possible method to rely on measuring the change of the particle structure and surface, based on which the Arrhenius equation can be modified to describe the apparent conversion (Liu et al. 2016). It is also a common approach to develop a detailed CFD model for the particle (Beckmann et al. 2017;Xue et al. 2019). The benefit of this method is that it can easily include the limiting effects of the environment. However, this kind of modeling is very computation heavy, and the result is usually too complex to use it as a part of a bigger boiler model. It is possible to solve this problem by developing a simplified submodel, which includes the relevant part of the results, but in a less computation-heavy way. Niemi and Kallio (2018) presented a method based on neural networks, which, after the time-consuming initialization, results in a very fast algorithm. However, these models still consider sphere-like particles or other simple geometries, and the devolatilization is also neglected.
Many papers discuss the devolatilization of solid fuels separately from the char combustion. Sunflower seed shells were investigated by López et al. (2014) and solid recovered fuels (SRF) by Conesa and Rey (2015), Çepelioğullar et al. (2016), and Radojevic et al. (2018); both will be covered in this paper later. However, these studies usually neglect the effect of the particle size and shape, which was summarized (but not modeled) by Gómez- Barea and Leckner (2010).
It is clear that there is no modeling method available in the literature that is solely capable of covering both the devolatilization and the char combustion of a nonspherical particle in a computation-efficient way.
This paper aims to connect the intrinsic (Szűcs et al. 2020) and apparent (Szűcs and Szentannai 2019) combustion models published by the same authors. The presented modeling method is capable of describing the whole apparent conversion (devolatilization and char combustion) of any fuel particle, regardless of its shape, size, and origin. In this work, the previous results are summarized, some important details are refined, and more insight is presented about the implementation. The method uses a new conversion function to modify the Arrhenius equation, which can be determined by modeling all aspects of the apparent combustion. The complexity of the output function can be changed according to how it will be utilized later. As an example, the intrinsic kinetic evaluation for an SRF and the whole apparent modeling procedure of a sunflower seed shell pellet are presented.

Intrinsic kinetics
The Arrhenius reaction equation is used for the intrinsic kinetics model (Eq. 1), which describes the reaction rate as the function of the temperature and the conversion (Green 2008). The pressure can affect the reaction rate as well, but it is not that significant in industrial boilers, as those operate close to the ambient pressure. Because of this, the pressure is beyond the scope of this paper.
Here, x is the conversion, which is the ratio of the actual and final reacted mass (Eq. 2).
The combustion of real fuels contains thousands of subreactions, which, in most cases, are impossible to describe, so it is common to collect them in reaction groups. These groups can be modeled using one Arrhenius equation and can be associated with the main fuel components, such as cellulose, hemicellulose, and lignin for biomasses, or various types of plastics. These reaction equations may be summarized using a c i mass share to describe the complex fuels (Eq. 3) (Conesa and Rey 2015).
Choosing the correct reaction function f (x) is an integral part of the evaluation. The most common is to consider n-thorder reactions, but others are available as well (Aboulkas and Harfi 2008).
The present authors published an in-depth intrinsic kinetic study (Szűcs et al. 2020), in which the validity of four different reaction functions to describe the combustion of an SRF sample were compared. The least-squares method was used to find the kinetics parameters (Várhegyi et al. 2011). The reaction model contains three pseudo reaction groups, two for the devolatilization of the two main components (biomass and plastics) and one for the combustion of the remaining char.
The least-squares method uses Eq. (4) to compare N i number of x m (t) measured conversion graphs containing N j points to an x c (t) calculated one, which results in an F fitness value. A genetic algorithm (McCall 2005) generates randomized kinetic parameters and calculates F for each of them. After numerous iterations, it concludes in the best parameters, which produce the smallest F.
(2) Table 2 summarizes the four reaction models that were used during the evaluation. The most complex is the distributed activation energy model (DAEM) developed by Anthony et al. (1975) and successfully used multiple times for complex solid fuels (Várhegyi et al. 2011;Cai et al. 2013bCai et al. , 2014Wu et al. 2014;Zhang et al. 2015;Lin et al. 2018). It assumes that every reaction groups consist of an infinite number of sub-reactions, and their activation energies have a specific distribution, which is described by a D(E) distribution function.
A sensitivity analysis was also performed based on the work of Cai et al. (2013a) to determine the impact of the various kinetic parameters on the fitness values. During the analysis, every parameter was lowered and increased by a maximum of 50% of the optimized value. The details of the implementation can be found in the original work (Szűcs et al. 2020).

Apparent kinetics
The apparent kinetic model includes the reaction groups determined during the intrinsic evaluation, which have two categories with their specific modeling technique. These are the decomposition-like devolatilization reactions and the char combustion reactions that are handled as surface reactions. In this model, they both contain a function of the apparent conversion, the mass-related reaction effectiveness factor ( m (x app ) ), which is an expanded version of the methodology developed by Gómez-Barea and Leckner (2010). Essentially, this function is a tool to scale up the intrinsic kinetics for a big particle. It shows how high the apparent reaction rate is compared to the intrinsic one, which is assumed to be the theoretical maximum.
The devolatilization reaction equations (first term on the right side of Eq. 9) are similar to Eq. 1, because they are independent of the oxygen transport. The purpose of m (x app ) here is to model the effects slowing the release of the volatiles. For example, before its fragmentation, the solid structure of the particle can act as an obstacle.
The char combustion reactions (second term on the right side of Eq. 9) are surface reactions, and the availability of the oxygen limits their rate. The oxygen concentration is also included as a linear scaling factor. The description is analog to the parallel resistances, as shows Eq. (10) in a simplified form, where k s and k d represent the two competing processes; whichever is the slower will control the whole conversion.
To understand m (x app ) , consider Eq. (11) as the surfacerelated reaction equation, which includes two parameters, the surface-related reaction rate constant and the changing surface of the particle. In the case of nonspherical particles, both of those are hard to measure.
To connect k s to the intrinsic kinetics (which are easier to determine), the surface dependency needs to be redeemed. Gómez-Barea and Leckner (2010) presented a general method by introducing the S m (x app ) mass-related reaction surface as the function of the apparent conversion (Eq. 12). Its purpose is to convert k s to a mass-related, but still apparent reaction rate constant k m (x app ) . By introducing the m (x app ) mass-related reaction effectiveness factor, and assuming, that S m0 is a reference S m (x app ) related to the pure TG measurement conditions, the connection can be described without the reaction surface.
This methodology is the key component of Eq. (9), as it makes it possible to summarize all of the reaction groups by (12) k m x app = k s S m x app = k s S m0 m x app = k r m x app .
Expanded n-th order (Várhegyi et al. 1996) dx c (t) Distributed activation energy model (DAEM) (Anthony et al. 1975) converting them to the same mass-related, but still apparent form.
In this work, Eq. (9) is used to model the apparent conversion of a sunflower seed shell pellet based on our previous work (Szűcs and Szentannai 2019). To calculate the apparent conversion (left side of Eq. 9), m (x app ) needs to be determined. First, an intrinsic reaction model (using the method presented in the previous section) was developed, from which determining m (x app ) is possible following the next steps, each aimed to derive one of the parameters in Eq. (9): 1. An intrinsic kinetics evaluation resulting in k r (T, x). 2. Measuring the apparent conversion graph ( x app (t) ) of the fuel sample during macroscopic combustion experiments. 3. Developing an appropriate model of the environment of the macroscopic measurements, which provides the temperature profile and the oxygen mass transport around the sample ( T(t) and k d (x)).

Finally, only m (x app ) remains an unknown variable in
Eq. (9), so it can be determined.
The intrinsic model of the pellet consists of two decomposition-like reactions ( n 1 = 2 ), one for the drying and one for the devolatilization, and one surface reaction for the char combustion ( n 2 = 1 ); more details are in the original work (Szűcs and Szentannai 2019).
After performing the previous procedure, by changing the environmental properties (step 3 in the previous list), the new Eq. (9) can model the sample's apparent conversion in the new combustion environment.
The parameter k d (x) is unique, as usually, the oxygen mass transfer does not depend on the conversion, but to add it to the chemical reaction rate constant, it is necessary to convert them to the same dimension.
To determine the environmental parameters, a CFD model was implemented in Comsol 5.3. The model consists of two main parts, one for the devolatilization, where the primary process is the heat transfer, and one for the char combustion focusing on the mass transfer. During setting up the model, the whole geometry was cut in half by a symmetry plane to reduce the computation cost.
The samples' devolatilization is assumed to be controlled by the heat-up because of their size, so the heat transfer model is essential to determine the apparent conversion. The interior of the particles was divided into multiple parts (Fig. 1, which differs from the mesh of the numerical solution), all of them having their own temperature graphs. A transient heat transfer model was used, where the initial temperature of the samples was 20 °C, while the interior of the oven was 850 °C. The model includes convection, and the corresponding velocity field relating to the exiting volatiles was calculated by a stationer fluid dynamics model using the k − turbulence approach. The radiation model used surface radiation combined with the discrete ordinates method. The internal heat conduction follows the Fourier law. The primary heat source is the flame associated with the instant ignition of the volatiles. It was modeled as a volume with a constant heat flow calculated from the stoichiometry of the fuel. The char combustion model includes three sub-models. A stationary heat transfer model calculates the temperature field, which consists of no convection, so the primary driving process is the temperature difference between the interior of the oven and the environment. A fluid dynamics model calculates the velocity field, in which the only driving force was the buoyancy effect; its input is the density profile derived from the temperature. The third sub-model is the species transport model, which contains diffusion and convection (calculated from the previous velocity field). Zero oxygen concentration was set as a boundary condition on the particles' surface, which assumes that the oxygen availability limits the reaction rate. It results in the oxygen flux around the particles from which k d (x) can be calculated based on the surface average of the oxygen flux.
Using all this information, m x app can be derived from Eq. (9). The first task is choosing the correct function form, which should be as simple as possible, especially if computation time is a severe issue. However, in this study, the aim was to introduce a general function (Eq. 13), which is not the simplest, but it is a good starting point for all kinds of solid fuels and can easily be simplified when needed. This function resembles the Arrhenius equation, which promotes a similar physical meaning as well. E η can be understood as a conversion limit, at which the apparent reaction rate starts to increase until it reaches the intrinsic reaction rate. Parallel to this, m (x app ) also increases from its initial condition m (0) (which cannot be 0) to 1. A η represents a scaling parameter, similarly to the preexponential factor.  13) is possible by the same optimum seeking algorithm used in the intrinsic kinetic evaluation. It compares the apparent conversion measured in the oven and the one calculated from Eq. (9). The details of the implementation and the final parameters are available in the original work of Szűcs and Szentannai (2019).
After finishing the calculation, it is possible to validate the assumption of zero oxygen boundary condition on the particle surface. For this, the reaction rates corresponding to the chemical reaction kinetics of the char combustion ( v r ) and the oxygen transfer ( v d ) should be compared. Both were calculated by assuming the other to be irrelevant (Eqs. 14 and 15, accordingly).

Experimental
The intrinsic kinetic evaluation of the SRF sample was performed using the results of TG measurements in air atmosphere at 5, 10, and 15 °C/min heating rates following the guidelines of Várhegyi et al. (2011Várhegyi et al. ( , 2018. The sample size was around 2 mg. The specifics of the TG apparatus are described by Bakos et al. (2018). The sample was an SRF processed by a Hungarian junkyard (its composition is in Table 3), and no other preparation was performed before the TG measurements, except grinding and homogenizing. This made it possible to maintain the sample's original state, e.g., the components were not separated, so the model covers the possible interactions as well. More details about the sample and the equipment are available in the work of Szűcs et al. (2020).
The sunflower seed shell pellets used in the apparent kinetic evaluation were 20 ± 5 mm long cylinders with a diameter of 6 mm, and their mass was approximately 1 g. The sample's composition is shown in Table 4. Similar to the SRF sample, before the TG measurements, the only preparations were grinding and homogenization. Due to simplifying the method, only one heating rate was used (10 °C/min). However, it is important to emphasize that in sensitive projects, it is essential to apply multiple heating rates to get accurate kinetic data. The same TG apparatus was used as before (Bakos et al. 2018).
The macroscopic combustion experiments were performed in a programmable oven described in more detail by Szentannai et al. (2015). In contrast to the TG measurements, the particles used in the macroscopic combustion experiments were not ground. Fourteen whole pellets were placed on a spoon-like rod and inserted into the preheated oven, where their mass change was recorded. More details are available elsewhere (Szűcs and Szentannai 2019). Table 5 show the best possible fitness values achieved during the intrinsic evaluation. The DAEM method provided the smallest F , but it is essential to mention that it required more time to finish the run (6.5 h vs. 5 min for the others on a regular PC). It means that if the reaction model is part of a bigger, already computation-heavy model, it is more beneficial to use a slightly less precise (and less time-consuming) reaction model. The results, including the kinetic parameters, can be found in more detail in Szűcs et al. (2020).

Figure 2 and
During the sensitivity analysis, it was found that the activation energies have the most significant impact on every reaction model's fitness value, followed by the preexponential factors. The parameters introduced by the more complex models have a less important role, which indicates that the increased precision is not connecting to the value of the parameters, but rather to the changed model structure itself. The detailed results are in the original work (Szűcs et al. 2020).  Figure 3 shows the temperature profile of the CFD model developed for the devolatilization. It indicates the temperature difference inside the particle, which is also visible in the temperature graphs of the cells (Fig. 4). The slowest heat-up (connecting to the interior cells) needed over 30 s, while the outer cells reached equilibrium in less than 10 s. Figure 5 shows the oxygen concentration calculated by the species transport model. The convective term mainly determines the overall profile of the oxygen distribution. However, the direct oxygen transport to the particles is covered by diffusion. It means that these processes are equally important; it is not recommended to neglect either of them in a proper model.
After combining all of the results, the graph in Fig. 6 was found to be the optimal m (x app ) , which provides a close resemblance to the measured conversion (Fig. 7). The value of m (x app ) changes from the initial condition to 1 between the conversion values 0.6 and 0.85. The initial low function value means that before this transition, the apparent reaction rate is much slower than the intrinsic kinetics. This changes at the end of the devolatilization, from where the reaction rate increases up to the intrinsic one. However, it does not translate to rapid char combustion, as the oxygen availability limits it. The narrow transition region suggests that one significant structural change   (Várhegyi et al. 1996) 2.28 × 10 −4 Distributed activation energy model (DAEM) (Anthony et al. 1975) 2.07 × 10 −4 Fig. 3 The temperature profile (K) of the devolatilization (early stage, 10 s after starting the measurement) 1 3 (e.g., the collapse of the pores) happens, which makes the limiting effect of the particle's shape and size irrelevant. Finally, plotting v r and v d reveals the relation between the controlling effect of the reaction kinetics and the oxygen availability for the char combustion (Fig. 8). The diffusion reaction rate is significantly lower than the chemical one, and the difference is a magnitude of 10 7 , which means the reaction consumes every available oxygen on the surface. It proves that it was correct to assume zero boundary condition in the char combustion model.

Conclusion
In this paper, a modeling procedure was introduced, which can describe the effect limiting the reaction rate originated from the continuously changing shape of big, nonspherical fuel particles. The method has two parts, the first is the determination of the intrinsic kinetics, and the second is to bridge it to the apparent kinetics of the real macroscopic particles. The intrinsic kinetics have a significant role in the model, so a detailed evaluation process based on the leastsquares method was presented to show the unique aspects related to waste-derived fuels. The most complex DAEM reaction model provides the most precise result; however, its computation demand can make it a less favorable choice.
The connection between the intrinsic and the apparent kinetics was developed as the mass-related reaction effectiveness factor m (x app ) , which represents a certain reaction depth, in which the reaction occurs. It also means the ratio of the actual reaction rate and the theoretical maximum associated with the conditions of a TG measurement. The applied method requires determining all parameters of a reaction equation expanded with the limiting processes present in a regular combustion chamber. For this, macroscopic combustion measurements and two CFD models were used. An Arrhenius-like differential equation is a proper initial form for m (x app ) , albeit it needs simplifying if used as a part of a computation-heavy model. Acknowledgement This work was performed in the frame of the FIEK_16-1-2016-0007 project, implemented with the support provided from the National Research, Development and Innovation Fund of Hungary, financed under the FIEK_16 funding scheme.
Funding Open access funding provided by Budapest University of Technology and Economics.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/. Time, s r , s −1 , Fig. 8 The reaction rates corresponding to the reaction kinetics and the oxygen availability of char combustion