Validation of MATLAB algorithm to implement a two-step parallel pyrolysis model for the prediction of maximum %char yield

Numerical modeling of biomass pyrolysis is becoming a cost and time-saving alternative for experimental investigations, also to predict the yield of the by-products of the entire process. In the present study, a two-step parallel kinetic model was used to predict char yield under isothermal condition. MATLAB ODE45 function codes were employed to solve a set of differential equations that predicts the %char at varying residence times and temperatures. The code shows how the various kinetic parameters and mass of pyrolysis products were determined. Nevertheless, the algorithm used for the prediction was validated with experimental data and results from past works. At 673.15 K, the numerical simulation using ODE45 function gives a char yield of 27.84%. From 573.15 K to 673.15 K, char yield ranges from 31.7 to 33.72% to 27.84% while experimental yield decreases from 44 to 22%. Hence, the error between algorithm prediction and experimental data from literature is − 0.26 and 0.22. Again, comparing the result of the present work with the analytical method from the literature showed a good agreement.


Introduction
Numerical simulations provide a more detailed insight in various aspects of pyrolysis processes under varying conditions. Pyrolysis being the first stage in any thermal treatment of biomass [1]. Studies have shown that its kinetics and products have a significant influence on the process [2]. However, char produced during pyrolysis of biomass not only reduces the amount of carbon emitted into the atmosphere, but it is also an eco-friendly alternative for activated carbon and other carbon materials [3]. Many models have no closed form solutions and hence it is a need to seek approximate solutions by means of numerical methods [4]. According to Atkinson et al. [5], it is possible to find closed form solutions for simple differential equations. An ordinary differential equation (ODE) is an equation that contains an independent variable, a dependent variable, and derivatives of the dependent variable. For some reasons, MATLAB does not include Euler functions. However, it has very sophisticated ones using Runge-Kutta algorithms. MATLAB has several different functions (built-ins) for the numerical solution of ordinary differential equations (ODE). The problems of solving an ODE are classified into initial-value problems (IVP) and boundary-value problems (BVP), depending on how the conditions at the endpoints of the domain are specified [6,7]. All the conditions of an initial-value problem are specified at the initial point. On the other hand, the problem becomes a boundary-value problem if the conditions are needed for both initial and final points [8,9]. An ordinary differential equation (ODE) is an equation that contains an independent variable, a dependent variable, and derivatives of the dependent variable. Many studies have not shown the solution to two-step parallel pyrolysis kinetics involving the prediction of pyrolysis products due to complex nature of the equations, and too many unknown variables. It has been established that the overall process of pyrolysis appears simple but the sequence of reactions is complex and involves both endothermic and exothermic processes whose thermodynamics and kinetics are poorly understood [10,11]. Under such complex phenomena, it is impossible to formulate a complete mathematical model of pyrolysis which will still be mathematically tractable. Yang et al. [12] presented that the major stage of biomass pyrolysis occurs between 250 and 450 °C. The two-step parallel reaction model of biomass pyrolysis has previously been used by other researchers [12][13][14][15][16][17][18] as reported by Sobamowo et al. [11]. Again, Babu and Chaurasia [19] studied numerically the effects of heating conditions, heating rates and order of the kinetic reactions on wood pyrolysis using the twostep model of Koufopanos et al. [20]. The model developed by Sobamowo et al. [11] it was based on Shafizadeh and Chin's model [13] which has been stated to be the most classical models for wood pyrolysis [21], the model was solved using the developed Runge-Kutta fifth-order for the simultaneous equations following the analytical method, only primary products where predicted. Most authors have focused on wood pyrolysis for primary reaction only, heating rate and effect on biomass particle, only very few reports on prediction of pyrolysis products for primary products especially [char, tar (bio-oil), and gas]. Sobamowo et al. [11] studied the pyrolysis of shrinking biomass particle in fixed bed gasifier and predicted the product yields (char, gas, and tar) using wood biomass. Presently no report on prediction of pyrolysis products under two-step parallel pyrolysis reaction scheme, which forms the novelty of this research.
This study has simplified the steps involved to implement a two-step parallel pyrolysis kinetics to predict the maximum char yield during primary and secondary reaction (char 1 , gas 1 , tar(bio-oil), char 2 and gas 2 ) using ODE45 MATLAB function. ODE45 is based on an explicit Runge-Kutta formula, the Dormand-Prince pair [22]. The numerical solver ODE45 combines a fourth order method and a fifth order method. This method varies the step size, choosing the step size at each step in an attempt to achieve the desired accuracy. In general, ODE45 is the best function to apply for most problems [23]. The purpose of study is to show the details of implementing the typical steps of numerical solution, so that it will be clear exactly what computations are being executed.

Materials and methods
In this study, a two-step parallel pyrolysis model was implemented through MATLAB ODE-45 algorithm to solve a set of differential equations (see Eqs. 2-7). The ODE 45 is function handle that evaluates the right side of the differential equations [24]. All solvers solve systems of equations in the form y′ = f(t, y) or problems that involve a mass matrix, M(t, y)y′ = f(t, y). The ode 45 solver can solve only equations with constant mass matrices [25]. The analytical model consisting of 6 partial differential equations was simulated in MATLAB 2019 with ODE 45 function which showed a plot of mass loss rate for the biomass and the various products. A set of differential equations for predicting parallel pyrolysis kinetics has been proposed by Di Blasi [26].

Pyrolysis kinetics
According to Di Blasi [26], the primary reaction is the thermal degradation of biomass into gas, tar, and char; eventually, when the temperature reaches a higher value, (see Fig. 1) tar decomposes into gas and char in secondary reactions. The first three reactions are endothermic, while the fourth and fifth are exothermic. The mechanism adopted for kinetic modeling in this study is presented in Fig. 1. The dependence on temperature is intrinsically introduced with the equilibrium constants k i of primary (k 1 , k 2 , k 3 ) and secondary (k 4 , k 5 ) reactions. They may be expressed with the Arrhenius correlation [27]: A is the frequency factor, E is the activation energy, R is the gas constant, and T is the temperature at which the reaction takes place.
According to the two-stage parallel reaction model shown above in Fig. 1, the wood biomass undergoes thermal degradation according to primary reactions (k 1 ; k 2 ; k 3 ) giving gas 1 , tar, and char 1 as products at T = 473.15 K where the biomass absorbs heat energy (endothermic reaction) from the surroundings to form products. However, tar may undergo secondary reactions (k 4 , k 5 ) at T > 473.15 K where heat energy is released (exothermic reaction) to produce gas 2 and char 2 . The yields of the different products ranging from 1 to 5, and the contemporary diminution of the biomass is a function of both time and temperature. Mass loss with time may be expressed as shown in Eqs. (2-7) according to Fantozzi et al. [27]: The energy required to carry out the above pyrolysis reactions can be expressed as the product of the mass-produced and the heat of reaction [27]: (3) dm gas1 dt = m Biomass k 1 .
(7) dm char2 dt = k 5 m tar . where, the index i refers to the five reactionsE and m were the energy required by, or released from, the ith reaction and the mass-produced, while h is the heat of reaction (negative for i = 1 to 3 and positive for i = 4 to 5). Since this study adopts a two-step parallel kinetics mechanism, where tar also known as bio-oil decomposes into gas 2 and char 2 . The above equation can be can be expressed in terms biochar and gas obtained from primary and secondary reactions, respectively as: Additionally, biochar yield was calculated using Eq. (9) below and expressed in mass percentage basis (%) according to Nurhidayah et al. [28] as: Initial mass of biomass of 0.5 kg was considered for the numerical simulation. The kinetic parameters sourced from past works (see Table 1), it shows the constant parameters used for the simulations of the kinetics and heat transfer during the pyrolysis of wood. Since, thermogravimetric analysis (TGA) was carried out on pyrolysis of wood to predict the kinetic parameters. The simulation in the present study was considers pyrolysis of wood based on availability of kinetic data from TGA for a two-step parallel reaction.
In Table 1, A represents the frequency factor, E is the activation energy, H is the heat of reaction and R is gas constant. The rate of reaction constants was computed in MATLAB 2019 using Arrhenius equation (Eq. 1) using the activation energy, frequency factor and gas constant from Table 1 for the primary and secondary reactions, respectively. These constant values were selected based on slow pyrolysis since the target by-product is biochar. The computed rate of reaction constants, initial mass of the biomass and reaction time of 5400-9000 s were entered into the set of differential equations (Eqs. 2-7). The set of differential equations was solved by an algorithm developed in this study using the MATLAB ODE 45 function. The following procedure has been designed to implement the algorithm.

Algorithms
These lines are steps undertaken to generate the codes for the model implementation and kinetic parameters used in the study.

Algorithm 1
This section calculates the kinetic parameters for the primary and secondary reactions at varying temperatures with a user-defined function.

Algorithm 3
This part contains the simulation codes using ODE 45 function, it is divided into sub-sections to showing the temperature variations. The flowcharts below describe the steps or algorithm 1 and 2, respectively in a graphical form, as well as the main program (see Fig. 2).

Results and discussion
The following rate constants at varying temperatures as computed in MATLAB is shown in Table 2. These rate of reaction constants were used to compute the residual mass of the wood biomass at time, t = 0-9000 s (see predicted residual mass of pyrolysis products in Fig. 2). It was observed that the rate of reaction increased with temperature. The rate constants in the present study as shown in Table 2 are in agreement with that of a multi-step thermogravimetric analysis and kinetics modeling of isothermal carbonization of wood in an inert atmosphere in the works of [29], which ranged from 0.00000122 to 0.0104 s −1 . The authors did not predict pyrolysis products, only kinetics and rate constants were reported. The following graphs shows the mass of pyrolysis product yield at different temperatures and time intervals (see Fig. 3). The graphs were obtained from the MATLAB ODE solver implemented in the software to solve the set of differential equations in Eqs. (2)(3)(4)(5)(6)(7). As seen in Fig. 1, the initial mass of the biomass is 0.5 kg which diminishes as pyrolysis takes place from initial time t (0 s) to the final time as well as the change in temperature. Additionally, as the mass of biomass decreases the pyrolysis products increases until tar is finally consumed at secondary reaction (see Table 3) which shows that the values of the product yield.
In Fig. 3, it was observed that, as the mass of biomass decrease with time, the pyrolysis products such as char, tar (biooil) and gas increases. However, gas 1 and gas 2 increased rapidly due to the effect of temperature. Therefore, at higher temperatures the mass of gas increased while char decreased. At 673.15 K, tar (bio-oil) is consumed completely, mass of tar ≅ 0 (see Fig. 4b, c, e, f ). This is in agreement with the two-step pyrolysis reaction scheme in Fig. 1. The negative sign in Table 3 indicates exothermic reaction which leads to an increase in the entropy of the surroundings. An exothermic reaction occurs when the temperature of the system increases due to the evolution of heat. This confirms the modified two-step pyrolysis reaction scheme in Fig. 1 which is in agreement with the works of Di Blasi [26] and Fantozzi et al. [27]. In Table 3 the %biochar and gas yield produced during the simulation is in tandem with the findings of Ojolo et al. [30] as shown in Fig. 4.
In Table 4, we observed tar (bio-oil) is decomposed to form gas 1 and gas 2 (i.e., the percentage of gas in superscript of a' and b') during primary and secondary reaction scheme (see two-step pyrolysis reaction scheme in Fig. 1). The findings of the present study agree strongly with that of Ojolo et al. [30] who reported that when temperature increases above 473 K, biomass begins to decompose into char, tar (bio-oil) and gaseous products until it reaches a temperature 723 K after 225 s of residence time. That the production rate of char is higher than that of tar and gas. This validates the values given in two-step pyrolysis reaction scheme of this study.

Validation of the algorithm
The final, but censorious step in modeling is validation [31]. This step is accomplished by comparing the mass loss rate obtained from the pyrolysis experiments to those predicted by numerical simulation. This methodology has been successfully applied to a number of non-charring [32,33] as well as charring [34,35] processes. The model developed in the present study is compared with the isothermal biomass loss experiments conducted by Lina Samuelsson at Lulea University [36,37] and the model developed by reaction model developed by Ojolo et al. [30]. Although, a comparison between the mass of pyrolysis products in the present study using numerical method and the concentration of products of during pyrolysis using analytical method by Ojolo et al. [30] showed very close resemblance in graphs plotted. The curves in the graphs showed huge similarity (see Fig. 4). However, the differences in the plots can be account for the residence time and number of pyrolysis products predicted.
The experimental procedure used by Samuelsson [36] for the pyrolysis of wood is discussed below. In the experiments, Norway spruce wood chips of three different thicknesses were prepared. The TGA setups did not include any instrumentation for gas analysis which limits us to predict solely the solid mass fraction as a function of temperature and time. The installed auto-sampler, together with the fast cooling of the furnace after finished runs make this instrument suitable for kinetic studies needing many repetitions at several heating rates. After each pyrolysis experiment the sample basket containing the wood chips was withdrawn from the furnace and placed in the purge zone to cool down before removing the char from the basket. The mass of the pre-dried sample chip and the resulting char was measured in an external balance (LPW-723, VWR). For all chip thicknesses the sample mass measured in the external balance varied less than 10 mg from the value measured in the macro-TGA. The experiment was performed at each isothermal temperature, i.e., at 574, 597, 615, 637 and 676 K. The minimum dwelling time for the medium and thick chips was 44 min and the longest time was 240 min. At all temperatures apart from 573.15 K, the model mass loss rates are higher than the experiments. At 673.15 K this temperature, the numerical simulation using ode45 function gives a char yield of 27.84% with the yield from the experiment by Samuelsson [36] been 22% (error = − 0.26). The maximum error from experimental results of past works is  To evaluate the accuracy and performance of the algorithm, predicted values was compared with experimental data from literature using the percentage error as a function of error using Eq. (10) according to [38,39]: A comparison of the rate of production of pyrolysis products measured in the experiments obtained from past work (secondary data) to that predicted by the model solved using ODE45 function was used to validate the model's parameters. This comparison was made using only primary %char yield (char 1 ), since no work has been done on prediction of secondary char 2 . This rate, also referred to as mass loss, multiplied by the effective heat of combustion of these products, yields the heat release rate, which is the single most important quantity characterizing material flammability [40]. Representative comparisons between the experimental work from literatures and simulated mass loss rates. All mass loss rates demonstrate a reasonable agreement with other studies.

Conclusion
MATLAB allows for a very simple implementation of the solution to differential equations involving two-step parallel pyrolysis of biomass and ODE45 functions as discussed in this study. Many researchers have not looked into the solution technique for a two-step parallel pyrolysis reaction to predict the by-products. In the present study, we use the results of the pyrolysis products predicted from the numerical simulation to compare with experiments from past works, a good agreement was achieved. Although, the kinetic equations used in this study are that of slow pyrolysis for wood biomass. The algorithm developed in this work is simplified, easy to use, and can be applied to any type of biomass undergoing parallel pyrolysis reaction if the kinetic parameters are known. Therefore, we have proposed this algorithm for the solution to pyrolysis models involving a two-step kinetic scheme which we have showed how it can be efficiently solved using MATLAB ODE45 function.