Universal Kinetic Model to Simulate Two-Step Biodiesel Production from Vegetable Oil

: To date, to simulate biodiesel production, kinetic models from di ﬀ erent authors have been provided, each one usually applied to the use of a speciﬁc vegetable oil and experimental conditions. Models, which may include esteriﬁcation, besides transesteriﬁcation simulation, were validated with their own experimental conditions and raw material. Moreover, information about the intermediate reaction steps, besides catalyst concentration variation, is either rare or nonexistent. Here, in this work, a universal mathematical model comprising the chemical kinetics of a two-step (esteriﬁcation and transesteriﬁcation) vegetable oil-based biodiesel reaction is proposed. The proposed model is universal, as it may simulate any vegetable oil biodiesel reaction from the literature. For this purpose, a mathematical model using the software MATLAB has been designed. Using the mathematical model, the estimation of mass variation with time, of both reactants and products, as well as glyceride conversion and homogeneous catalyst concentration variation (instead of only alcohol / catalyst solution) are allowed. Moreover, analysis of the inﬂuence of some important variables a ﬀ ecting the reaction kinetics of biodiesel production (e.g., catalyst concentration), along with comparison and model validation with data from di ﬀ erent authors may be carried out. In addition, Supplementary material with a collection of 290 rate constants, derived from 55 di ﬀ erent experiments using di ﬀ erent vegetable oils and conditions is provided.


Introduction
In the global effort to generate non-contaminant renewable fuels, biofuels comprising biodiesel, bioethanol, biobutanol and biogas, among others, play a significant role. Nowadays, many researchers work to improve processing efficiency and the integration of cellulose (to date, expensive) and organic waste recycling, as raw materials, in economically sound processes [1][2][3].
Biodiesel is a green fuel that can be used in diesel engines instead of fossil fuels. It is mainly produced through transesterification of vegetable oils/fats or animal fats, including further purification steps. Biodiesel is then used in diesel engines as a substitute for diesel fuel, straight or mixed with diesel fuel. Furthermore, it is well known that, during the process of combustion, biodiesel exhaust emissions, including carbon monoxide, hydrocarbons and particulates, are lower than those of diesel fuel, providing zero emissions of sulfur, lead and other heavy metals, besides the added advantage of its powerful lubricating properties [4].
The debate on food versus fuels focuses the discussion on the influence of biofuel prices on food prices and shortages, that are especially significant during a food crisis [5]. However, the crude oil price seems to be the main cause of food price increases [6]. Moreover, traditionally, land has been used to grow nonfood crops, e.g., tobacco or crops to produce alcoholic beverages, so its use to produce energy should be considered at the same level of acceptance. In any case, if land is not used to produce energy, food crop residues may still be used to provide biofuels, thus minimizing the polluting effect of these wastes [7,8].
In the search for new alternative low-cost raw materials, the reuse of edible frying oil may become an interesting raw material. Although this oil exhibits a high amount of free fatty acids (FFA), it may be used to produce biodiesel with minor pretreatment [9]. Other affordable commodities, which may be used as feedstock to produce biodiesel, are alternative vegetable oils, e.g., Jatropha curcas oil, also characterized by a high content of FFA. This represents a key problem for the common alkaline transesterification process because the alkaline catalyst may react with FFA, thus producing soap (saponification reaction), with consequent reduction of biodiesel yield and compounding separation of esters, glycerol and washing water. Soap formation also increases viscosity and leads to the formation of a gelatinous final product. In conclusion, the use of a raw material with an FFA content above 0.5% is not recommended for alkali-catalyzed transesterification reactions; therefore, a previous step consisting in FFA esterification with an acid catalyst (typically, H 2 SO 4 ) is required. This previous step is followed by a second one, consisting in the transesterification of triacylglycerols by alcoholysis, usually with a basic catalyst.
During esterification, FFAs react with an alcohol to produce fatty acid ester and water. Subsequently, in the transesterification reaction, triacylglycerols react with an alcohol to produce other esters of the same fatty acids and glycerol. This reaction can be written as three consecutive reversible reactions with intermediate formation of diacylglycerols and monoacylglycerols. Finally, biodiesel is separated from glycerol by decantation or centrifugation [10]. Methanol is the most commonly used alcohol and, in this case, the reaction is known as methanolysis. There are physical and chemical properties that affect both esterification and transesterification reactions, i.e., the source of triacylglycerols, reaction temperature, type and amount of catalyst, type and amount of alcohol, stirring time, etc. [11]. Although this is a well-known process, the ratio of reactants affects the process in terms of conversion efficiency and this factor differs depending on the raw material chemical composition.
There are several kinetic studies about biodiesel production from vegetable oils under different reaction conditions. Esterification using H 2 SO 4 as a catalyst has been carried out with methanol and either coconut or sunflower oil, and with ethanol and sunflower oil, at different temperatures and concentrations [12,13]. In the same way, methanolysis has been carried out using KOH or NaOH as catalyst, using rapeseed oil [14], sunflower oil [15,16], Pongamia pinnata oil [17], Brassica carinata oil [18], palm oil [19] or soybean oil [20,21] also at different temperatures and concentrations. However, only a few kinetic studies about enzyme-or lipase-assisted transesterification have been reported [22][23][24]. More recently, the use of new catalysts, e.g., barium cerate [25], metal oxide mixed Sr-Ce [26], Li/NaY (zeolite) [27], lithium-based chicken bone (Li-Cb) composite [28], calcium oxide from seashell [29], graphene [30] or solid acidic ionic liquid polymer [31] have been included in kinetic studies. Despite the potential of microalgal and microbial oil for biodiesel production [8], kinetic models have not yet been developed.
After a thorough analysis of previous studies about the kinetics of vegetable oil esterification and transesterification [12][13][14][15][16][17][18][19][20][21], it may be noticed that proposed mathematical models are generally similar and based on chemical stoichiometry, providing mass variation along reaction time. The input variable is the reaction rate constant, which depends on reaction conditions. Therefore, to work with these models, knowledge of these constants to further validate the model is mandatory. Current studies include kinetic studies on the synthesis of fuel additives from glycerol, biodiesel/hexanol blends or by means of neural networks [32][33][34].
The present work aims to provide a universal mathematical model, to simulate a two-step homogeneously catalyzed biodiesel reaction, that may be of application to any vegetable oil and reaction conditions. Moreover, intermediate conversion steps will be provided, besides any concentration variation of reactants and products. The present study is based on the characterization of reaction kinetics, while the objective is to show, through mathematical modeling, the influence of the most important variables affecting the reaction kinetics of biodiesel production, which are collected into constants of reaction rate. In this sense, models are an indisputable tool to both assist decision-making in the production process of biodiesel and predict reaction yield. However, to collect information about experiment conditions, previous laboratory tests must be carried out. To the best of our knowledge, a universal kinetic model that may be applied to any vegetable oil and reaction conditions is missing. Moreover, catalyst concentration variations are seldom mentioned, thus making difficult the comparison between models or model validation with data from literature.

Analysis of Process Variables
In this work, a kinetic model based on reaction stoichiometry is proposed. One mole of triacylglycerol (TAG) and three moles of alcohol (A) react to produce three moles of fatty acid esters (FAE) and 1 mole of glycerol (Gly). The general reaction is named transesterification.
To reduce the order of the model equation, while describing the intermediate steps, the previous equation is divided into three consecutive reversible reactions. At each step, one mole of FAE, diacylglycerols (DAG), monoacylglycerols (MAG) and Gly is released.
In the case of a raw material with a large presence of FFA, a previous esterification reaction is needed. In this situation, one mole of FFA reacts with one mole of A to produce one mole of FAE, plus water (W).
While the catalysts used for esterification are acids (H+), the preferred transesterification catalysts are alkalis (OH-). In fact, during biodiesel production, esterification is carried out under acidic medium (H 2 SO 4 ) while transesterification usually takes place in a basic medium (NaOH or KOH); both processes being completely independent. Moreover, both reactions cannot take place simultaneously. This ends up with two processes that cannot meet either physically or mathematically.
The reaction rate varies proportionally with the involved species concentrations; the proportionality constant is known as the speed constant, k i . This constant depends on the type of oil, alcohol type, temperature, presence and type of catalyst. It is important to note that the rate constant, k, only depends on the catalyst, whatever the concentration may be.
Although stirring speed is usually not considered among variables, it should be included among parameters affecting processing. Variations of the amount of oil, alcohol, catalyst and temperature affect the reaction rate, although the rate constant remains unchanging.
Previous kinetic studies proposed that oil and alcohol type, alcohol-to-oil molar ratio, catalyst type and concentration, and reactor stirring speed were parameters to take into consideration [12][13][14][15][16][17][18][19][20]. Then, performing the reaction at different temperatures provided information to further calculate activation energy and, thus, rate constants.
A review of the literature shows that, when mathematical models are proposed and validated with experimental data, the catalyst concentration is seldom or never included in the reaction rate equation. To gain knowledge about the implications of this decision, and considering first order (with respect to OH − ) transesterification reaction, it may be found that by doubling the catalyst concentration, reaction rate would also be doubled. That is, catalyst concentration may increase or decrease reaction rate and, therefore, reaction time, which is really important in the production of biodiesel. Furthermore, as the rate constant has usually been considered independent of reagent concentration, then, alcohol-to-oil molar ratio should also be independent. However, a dependence between rate constant and alcohol-to-oil molar ratio may be found, as the last one depends on catalyst type and concentration [11].

Mathematical Analysis
Using previous equations and in agreement with the above conditions and hypothesis, a series of expressions based on mass balance and their derivatives with respect to time are proposed. These expressions allow the inclusion of the change in concentration of both reagents and final products with time. These expressions can be represented by the following differential Equations (1): where, k 1 , k 3 , k 5 and k 7 are reaction effective rate constants and k 2 , k 4 , k 6 and k 8 are reverse reaction constants. The chemical elements located between brackets show the concentration of the different species.

Mathematical Model
To solve the set of differential Equation (1) and mathematical calculations, the software MATLAB 5.3 Academic Use (The MathWorks Inc., Natick, MA, USA) was used and toolbox "Ordinary Differential Equations" (ODE) solver was selected. For this purpose, allocation of states representing Equation (1) was performed. This tool allows ordinary differential equations to be converted into a system of equations of state, where each state corresponds to the concentration of each reagent or product involved in the kinetics. Thus, the model is transformed into a third-degree equation system (Expression (2)), depicting the same number of equations and unknown variables.
−k 1 x 1 x 2 x 9 + k 2 x 3 x 4 x 9 = 0; −k 1 x 1 x 2 x 9 + k 2 x 3 x 4 x 9 − k 3 x 2 x 5 x 10 + k 4 x 3 x 6 x 10 − k 5 x 2 x 6 x 10 + k 6 x 3 x 7 x 10 − k 7 x 2 x 7 x 10 + k 8 x 3 x 8 x 10 = 0; k 1 x 1 x 2 x 9 − k 2 x 3 x 4 x 9 + k 3 x 2 x 5 x 10 − k 4 x 3 x 6 x 10 +k 5 x 2 x 6 x 10 − k 6 x 3 x 7 x 10 + k 7 x 2 x 7 x 10 − k 8 x 3 x 8 x 10 = 0; k 1 x 1 x 2 x 9 − k 2 x 3 x 4 x 9 = 0; −k 3 x 2 x 5 x 10 + k 4 x 3 x 6 x 10 = 0; k 3 x 2 x 5 x 10 − k 4 x 3 x 6 x 10 − k 5 x 2 x 6 x 10 + k 6 x 3 x 7 x 10 = 0; k 5 x 2 x 6 x 10 − k 6 x 3 x 7 x 10 − k 7 x 2 x 7 x 10 + k 8 x 3 x 8 x 10 = 0; k 7 x 2 x 7 x 10 − k 8 x 3 x 8 x 10 = 0; − k 1 x 1 x 2 x 9 + k 2 x 3 x 4 x 9 = 0; − k 3 x 2 x 5 x 10 + k 4 x 3 x 6 x 10 − k 5 x 2 x 6 x 10 + k 6 x 3 x 7 x 10 − k 7 x 2 x 7 x 10 +k 8 x 3 x 8 x 10 = 0; (2) where allocation of states corresponds to Equation (3): To simulate the previous model, a rigid system solver (ODE15 and ODE45, included in ODE solver, from MATLAB 5.3 Academic Use) based on a numerical differentiation formula of variable order and multiple width or step was used. ODE15 function is recommended to solve rigid differential equations, while ODE45 function is time consuming due to the rigidity of ordinary differential equations. This rigidity is due to the presence of some terms in equations that can lead to a rapid and unstable variation of solution, unless an extremely small step size is selected. In general, function ODE45 is more accurate, although it consumes more computing time to perform numerical integration. To develop the mathematical model, experimental data from literature were used.

Results and Discussion
The proposed model may be positively validated with any previous work from literature. Results from proposed model simulation, with data from literature, are shown below.

Rate Constants
As mentioned above, as a first step prior to experimentation, rate constant determination was required. In this work, a set of experimental data from different authors corresponding to rate constants, depending on experimental conditions (the most important variables affecting kinetics of biodiesel production) were studied. A small collection of rate constants for different vegetable oils with their respective variables is provided as Supplementary material. These constants are nominated as k 1 , k 2 , k 3 , k 4 , k 5 , k 6 , k 7 and k 8 , following the criteria established in the proposed model and constitute a small database for virtual experimentation with this kind of mathematical model.

Concentration of Reactants and Products
To validate the model, three experiments provided by different authors from the literature, with different conditions (experiment no. 1 [15], experiment no. 2 [16] and experiment no. 3 [12]), have been analyzed and applied to the proposed model.

Model Validation with Experiment No. 1
Kinetics of sunflower oil methanolysis at 60 • C, 400 rpm stirring speed and a catalyst concentration (NaOH) of 0.5% (w/w) are shown in Table 1 [15]. For model validation purposes, simulation results from this experiment, together with experimental data from Table 1, are shown in Figure 1. As may be seen, the model suited experimental data, providing information that corresponded to the simulated reality (see Section 3.7 for statistical analysis). Figure 1 exhibits the concentration vs. time of each compound, showing that reaction is completed after the first half-hour. Table 1. Kinetics of sunflower oil methanolysis at 60 • C, 400 rpm stirring speed and a catalyst concentration (NaOH) of 0.5% (w/w) [15]; t: time, TAG: triacylglycerol concentration, DAG: diacylglycerol concentration, MAG: monoacylglycerol concentration, FAME: fatty acid methyl ester concentration, Gly: glycerol concentration, A: alcohol concentration.

Validation of the Model with Experiment No. 2
In this section, the aims were both to validate the model with other experimental data and to check the effect of the catalyst concentration variation in the proposed model. In this sense, Table 3 shows kinetics from methanolysis with sunflower oil at 25 • C, 600 rpm stirring speed and a catalyst concentration (KOH) of 1.5% (w/w), according to Vicente et al. [16].  For model validation purposes, simulation results are shown in Figure 2, together with experimental data shown in Table 3. It may be inferred that the model appears to fit to the experimental data from experiment no. 2, provided by [16], thus validating the proposed model (see Section 3.7 for statistical analysis).  Table 4 shows the effect of catalyst concentration (KOH) in the production of sunflower oil biodiesel, using methanol at 25 °C and 600 rpm stirring speed, according to results from [16]. Table 5 shows the rate constants used for simulation. Simulation results and experimental data (from Table  4) are shown in Figure 3. It may be inferred that the model appears to fit to experimental data from experiment no. 2, provided by [16], thus validating the proposed model (see Section 3.7 for statistical analysis). Figure 3 shows the effect of catalyst concentration variation included in the proposed model. Table 4. Effect of catalyst concentration (KOH) in the production of sunflower oil fatty oil methyl ester (mol l −1 ) and methanol (25 °C y 600 rpm) [16].

Time (min)
Catalyst Concentration % (w/w) 0. 5 Table 4 shows the effect of catalyst concentration (KOH) in the production of sunflower oil biodiesel, using methanol at 25 • C and 600 rpm stirring speed, according to results from [16]. Table 5 shows the rate constants used for simulation. Simulation results and experimental data (from Table 4) are shown in Figure 3. It may be inferred that the model appears to fit to experimental data from experiment no. 2, provided by [16], thus validating the proposed model (see Section 3.7 for statistical analysis). Figure 3 shows the effect of catalyst concentration variation included in the proposed model. Table 4. Effect of catalyst concentration (KOH) in the production of sunflower oil fatty oil methyl ester (mol L −1 ) and methanol (25 • C y 600 rpm) [16].

Time (min)
Catalyst Concentration % (w/w)  Table 5. Reaction conditions and rate constants used in sunflower oil biodiesel simulation, according to conditions provided by [16].

Parameters Values
Temperature  Table 5. Reaction conditions and rate constants used in sunflower oil biodiesel simulation, according to conditions provided by [16].

Validation of the Model with Experiment
No. 3 Table 6 shows results from M. L. Pisarello et al. [12] after sunflower oil esterification with ethanol at 70 °C and a concentration of catalyst (H2SO4) of 0.2% (v/v). On the other hand, Table 7 shows rate constants used for simulation, according to data conditions provided in Table 6.

Validation of the Model with Experiment
No. 3 Table 6 shows results from M. L. Pisarello et al. [12] after sunflower oil esterification with ethanol at 70 • C and a concentration of catalyst (H 2 SO 4 ) of 0.2% (v/v). On the other hand, Table 7 shows rate constants used for simulation, according to data conditions provided in Table 6.  Table 7. Reaction conditions and rate constants used in the simulation according to data provided by [12]; n.d.: no data.
For model validation purposes, simulation results from this experiment are shown in Figure 4, together with experimental data from Table 6. It may be seen that the theoretical data calculated from the model (representing concentration variation of W and FFA with time) are in good agreement with experimental data from [12]. However, an exception is provided by fatty acid ethyl esters (FAEE), for which the comparison between last simulated and experimental values seem to differ: this may be explained by an experimental calculation error (see Section 3.7 for statistical analysis).
Energies 2020, 13, x FOR PEER REVIEW 10 of 14 For model validation purposes, simulation results from this experiment are shown in Figure 4, together with experimental data from Table 6. It may be seen that the theoretical data calculated from the model (representing concentration variation of W and FFA with time) are in good agreement with experimental data from [12]. However, an exception is provided by fatty acid ethyl esters (FAEE), for which the comparison between last simulated and experimental values seem to differ: this may be explained by an experimental calculation error (see Section 3.7 for statistical analysis).

Simulation to Analyze the Effect of Catalyst Concentration
Considering the proposed model, change in concentration of transesterification reactants and products, after 1, 5 and 10 min of reaction time, have been simulated. Results are shown in Figure 5, while simulation conditions are listed in Supplementary material, considering an initial concentration of methanol and triacylglycerols of [A] = 12 mol/l and [TAG] = 2 mol/l, respectively. As expected, simulation shows that, for the same time instants, the higher the amount of catalyst, the higher the reaction rate. It may also be inferred that, for a catalyst concentration around 1% (w/w),

Simulation to Analyze the Effect of Catalyst Concentration
Considering the proposed model, change in concentration of transesterification reactants and products, after 1, 5 and 10 min of reaction time, have been simulated. Results are shown in Figure 5, while simulation conditions are listed in Supplementary material, considering an initial concentration of methanol and triacylglycerols of [A] = 12 mol/L and [TAG] = 2 mol/L, respectively. As expected, simulation shows that, for the same time instants, the higher the amount of catalyst, the higher the reaction rate. It may also be inferred that, for a catalyst concentration around 1% (w/w), the conversion of [TAG]

Statistical Analysis of the Model
A summary of the statistics of the experiments is shown in Tables 8 and 9. As may be seen from Table 9, model and experimental values are in good agreement. The model goodness of fit has been calculated by means of the coefficient of determination (R-squared) that shows the proportion of

Statistical Analysis of the Model
A summary of the statistics of the experiments is shown in Tables 8 and 9. As may be seen from Table 9, model and experimental values are in good agreement. The model goodness of fit has been calculated by means of the coefficient of determination (R-squared) that shows the proportion of variation of results that can be explained by the model. In other words, a value of 1.0 indicates a highly reliable model, while a value close to 0.0 means a failing model. In this sense, the correlation between experimental data and simulated values (Figures 1-4) results in 0.92184506, 0.99999914, 0.99232725 and 0.98553624, respectively. As may be seen, R 2 value is, in all cases, above 92%, thus showing the goodness and reliability of the model.

Conclusions
Mathematical models are versatile tools for a low-cost study, reducing raw material consumption and experimentation instrument needs, among others. In this sense, models for biodiesel production are increasingly widespread. Models found in literature consider the use of a solution of catalyst and alcohol; therefore, experiments with separated reagent concentration variations are not modeled. In this work, as a novelty, the proposed model considers the presence and variations of homogenous catalysts and alcohol, separately, allowing further simulations.
Comparing studies from literature with the proposed model, it may be seen that the reaction depends on a number of variables, including the type of vegetable oil, reaction stirring speed, type and concentration of catalyst, type of alcohol, alcohol-to-oil molar ratio, reaction time and temperature. These parameters are crucial and differentiating elements in the transformation of oils into biodiesel.
To increase the usefulness of kinetic models from literature, a universal model of biodiesel production that allows implementing control strategies is provided. In this sense, the model could provide information about best reaction time to achieve maximum biodiesel production, optimum catalyst concentration for different reaction times, etc., allowing automated production, thus reducing production costs.
The biggest drawback presented by previous models found in literature is the need, before starting simulation, of the rate constant determination, which means that experimental work is needed. For this reason, a collection of rate constants, derived from different experiments using different vegetable oils and conditions is provided as Supplementary material. This collection allows simulation of a potential production plant, thus helping to find out production deviations.