Pseudo-Homogeneous and Heterogeneous Kinetic Models of the NaOH-Catalyzed Methanolysis Reaction for Biodiesel Production

Abstract

Methanolysis of vegetable oils in the presence of homogeneous catalysts remains the most important process for producing biodiesel. However, there is still a lack of accurate description of the reaction kinetics. This is in part due to the complexity of the reacting system in which a large number of interconnected reactions take place simultaneously. In this work, attention is focused on the biphasic character of the reaction medium, formed by two immiscible liquid phases. The behavior of the phases is investigated regarding their physicochemical properties, mainly density and mutual solubility of the components, as well as composition. In addition, two kinetic models with different level of complexity regarding the biphasic character of the reaction medium have been developed. It has been found that a heterogeneous model considering the presence of the two phases and the distribution of the several compounds between them is indispensable to get a good description of the process in terms of oil conversion and products yields. The model captures the effects of the main variables of an isothermal batch methanolysis process: methanol/oil molar ratio, reaction time and catalyst concentration. Nevertheless, some adjustment is still required as concerns modelling of the saponification reactions and catalyst deactivation.

1. Introduction

Triglyceride transesterification is the basis of biodiesel synthesis from vegetable oils. Methanol is the most commonly used transesterification reagent; in this case, the reaction is also known as methanolysis; the fatty acid methyl esters (FAMEs or MEs) formed constitute biodiesel. As for the catalyst, in most industrial biodiesel production processes, basic compounds such sodium methoxide or hydroxide are used. These compounds stand out for being very active due to their strong basicity and solubility in methanol, and for this reason, they are categorized as homogeneous catalysts [1].

Knowing the kinetics of the methanolysis reaction is key for obvious reasons related to the design, scale-up, operation, and control of industrial reactors. Several researchers have performed kinetic studies on the vegetable oil transesterification reaction [2,3]. In addition, some kinetic modelling approaches have been reported. A pseudo-homogeneous model was proposed by Noureddini and Zhu [4] that included the three consecutive and reversible transesterification reactions corresponding to the complete conversion of triglycerides into biodiesel and glycerol. It was found that a sigmoid curve described the evolution of the oil conversion with reaction time. This behavior was ascribed to the existence of two distinct regimes during the progress of the reaction. Initially, the reaction takes place under a mass-transfer-limited regime, followed by a regime controlled by the intrinsic reaction kinetics. Pseudo-homogeneous models are characterized by disregarding the biphasic nature of the reacting system which, actually, consists in an emulsion formed by droplets containing methanol and glycerol dispersed in a continuous phase constituted by the partially converted triglycerides and fatty acid methyl esters. This model has been adopted in several previous woks [5,6].

Stamenkovic et al. [7], on the other hand, developed a heterogeneous model considering the presence of two immiscible liquid phases and three reaction regimes. The first and slowest regime corresponded, like in the previous case, to a mass-transfer-limited step due to the biphasic nature of the system that can be avoided by operating the reactor at a sufficiently high temperature. Successively, the process was controlled by intrinsic chemical kinetics and thermodynamic equilibrium. A drastic reduction of the emulsion droplet size with reaction time was noticed.

A variety of kinetic studies can be found in the literature in which important simplifications were adopted, such as considering the reactions irreversible [8,9,10] or neglecting the presence of partially converted intermediates (i.e., diglycerides and monoglycerides) [11]. In addition, many kinetic models reported in the literature do not take into consideration explicitly a very important variable such as the catalyst concentration, though some relevant studies do include it in the reaction rate equations [12,13]. Saponification reactions of glycerides and methyl esters have been also included in other models [14,15,16].

Clearly, the heterogeneous character of the reaction medium, in which one of the reacting phases is dispersed in the other one forming an emulsion [17,18], has important consequences. In particular, the interfacial surface area, which depends on the size of the dispersed phase droplets, and the rate of mass transfer through the liquid–liquid interphase become essential for developing a model that properly describes the reaction kinetics [19,20].

In the present work, the kinetics of the methanolysis of sunflower oil catalyzed by NaOH is investigated with the aim of developing models with different level of complexity that can be suitable for different levels of description of the process. In this regard, pseudo-homogeneous and heterogeneous models considering the distribution of the components among the phases and the consumption of the catalyst by saponification reactions are developed.

Although the use as feedstock for biodiesel production of low-cost non-edible vegetable oils is highly desirable [21], in this work, refined sunflower oil was used. As the aim of the study was to develop an adequate kinetic model, refined sunflower oil is particularly suitable for this purpose both because its composition is close to that of relevant non-edible oils such as jatropha curcas regarding the fatty acids profile, containing mainly glycerides of the linoleic and oleic acids, as well as low acidity [22], and because it is easily available with the required purity to perform a kinetic investigation without uncontrolled interferences.

2. Materials and Methods

2.1. Transesterification Reactions

Transesterification reactions were carried out in a batch 1 L jacketed glass reactor fitted with a reflux condenser, a nitrogen gas inlet, a stainless-steel stirrer comprising a turbine, a thermocouple probe, and a heated circulating water bath for reaction temperature control. Kinetic runs were conducted at 323 K and atmospheric pressure with refined sunflower oil (acid value of 0.07 mg KOH/g; density of 0.902 g/mL at 323 K). Selected reaction temperature is not far from the boiling point of methanol (337.8 K at 1 atm) while allowing to limit methanol losses that could otherwise affect the reactants and products concentrations and then the kinetics. HPLC grade methanol and NaOH (Sigma Aldrich, St. Louis, MO, USA, 99.9% P.A.) were used. Experiments were performed at constant reaction volume, with the proper amount of the reactants as to get initial methanol:oil molar ratios between 6:1 and 24:1. The catalyst (NaOH) concentration was within the 0.1–0.4 wt.% range, referred to the oil weight. Two kinds of samples were withdrawn at various intervals during the course of the reaction: samples for global composition analysis and samples for liquid phases characterization. As for the kinetic study, a first set of samples (around 0.5 g each) were taken using a recirculation loop comprising a diaphragm-type metering pump and a stainless-steel three-way ball valve. This liquid was poured upon 0.1 mL of a 0.6 N acetic acid solution in tetrahydrofuran for catalyst neutralization (THF, Scharlau, Barcelona, Spain, HPLC grade), and immediately diluted in 10 g of additional THF. As concerns the samples for liquid phases analysis, about 40 g of the reaction medium were directly taken through one of the reactor’s ports and immediately centrifuged at 6000 rpm for 15 min in order to separate the liquid phases and quench the reaction. The two liquid phases obtained were further diluted in THF. This process was typically repeated on about 7–8 samples for each reaction carried out.

2.2. Analytical Methods

Samples dissolved in THF were filtered on a 0.22 μm Teflon filter and analyzed by means of Size Exclusion Chromatography (SEC) for determining their composition. This technique allowed the determination of the reactants, both methanol and triglycerides from sunflower oil, as well as the main products, namely, biodiesel and glycerol, and the intermediate diglycerides and monoglycerides. The analysis system consisted of a Waters 510 HPLC pump, a Rheodyne 7725i manual injector and a Waters model 410 differential refractive index (RI) detector [23]. The mobile phase was HPLC grade THF at 1.0 mL/min flow rate. A configuration of three GPC columns connected in series was used: Styragel^® HR0.5, HR1 and HR2. Sample injection volume was 50 μL, and all the analyses were carried out at room temperature. The experimental errors were ±2% for the oil conversion and FAMEs yield, ±4% and ±6% for diglycerides and monoglycerides yield, respectively.

The miscibility of the oil-methanol-biodiesel system was studied by titrations with biodiesel of methanol-oil mixtures. The end point was established by turbidity analyses [22,24]. Briefly, stirring of methanol-oil mixtures generates a highly turbid emulsion; after that, biodiesel was added until a single transparent homogeneous phase was obtained. Biodiesel produced at the laboratory from refined sunflower oil was used in these experiments.

Density values were obtained by a tensiometer Krüss K100, using a DE-0601 adapter and viscosity was measured by a Viscostar Plus L. Densities of the pure components (oil, methanol, glycerol and biodiesel) were measured both at room temperature and at 323 K.

3. Kinetic Modelling

3.1. Reaction System

Vegetable oils methanolysis comprises three consecutive reversible steps represented by Equations (1)–(3) [25]. In the first one, triglycerides (TG) conversion into diglycerides (DG) occurs, followed by the conversion of diglycerides into monoglycerides (MG). The last step permits to produce glycerol (GL) from MG. In each of these reactions, a molecule of methanol (MeOH) is consumed whereas a molecule of fatty acid methyl ester (ME), or biodiesel, is formed.

$$\mathrm{TG}+\mathrm{MeOH} \begin{array}{c}{\mathrm{k}}_1\\ \rightleftarrows \\ {\mathrm{k}}_{-1}\end{array} \mathrm{DG}+\mathrm{ME}$$

(1)

$$\mathrm{DG}+\mathrm{MeOH} \begin{array}{c}{\mathrm{k}}_2\\ \rightleftarrows \\ {\mathrm{k}}_{-2}\end{array} \mathrm{MG}+\mathrm{ME}$$

(2)

$$\mathrm{MG}+\mathrm{MeOH} \begin{array}{c}{\mathrm{k}}_3\\ \rightleftarrows \\ {\mathrm{k}}_{-3}\end{array} \mathrm{GL}+\mathrm{ME}$$

(3)

According to the overall reaction stoichiometry, three methanol molecules are necessary for the complete conversion of a TG molecule to obtain three ME molecules and GL. Some kinetic studies take into account the overall conversion as a single reaction, the so-called Shunt reaction [2,25].

When alkaline hydroxides are used as catalysts, it is necessary to consider the saponification reactions that lead to the irreversible formation of soaps and then to the decrease of the catalyst concentration, which is obviously undesirable. Saponification reactions can involve free fatty acids (FFAs) and acylglycerols as represented by Equations (4)–(7).

$$\mathrm{FAA}+\mathrm{NaOH} \stackrel{{ \mathrm{k}}_{\mathrm{sap} }}{\to }{ \mathrm{RCOONa}}_{}+{\mathrm{H}}_2\mathrm{O}$$

(4)

$$\mathrm{TG}+\mathrm{NaOH} \stackrel{{ \mathrm{k}}_{\mathrm{sap} }}{\to }{ \mathrm{RCOONa}}_{}+\mathrm{DG}$$

(5)

$$\mathrm{DG}+\mathrm{NaOH} \stackrel{{ \mathrm{k}}_{\mathrm{sap} }}{\to }{ \mathrm{RCOONa}}_{}+\mathrm{MG}$$

(6)

$$\mathrm{MG}+\mathrm{NaOH} \stackrel{{ \mathrm{k}}_{\mathrm{sap} }}{\to }{ \mathrm{RCOONa}}_{}+\mathrm{GL}$$

(7)

It can be seen that the methanolysis of triglycerides constitutes a complex reaction system in which various chemical reactions can take place simultaneously. In addition, two immiscible liquid phases coexist that will contain different concentrations of the several chemical species involved in the process according to their affinity towards the oil and the alcohol. The distribution of several compounds, including the catalyst, among the phases was investigated by Zhou and Boocock [26]. Chiu et al. [27] used a MeOH distribution coefficient, defined as the ratio of the mass fraction in each phase. The distribution of the components among the phases in mixtures of biodiesel, glycerol, methanol and water has been also reported [28]. All of these studies were concerned with the reactants and final products. As for the intermediate products, Devender et al. stated that presence of monoglycerides greatly affected the compounds partitions [29].

3.2. Kinetic Models Formulation

Two kinetic models have been developed. The simplest one is the pseudo-homogeneous model, so named because no distinction among phases is made. In this case, it is considered that the reactions take place in the entire reactor volume, which behaves like a homogeneous system. Other models can be elaborated with similar complexity level by considering that the transesterification reactions occur exclusively in one of the liquid phases (lipidic or alcoholic).

A more realistic approach is to consider that the reactions take place at the liquid–liquid interface. However, this requires knowing the interfacial area, which is directly related to the droplet size distribution of the dispersed phase, that is, the alcoholic one. However, the drops size depends on the mixing conditions and undergoes variations throughout the reaction that cannot been easily measured. To overcome this difficulty retaining the heterogeneous character of the reacting system, an alternative model is proposed in this work considering that the overall reaction is the result of the contributions from each individual phase. It has been assumed that the interfacial mass transfer is not a limiting factor provided that an efficient mixing system is used. In addition, it has been considered that the distribution of the several compounds between the two phases is fast enough as to reach partitioning equilibrium, so that the concentrations can be related through the corresponding partition coefficients.

In this work, the partition coefficient of each individual compound ($\mathrm{K}\left(\mathrm{i}\right)$) involved in the reaction is defined as the ratio between the mass fractions in the lipidic (${\mathrm{w}}_{\mathrm{Lip}}\left(\mathrm{i}\right)$) and alcoholic (${\mathrm{w}}_{\mathrm{Ol}}\left(\mathrm{i}\right)$) phases, so being ${\mathrm{m}}_{\mathrm{j}}\left(\mathrm{i}\right)$ the mass of species i the phase j:

$$\mathrm{K}\left(\mathrm{i}\right)=\frac{{\mathrm{w}}_{\mathrm{Lip}}\left(\mathrm{i}\right)}{{\mathrm{w}}_{\mathrm{Ol}}\left(\mathrm{i}\right)}=\frac{\raisebox{1ex}{${\mathrm{m}}_{\mathrm{Lip}}\left(\mathrm{i}\right)$}\!\left/ \!\raisebox{-1ex}{$\sum {\mathrm{m}}_{\mathrm{Lip}}\left(\mathrm{i}\right)$}\right.}{\raisebox{1ex}{${\mathrm{m}}_{\mathrm{Ol}}\left(\mathrm{i}\right)$}\!\left/ \!\raisebox{-1ex}{$\sum {\mathrm{m}}_{\mathrm{Ol}}\left(\mathrm{i}\right)$}\right.}$$

(8)

These coefficients, like the kinetic constants, depend on the temperature and should be determined at the prevailing reaction conditions.

The general expression for the mass balance of each individual species i can be written as follows:

$$\frac{{\mathrm{dN}}_{\mathrm{i}}}{\mathrm{dt}}={\displaystyle \sum }_{\mathrm{j}=1}^2{\mathsf{\nu }}_{\mathrm{ik}}·{\mathrm{R}}_{\mathrm{ik}}·{\mathrm{V}}_{\mathrm{j}}$$

(9)

In this equation, ${\mathrm{N}}_{\mathrm{i}}$ stands for the number of moles of species i; ${\mathsf{\nu }}_{\mathrm{ik}}$ and ${\mathrm{R}}_{\mathrm{ik}}$ denote the stoichiometric coefficient and reaction rate, respectively, of species i in reaction k. A term adopts negative values when i is a reactant and positive values when i is a reaction product. ${\mathrm{V}}_{\mathrm{j}}$ stands for the total reaction volume of the pseudo-homogeneous model, and the specific volume of each phase (j = 1 for the alcoholic phase and j = 2 for the lipid phase) in the heterogeneous model, which has been calculated following Equation (10):

$${\mathrm{V}}_{\mathrm{j}}={\displaystyle \sum }_{\mathrm{i}}\frac{{\mathrm{m}}_{\mathrm{j}}\left(\mathrm{i}\right)}{{\mathsf{\rho }}_{\mathrm{i}}}$$

(10)

Accordingly, the individual balances written in terms of the molar concentrations of each species [i] are as follows:

$$\begin{array}{ll}\raisebox{1ex}{$\mathrm{dTG}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.={\displaystyle \sum _{\mathrm{j}=1}^2}& (-{\mathrm{k}}_1·\left({\left[\mathrm{TG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}1}}·{\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}-{\mathrm{k}}_{\mathrm{sap}}·{\left[\mathrm{TG}\right]}_{\mathrm{j}}\\ & ·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}·{\mathsf{\delta }}_{\mathrm{sap}})\end{array}$$

(11)

$$\begin{array}{cc}\hfill \raisebox{1ex}{$\mathrm{dDG}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.=& {\displaystyle \sum _{\mathrm{j}=1}^2}({\mathrm{k}}_1·\left({\left[\mathrm{TG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}1}}·{\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}-{\mathrm{k}}_2\hfill \\ & ·\left({\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}2}}·{\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}+{\mathrm{k}}_{\mathrm{sap}}\hfill \\ & ·({\left[\mathrm{TG}\right]}_{\mathrm{j}}-{\left[\mathrm{DG}\right]}_{\mathrm{j}})·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}·{\mathsf{\delta }}_{\mathrm{sap}})\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill \raisebox{1ex}{$\mathrm{dMG}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.=& {\displaystyle \sum _{\mathrm{j}=1}^2}({\mathrm{k}}_2·\left({\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}2}}·{\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}-{\mathrm{k}}_3\hfill \\ & ·\left({\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}3}}·{\left[\mathrm{GL}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}+{\mathrm{k}}_{\mathrm{sap}}\hfill \\ & ·({\left[\mathrm{DG}\right]}_{\mathrm{j}}-{\left[\mathrm{MG}\right]}_{\mathrm{j}})·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}·{\mathsf{\delta }}_{\mathrm{sap}})\hfill \end{array}$$

(13)

$$\begin{array}{cc}\hfill \raisebox{1ex}{$\mathrm{dME}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.={\displaystyle \sum _{\mathrm{j}=1}^2}& ({\mathrm{k}}_1·\left({\left[\mathrm{TG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}1}}·{\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}+{\mathrm{k}}_2\hfill \\ & ·\left({\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}2}}·{\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}+{\mathrm{k}}_3\hfill \\ & ·\left({\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}3}}·{\left[\mathrm{GL}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}})\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill \raisebox{1ex}{$\mathrm{dGL}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.={\displaystyle \sum _{\mathrm{j}=1}^2}& ({\mathrm{k}}_3·\left({\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}3}}·{\left[\mathrm{GL}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}+{\mathrm{k}}_{\mathrm{sap}}·{\left[\mathrm{MG}\right]}_{\mathrm{j}}\hfill \\ & ·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}·{\mathsf{\delta }}_{\mathrm{sap}})\hfill \end{array}$$

(15)

$$\begin{array}{cc}\hfill \raisebox{1ex}{$\mathrm{dMeOH}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.=& {\displaystyle \sum _{\mathrm{j}=1}^2}(-{\mathrm{k}}_1·\left({\left[\mathrm{TG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}1}}·{\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}-{\mathrm{k}}_2\hfill \\ & ·\left({\left[\mathrm{DG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}2}}·{\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}-{\mathrm{k}}_3\hfill \\ & ·\left({\left[\mathrm{MG}\right]}_{\mathrm{j}}·{\left[\mathrm{MeOH}\right]}_{\mathrm{j}}-\frac{1}{{\mathrm{K}}_{\mathrm{C}3}}·{\left[\mathrm{GL}\right]}_{\mathrm{j}}·{\left[\mathrm{ME}\right]}_{\mathrm{j}}\right)·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}} )\hfill \end{array}$$

(16)

$$\raisebox{1ex}{$\mathrm{dNaOH}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.={\displaystyle \sum }_{\mathrm{j}=1}^2\left(-{\mathrm{k}}_{\mathrm{sap}}·({\left[\mathrm{TG}\right]}_{\mathrm{j}}+{\left[\mathrm{DG}\right]}_{\mathrm{j}}+{\left[\mathrm{MG}\right]}_{\mathrm{j}})·{\left[\mathrm{NaOH}\right]}_{\mathrm{j}}·{\mathrm{V}}_{\mathrm{j}}\right)·{\mathsf{\delta }}_{\mathrm{sap}}$$

(17)

Operator ${\mathsf{\delta }}_{\mathrm{sap}}$ allows considering (${\mathsf{\delta }}_{\mathrm{sap}}$ = 1) or not (${\mathsf{\delta }}_{\mathrm{sap}}$ = 0) the saponification reactions.

In this work, refined sunflower oil of very low acidity (0.07 mg KOH/g) has been used; it is assumed that FFA neutralization takes place instantaneously. In order to reduce the number of adjustable parameters, the same kinetic constant (${\mathrm{k}}_{\mathrm{sap}}$) is applied to all saponification reactions, among which, TG saponification obviously stands out due to predominance of that reactant. Equilibrium constants (${\mathrm{K}}_{\mathrm{C}1}$, ${\mathrm{K}}_{\mathrm{C}2}$, ${\mathrm{K}}_{\mathrm{C}3}$) were taken from the literature [4], and the partition coefficients ($\mathrm{K}\left(\mathrm{i}\right)$) have been considered constant throughout the reaction.

3.3. Kinetic Parameters Estimation

To integrate Equations (11)–(17), it is necessary to calculate the concentration of species i in both phases. Thus, the following iterative method has been implemented [30]: suppose the mass of the lipidic (${\mathrm{m}}_{\mathrm{Lip},\mathrm{s}}$) and alcoholic (${\mathrm{m}}_{\mathrm{Ol},\mathrm{s}}$) phases are known; then, Equations (18) and (19) can be solved for each species i. Furthermore, Equations (20) and (21) provide new values for each phase masses, and the condition given by Equation (22) is applied, giving rise to a new iteration or the converged solution of the problem.

$${\mathrm{m}}_{\mathrm{Lip}}\left(\mathrm{i}\right)=\frac{\mathrm{m} \left(\mathrm{i}\right)}{1+\frac{{\mathrm{m}}_{\mathrm{Ol},\mathrm{s}}}{{\mathrm{m}}_{\mathrm{Lip},\mathrm{s}}·\mathrm{K}\left(\mathrm{i}\right)}}$$

(18)

$${\mathrm{m}}_{\mathrm{Ol}}\left(\mathrm{i}\right)=\mathrm{m}\left(\mathrm{i}\right)-{\mathrm{m}}_{\mathrm{Lip}}\left(\mathrm{i}\right)$$

(19)

$${\mathrm{m}}_{\mathrm{Lip}}={\displaystyle \sum }_{\mathrm{i}}{\mathrm{m}}_{\mathrm{Lip}}\left(\mathrm{i}\right)$$

(20)

$${\mathrm{m}}_{\mathrm{Ol}}={\displaystyle \sum }_{\mathrm{i}}{\mathrm{m}}_{\mathrm{Ol}}\left(\mathrm{i}\right)$$

(21)

$$\mathrm{if} \left|\frac{{\mathrm{m}}_{\mathrm{Ol},\mathrm{s}}-{\mathrm{m}}_{\mathrm{Ol}}}{{\mathrm{m}}_{\mathrm{Ol},\mathrm{s}}}\right| \mathrm{or} \left|\frac{{\mathrm{m}}_{\mathrm{Ol},\mathrm{s}}-{\mathrm{m}}_{\mathrm{Ol}}}{{\mathrm{m}}_{\mathrm{Ol},\mathrm{s}}}\right|>{\mathsf{\epsilon } \mathrm{then} \mathrm{m}}_{\mathrm{Ol},\mathrm{s}}={\mathrm{m}}_{\mathrm{Ol}} ;{ \mathrm{m}}_{\mathrm{Lip},\mathrm{s}}={\mathrm{m}}_{\mathrm{Lip}}$$

(22)

Integration of Equations (11)–(17) was carried out using Compaq Visual Fortran V6.6 that includes the IMSL subroutines library (IMSL Fortran 90 MP Library). The ordinary differential equations were solved using the DIVPRK subroutine that applies the 5th and 6th orders adaptive Runge-Kutta-Verner method. This allowed obtaining the number of moles of each compound at a given reaction time (${\mathrm{N}}_{\mathrm{i}}$) and, given the initial number of moles of oil (${\mathrm{N}}_{\mathrm{TG}0}$), calculating the oil conversion (${\mathrm{X}}_{\mathrm{TG}}$, Equation (23)) and the yields of biodiesel (${\mathrm{Y}}_{\mathrm{BD}}$, Equation (24)), diglycerides, monoglycerides and glycerol (${\mathrm{Y}}_{\mathrm{DG}}$, ${\mathrm{Y}}_{\mathrm{MG}}$, and ${\mathrm{Y}}_{\mathrm{GL}}$, respectively, Equation (25)):

$${\mathrm{X}}_{\mathrm{TG}}=\frac{{\mathrm{N}}_{\mathrm{TG}0}-{\mathrm{N}}_{\mathrm{TG}}}{{\mathrm{N}}_{\mathrm{TG}0}}$$

(23)

$${\mathrm{Y}}_{\mathrm{BD}}=\frac{{\mathrm{N}}_{\mathrm{BD}}}{3·{\mathrm{N}}_{\mathrm{TG}0}}$$

(24)

$${\mathrm{Y}}_{\mathrm{i}}=\frac{{\mathrm{N}}_{\mathrm{i}}}{{\mathrm{N}}_{\mathrm{i}0}};\hspace{1em}\left(\mathrm{i}=\mathrm{DG}, \mathrm{MG}, \mathrm{GL}\right)$$

(25)

Kinetic parameters estimation was carried out using a complex method of direct search furnished by the DBCPOL optimization subroutine (IMSL Fortran 90 MP Library). An objective function given by the sum of the differences between the experimental values of the oil conversion and products yields and the values predicted by the model was established as follows:

$$\mathrm{F}={\displaystyle \sum }_{\mathrm{i}=1, \mathrm{NR}}{\displaystyle \sum }_{\mathrm{j}=1, \mathrm{N}}{\displaystyle \sum }_{\mathrm{k}=1, \mathrm{NV}}\frac{1}{\mathrm{N}}·\frac{{\left({\mathrm{V}}_{\mathrm{mod}}\left(\mathrm{i}, \mathrm{j}, \mathrm{k}\right)-{\mathrm{V}}_{\exp }\left(\mathrm{i}, \mathrm{j}, \mathrm{k}\right)\right)}^2}{{\mathrm{V}}_{\exp }{\left(\mathrm{i}, \mathrm{j}, \mathrm{k}\right)}^2}·{\mathrm{w}}_{\mathrm{k}}$$

(26)

where ${\mathrm{V}}_{\exp }\left(\mathrm{i}, \mathrm{j}, \mathrm{k}\right)$ and ${\mathrm{V}}_{\mathrm{mod}}\left(\mathrm{i}, \mathrm{j}, \mathrm{k}\right)$ are the values of the k different variables obtained experimentally and throughout the model, respectively, at each sampling point (j) for a given reaction (i). NR stands for the number of experimental runs performed, NV for the number of variables, and N for the number of experimental data for each run. Finally, ${\mathrm{w}}_{\mathrm{k}}$ is a weighting factor for each variable introduced to balance its contribution to the total error.

The mass balance equations can be rewritten as functions of the oil conversion instead of the reaction time applying the following relationships [16]:

$$\frac{{\mathrm{dN}}_{\mathrm{i}}}{\mathrm{dX}}=\frac{\raisebox{1ex}{${\mathrm{dN}}_{\mathrm{i}}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.}{\raisebox{1ex}{$\mathrm{dX}$}\!\left/ \!\raisebox{-1ex}{$\mathrm{dt}$}\right.}$$

(27)

$$\frac{\mathrm{dX}}{\mathrm{dt}}=-\frac{1}{{\mathrm{TG}}_0}·\frac{\mathrm{dTG}}{\mathrm{dt}}$$

(28)

Therefore, it is possible to estimate the ratios ${\mathrm{k}}_2$/${\mathrm{k}}_1$ and ${\mathrm{k}}_3$/${\mathrm{k}}_1$, which are used to estimate the kinetic constant ${\mathrm{k}}_1$ and the saponification constant, ${\mathrm{k}}_{\mathrm{sap}}$.

4. Results and Discussion

4.1. Relevance of the Heterogeneous Character of the Reaction System

As mentioned before, the reacting medium is a biphasic liquid–liquid system. This feature has obvious consequences as concerns the reaction rate because of the possible influence of mass transport limitations. Once the reaction starts, methanol and TG contents of the alcoholic and lipidic phases begin to decrease whereas those of ME and GL increase. These are the main components of the system which have different mutual affinities. In this study was observed that whereas the ME-MeOH, ME-TG and MeOH-GL binary mixtures are completely miscible, the GL-ME and GL-TG ones are virtually immiscible. At the reaction start, the alcoholic and the lipidic phases do not contain a single reactant because methanol and sunflower oil are not fully immiscible. It has been determined that mixtures of methanol and the oil generate an alcoholic phase containing 2.8 wt.% oil and a lipidic phase containing 5.9 wt.% of methanol at room temperature.

Previous studies [29,31] have focused on the behavior of GL-MeOH-ME ternary mixtures because they contain the most abundant compounds at very high oil conversions. This mixture is of interest as concerns the separation processes required to isolate and purify biodiesel. However, during the course of the reaction, it is also interesting to know the behavior of mixtures containing oil (TG). In this regard, Figure 1 shows the equilibrium diagram with the solubility curve of the ternary TG-ME-MeOH system.

It can be seen that any mixture composition located under the curve will give rise to two phases. In contrast, compositions above the curve lead to a single phase, which corresponds to mixtures very rich in biodiesel. Interestingly, this means that biodiesel can be regarded as a co-solvent of the methanolysis reactants. The curve is closer to the axis on the MeOH side; this is because the solubility of the oil in methanol is much lower than that of the alcohol in the oil.

Besides the main components mentioned before, reaction intermediates such as monoglycerides and diglycerides should be taken also into account as they can reach at intermediate oil conversion levels non-negligible contents of c.a. 10 wt.%. Figure 2 shows the experimental results for the phases’ composition evolution with TG conversion.

For an oil conversion of 0.5, TG content in the lipidic phase is about 30 wt.% whereas those of ME and MeOH reach c.a. 50 wt.% and roughly 10 wt.%, respectively (Figure 2). Among the intermediate acylglycerols, DG predominate over MG, as expected for a reaction in series at intermediate conversions of the main reactant. Glycerol content of the lipidic phase is extremely low, though it reaches about 2 wt.% at very high TG conversions. As for the alcoholic phase, clearly MeOH is the most abundant compound throughout the reaction due to the very high methanol excess used compared to the stoichiometric value (MeOH/oil ratio of 3). Glycerol content is of course relevant (about 16 wt.% at 0.5 oil conversion) and that of ME reaches c.a. 6 wt.%. It can be concluded that the lipidic phase is much more affected than the alcoholic one by the presence of reaction products. As a matter of fact, the presence of triglycerides and acylglycerols in the alcoholic phase is almost negligible.

Partition coefficients were obtained experimentally from the analysis of samples with different oil conversion levels corresponding to methanolysis reactions conducted at 323 K and methanol/oil molar ratios between 6 and 12. These coefficients were estimated from the analysis of reaction samples because it is necessary to take into account the presence of intermediate products (diglycerides, and specially monoglycerides), which have emulsifying properties and can affect the distribution of products. Results are shown in Figure 3a,b for the main and minority compounds, respectively. According to the definition given by Equation (8), the partition coefficients of the main components in the lipidic phase have values higher than one, whereas the most abundant in the alcoholic phase have values within the 0–1 range. It was found that the oil conversion level does not affect the partition coefficient values for most compounds. For TG, DG and MG, the following mean values have been adopted: 92, 13.5 and 0.67, respectively. It is apparent that the values decrease as the polarity and affinity towards the alcohol increase.

It should be noted that at the beginning of the reactions (very low oil conversions), the experimental error is higher due to the low concentration of the products. ME partition coefficient is slightly affected by the excess of methanol in the reaction mixture. A linear relationship can be established between the partition coefficient and the molar methanol/oil ratio; as the MeOH in excess increases, the biodiesel partition coefficient decreases. It is noteworthy the very low values of the mass partition coefficients for MeOH and specially GL, whereas for MG it is close to 1, meaning that monoglycerides show similar affinities for both phases. This is the reason by which monoglycerides are capable of stabilizing the emulsion and are undesired in the final mixture as they complicate the separation of the lipidic and alcoholic phases. In this work, the catalyst is previously dissolved in the methanol, so it will be supposed that both MeOH and NaOH have the same partition coefficient.

The change of composition of the liquid phases during the reaction is also evidenced by significant changes in some of their physico-chemical properties. It has been found that both phases experience changes in their viscosity. In this regard, the presence of minority components such as diglycerides and monoglycerides can have a significant influence on the liquid phases’ viscosity. Furthermore, there are evident changes in the density. As for the reactants, the oil is significantly heavier (0.902 g/mL at 323 K) than methanol (0.764 g/mL at 323 K). As the methanolysis reaction progresses, the lipidic phase composition evolves increasing the biodiesel content at the expense of oil, which results in a decrease of the density due to the lower density of biodiesel (0.863 g/mL at 323 K) compared to that of the oil. At the same time, the alcoholic phase becomes heavier because methanol is being replaced by glycerol which has a much higher density (1.244 g/mL at 323 K). Consequently, an oil conversion exists at which the alcoholic phase becomes heavier than the lipidic one. Changes in viscosity and density can affect relevant emulsion properties such as the interfacial surface area. Previous works on immiscible liquids behavior [32] demonstrated that inversion between the dispersed and continuous phases can take place aided by enough mixing time and the presence of suitable surfactants. The general trend is that the higher is the continuous phase volume, the longer the time required to attain the inversion. As for the system of interest in this study, the oil conversion at which the alcoholic phase becomes the heaviest one depends on the initial methanol/oil molar ratio. When the ratio was 6:1, the change was observed a few minutes after the reaction has started. On the other hand, when the reaction is carried out with a larger excess of methanol, the change requires longer reaction times. It can be expected that such changes in the liquid phases’ properties will have an important impact on the emulsion behavior. Diglycerides and monoglycerides seem to play a relevant role due to their surfactant properties [33], as well as the mixing conditions. The nature of the dispersed and that of the continuous phases depend in a complex way on the initial composition of the system and the operating conditions [32].

4.2. Kinetic Models Performance

Previous pseudo-homogeneous [4,34] or hybrid approaches [7] considered only a single methanol/oil ratio of 6:1. These models failed in describing the evolution of the reaction system upon changing that ratio. Vicente et al. [34] considered that the conversion of monoglycerides into glycerol is an irreversible reaction. In our case, the kinetic constant, though low, is not so small as to neglect the reverse step.

The best fits for ${\mathrm{k}}_1$ were 0.014 L²/(mol²·s) for the pseudo-homogeneous model and 0.19 L²/(mol²·s) for the heterogeneous model. In both cases, wide ranges of methanol/oil molar ratios (6:1 to 24:1) and NaOH concentrations (0.005 to 0.5 wt.%) were considered. It should be noted that, according, to the pseudo-homogeneous model, the concentrations of the several compounds become diluted because the whole volume of the liquid system is accessible to reactants and products. In contrast, as for the heterogeneous model, the concentration of methanol in the lipidic phase and that of oil in the alcoholic phase are very low due to poor mutual solubility. These distinct features of the models can explain the large difference (almost an order of magnitude) between the values of the kinetic constants. The same happens to ${\mathrm{k}}_{\mathrm{sap}}$, whose values are 0.00027 and 0.0028 L/(mol·s) for the pseudo-homogeneous and the heterogeneous models, respectively.

Figure 4 shows the results of the reactions carried out at 323 K, with 0.2 wt.% of NaOH and methanol/oil molar ratios of 6:1, 9:1, 12:1, and 24:1. The Figure also includes the best fittings of the pseudo-homogeneous and heterogeneous kinetic models to the experimental results. The pseudo-homogeneous model (Figure 4a) overestimates the oil conversions for the highest methanol/oil ratios. Even for the lowest methanol/oil ratio, the model does not capture the evolution of the system at short reaction times. In contrast, the heterogeneous model (Figure 4b) provides a better description of the system evolution. Particularly, the model adapts well to the range of oil conversions achieved with each methanol/oil ratio for long reaction times.

In previous works, the pseudo-homogeneous model has been applied to cases in which the initial methanol/oil ratio remained constant. For the sake of comparison, in this study, parameters estimation was also been carried for each methanol/oil molar ratio separately. The results show a clear trend towards lower ${\mathrm{k}}_1$ values as the excess of methanol increased, i.e., ${\mathrm{k}}_1$ was 0.030, 0.013, 0.008 and 0.006 L²/(mol²·s) for initial methanol/oil molar ratios of 6, 9, 12 and 24:1, respectively. When all the results are considered, a somewhat intermediate value of 0.014 L²/(mol²·s) is obtained, as indicated before. The same effect is observed on the saponification kinetic constant. As far as only one set of results corresponding to the same methanol/oi ratio is considered, the quality of the individual fit improves notably, but the model loses general validity and prediction capability.

From the foregoing discussion, it can be concluded that the pseudo-homogeneous model is only capable of describing the methanolysis of triglycerides if a set of specific kinetic parameters is used for each methanol/oil molar ratio. This model is not capable of capturing the effects associated to the use of different methanol/oil ratios. The proportion between these two immiscible reactants is a key parameter arising from the heterogeneous character of the methanolysis reacting system, so a heterogeneous model is required to describe it. In this regard, the heterogeneous model provides the following ratios between the kinetic constants at 323 K: ${\mathrm{k}}_2$/${\mathrm{k}}_1$ = 0.979 and ${\mathrm{k}}_3$/${\mathrm{k}}_1$ = 0.12. It is clear that the second and specially the third transesterification steps are slower than the first one.

Figure 5 shows the evolution of the contents of the several components of the reaction mixture with reaction time at 323 K with 0.2 wt.% NaOH and methanol/oil molar ratios of 6:1 (Figure 5a) and 12:1 (Figure 5b). It can be seen that the heterogeneous model fits well the results for the most abundant compounds, triglycerides, methyl esters and glycerol under different conditions of methanol excess. It also captures the maxima that characterize the evolution with reaction time of the diglycerides and monoglycerides contents.

The combined effects of the methanol/oil molar ratio and the oil conversion on the evolution of the intermediate compounds contents according to the heterogeneous kinetic model are shown in Figure 6. The two maxima are clearly seen. Both become more pronounced as the methanol/oil ratio decreases, i.e., values between 9:1 and 4:1, likely due to the slowing down of the second and third methanolysis steps. The yields are almost insensitive to the methanol excess for methanol/oil ratios above 12:1.

Another advantage of the heterogeneous kinetic model is the possibility of predicting the evolution of the composition of each phase. Figure 7 shows the results corresponding to the methanol and biodiesel contents in the lipidic and alcoholic phases for reactions performed at 323 K with 0.1 wt.% NaOH and different initial methanol/oil molar ratios. It can be seen that the model describes very well the biodiesel content. It is remarkable the good description provided of the very low contents of methanol in the lipidic phase and that of biodiesel in the alcoholic one. It is observed that the methanol/oil ratio has little effect on the biodiesel content of the lipidic phase. The non-negligible solubility of methanol in the lipidic phase, particularly if biodiesel is present, impedes to reach methyl esters contents closer to 100 %. Obviously, the fact that the reactions are reversible and the presence of intermediate transesterification products also contribute to the maximum biodiesel content in the lipidic phase. In this case, it should be taken into account also that the catalyst content is very low, which makes difficult reaching higher biodiesel contents due to the negative influence of the saponification reactions.

As for the decrease of the methanol content in the alcoholic phase, the model captures well the trend associated to the effect of the methanol/oil ratio although there is an important scatter of the experimental results that affects the fitting quality. The experimental values of methanol content in the lipidic phase and that of biodiesel content in the alcoholic phase are very similar for methanol/oil ratios between 6:1 and 12:1. This result supports the assumption made that the partition coefficients are independent of the methanol/oil ratio.

Figure 8 shows the experimental and kinetic modelling results corresponding to the effect of the initial catalyst concentration. As can be seen, the pseudo-homogeneous model overestimates the oil conversion, especially at low catalyst contents. Therefore, just considering the saponification reactions is not sufficient to obtain a good description of the reacting system evolution. The heterogeneous kinetic model, on the other hand, fits well the results obtained at catalyst concentrations above 0.1 wt.% NaOH; however, for the lowest catalyst concentration, the discrepancies remain important. Clearly, the model provides results as if the catalyst concentration was higher than expected. In this regard, one of our assumptions was to assign the same kinetic constant to all the saponification reactions in order to reduce the number of parameters. However, it seems that this simplification has to be reconsidered. In addition, catalyst deactivation can be also an issue, as previously noted [35].

5. Conclusions

In this work, triglyceride methanolysis with a homogeneous catalyst has been studied experimentally and by means of kinetic modelling. Two levels of detail have been considered in the description of the reaction medium nature as concerns its biphasic character.

The complex behavior of the reaction system has obvious consequences as regards the methanolysis kinetics. For this reason, a pseudo-homogeneous model is not able to describe adequately the process, since it disregards the biphasic character of the reaction medium. This simple model, however, can be of some utility if the initial methanol/oil molar ratio is not changed. However, this parameter is of high practical importance due to the reversible nature of the three steps involved in the methanolysis process. Conversely, a heterogeneous model is capable of capturing the effects of the main reaction parameters on the oil conversion and the composition of the phases because it considers the presence of the two liquid phases, the distribution of the several compounds between them, and saponification reactions. Good results have been obtained for the main and minority components. The quality of the fittings is lower at the lowest catalyst concentrations considered. Further adjustment is required, for example, considering different saponification kinetic constants and catalyst deactivation phenomena.