Mathematical Modeling of the Transesterification Reaction by Finite Elements: Optimization of Kinetic Parameters Using the Simplex Sequential Method

Abstract

The use of biodiesel as fuel has been presented as a viable alternative in the search for renewable energies. It can be produced from the transesterification reaction of vegetable oils with methanol in the presence of sodium or potassium hydroxide as catalysts. In the present research, the transesterification reaction of soybean oil was modeled considering the three steps of the direct and reverse reactions following a second order general kinetics by the finite element method using the COMSOL Multiphysics^® software. The values of rate constants were determined using the simplex optimization method coupled with the desirability functions. The optimized rate constants for the forward reactions were 0.250 × 10⁻⁶ (k_1f), 1.137 × 10⁻⁶ (k_2f) and 3.134 × 10⁻⁶ (k_3f); and for the reverse reactions were 0.202 × 10⁻⁶ (k_1r), 0.884 × 10⁻⁶ (k_2r) and 0.219 × 10⁻⁶ (k_3r) all in m³mol⁻¹s⁻¹. The kinetic model proposed for the reaction can be simulated by the finite element method (FEM) under realistic conditions.

Keywords:
biodiesel; simulation; rate constants

Introduction

Biodiesel is classified as a fuel obtained from renewable sources and its production has as main objective the gradual replacement of petroleum diesel, encouraging the research and production of this biofuel.¹ It is one of the main biofuels used in the global transport fleet with increasing consumption in the last decade. United States, Brazil and some European countries are mainly responsible for this growing market share.²

For the industrial production of biodiesel, kinetic parameters and reaction rate constants must be established, as they are fundamental for the design of production reactors. Regardless of the type of oil, reaction temperature and the molar ratio of methanol/oil, the mechanism involves three consecutive reversible reactions.³ Fatty acid methyl esters (FAME) and glycerol (GL) are products of the transesterification of triglyceride (TG) of fats and oils with alcohols as methanol (MeOH) or ethanol (EtOH) in the presence of an acid or an alkaline catalyst. Three consecutive stepwise reversible reactions occur with intermediate formation of diglycerides (DG) and monoglycerides (MG).⁴

Although the kinetic model remains unchanged regardless of the feedstock, temperature or molar ratio oil and alcohol, the values of the rate constants change, and all these variables are evaluated at the same time, thus the simulation becomes it more viable. The multiphysics platform simulation based on the finite element method (FEM) is an analytical methodology that has been highlighted in science areas.^5,6

FEM is used to solve differential equations using a variational formulation or a residual weight procedure, where a partial differential equation is transformed into a system of ordinary differential equations, when the problem is time dependent. Several authors have used the simulation by finite elements, coupled to optimization methods and desirability functions, to determine parameters such as diffusion coefficients, laminar and Biot number.^5,7,8

One of the most used optimization methods is the simplex which is a recurring process that tends to bring the response to an optimal value through the reflection of specific points. Once in the vicinity of the optimal value, the simplex can undergo contraction in order to determine a more accurate position.^9-11

A variation of the simplex optimization method, called modified simplex, was developed by Nelder and Mead.¹² It can change size and shape, adjusting better in the search for optimum response. In order to improve and bring the method quickly, the super modified simplex method was developed.^10,11,13 In this algorithm, the simplex can be submitted to five operations: reflection, expansion, contraction, contraction with change of direction and massive contraction, whose purpose is to improve and bring more quickness to the method.¹⁰

The objective of this study was to model and simulate transesterification reaction using the finite element method associated with the simplex optimization to determine the forward and reverse rate constants in the three-step consecutive reversible reactions.

Experimental

Reagents and materials

For transesterification reactions was used commercial soybean oil (Soya, Brazil) purchased from the local market. All the reagents (methanol (Anidrol, Brazil), sodium hydroxide (Fmaia, Brazil), acetic acid (Ecibra, Brazil) and anhydrous sodium sulfate (Anidrol, Brazil)) are analytical grade. Deuterated chloroform containing 0.05% (v/v) tetramethylsilane (TMS, Sigma Aldrich, USA) was used for high field proton nuclear magnetic resonance (HF-¹H NMR) analysis.

Transesterification reaction

Transesterification reactions were carried out by homogeneous basic catalysis mixing 500 g of soybean oil, 140 mL of methanol (1:6 molar ratio of oil:methanol) and 0.75% m/m sodium hydroxide as catalyst, previously dissolved in methanol. The reaction was conducted with mechanical stirring maintained at 150 rpm for 62 min and temperature controlled at 20 ± 2 °C. All the reactions were carried out three times.

Sampling of transesterification reactions with HF- ¹H NMR

During the development of the reaction, were acquired 16 aliquots of 7.0 mL of the reaction medium for each reaction in the times of 3, 7, 10, 14, 18, 22, 26, 30, 34, 38, 42, 46, 50, 54, 58 and 62 min for analysis in the HF-¹H NMR.

Each aliquot was added to a tube containing glacial acetic acid to neutralize the catalyst. After homogenization, the mixture was washed with distilled water and centrifuged for three times. Always recovering the organic phase and discarding the aqueous phase. The organic phase was dried on Na₂SO₄ and 50 µL of organic phase were dissolved in 500 µL of deuterated chloroform (CDCl₃) for HF-¹H NMR spectra acquisition.

An NMR spectrometer (400 MHz, Bruker) was used. The acquisition parameters of the HF-¹H NMR spectra were single pulse, spectral width of 8012 Hz, relaxation delay of 1 s, 16 scans, acquisition time of 4.089 s and pulse width of 90°. The spectral regions studied were −2 to 10 ppm processed in the software Matlab^® 2016b (MathWorks, USA).¹⁴

Determination of conversion rates

The molar percentage conversion rates can be obtained by equations of glycerides (triglycerides, diglycerides and monoglycerides) and fatty acid methyl esters compounds present in the reaction medium during the transesterification reactions. Thus, the molar percentage conversion rates were determined based on the signal areas showed in the HF-¹H NMR spectra. The determination of the peak areas found in the HF-¹H NMR spectra was performed through a Matlab routine. The methanol and glycerol conversion rates were determined, through the material balance of the reaction, from the previously calculated concentrations and the initial concentrations of triglyceride and methanol.^15-17

Simulation and chemistry of the transesterification reaction

The simulation of reaction was performed using the COMSOL Multiphysics^® software¹⁸ and the default physical interface “Reaction Engineering (re)”. In this study, the three-step reversible reactions were forward and reverse reactions to follow second-order overall kinetics (equations 1-4).^4,19

Overall reaction:

(1) $\mathrm{TG}+3\mathrm{MeOH}\stackrel{\mathrm{NaOH}}{\rightleftharpoons }\mathrm{GL}+3\mathrm{FAME}$

Stepwise reactions:

(2) $\mathrm{TG}+\mathrm{MeOH}\underset{k_{1r}}{\overset{k_{1f}}{\rightleftharpoons }}\mathrm{DG}+\mathrm{FAME}$

(3) $\mathrm{DG}+\mathrm{MeOH}\underset{k_{2r}}{\overset{k_{2f}}{\rightleftharpoons }}\mathrm{MG}+\mathrm{FAME}$

(4) $\mathrm{MG}+\mathrm{MeOH}\underset{k_{3r}}{\overset{k_{3f}}{\rightleftharpoons }}\mathrm{GL}+\mathrm{FAME}$

The reaction scheme for methanolysis of vegetable oil is presented in equation 1. For overall reaction transesterification requires TG with MeOH, in the presence of an alkaline catalyst, yields FAME and GL. The stoichiometry of reaction requires three mol of MeOH and one mol of TG to give three mol of FAME and one mol of GL. For stepwise reactions MG and DG the intermediates are formed. Firstly, due to the TG and MeOH molecules reaction a DG and FAME are formed. Sequentially, this DG molecule reacts with a MeOH molecule and a MG molecule and another FAME molecule is formed. Finally, on the last step MG reacts with MeOH to form another FAME molecule and a GL molecule. That is, three FAME molecules on the whole reaction process.^4,19

The reaction steps are in equation 2-4, where k_nf are the forward rate constants and k_nr are the reverse rate constants. The three consecutive reversible reactions equations were defined under global definitions in COMSOL software;¹⁸ and all necessary physical values, like reaction rate constants, initial concentrations were described as parameters in all cases.

The three reversible reactions with rate constant expressions of the reaction (r_n), expressed in mol m⁻³ s⁻¹, are represented by equations 5-7.

(5) $r_{1f}=k_{1f}{c}_{\mathrm{TG}}{c}_{\mathrm{MeOH}}\mathrm{and}{r}_{1r}=k_{1r}{c}_{\mathrm{DG}}{c}_{\mathrm{FAME}}$

(6) $r_{2f}=k_{2f}{c}_{\mathrm{DG}}{c}_{\mathrm{MeOH}}\mathrm{and}{r}_{2r}=k_{2r}{c}_{\mathrm{MG}}{c}_{\mathrm{FAME}}$

(7) $r_{3f}=k_{3f}{c}_{\mathrm{MG}}{c}_{\mathrm{MeOH}}\mathrm{and}{r}_{3r}=k_{3r}{c}_{\mathrm{GL}}{c}_{\mathrm{FAME}}$

where f is forward, r is reverse rate constant expressions of the reaction, c is theoretical concentration of triglyceride (TG), diglyceride (DG), monoglyceride (MG), glycerol (GL), methanol (MeOH) and fatty acid methyl esters (FAME) for each time (min) for the HF-¹H NMR analysis.

In the transesterification reaction, vegetable oil and methanol are not miscible, so the reaction system consists of two phases at the initial stage. After this stage the mass transfer controls the kinetics of the reaction. As FAME are formed in the reaction system, they act as a system cosolvent solubilizing the vegetable oil and the methanol, thereby a single phase is formed.^4,15

As a boundary condition for the modeling of the kinetics of the transesterification reaction, it was considered the statements of Vicente et al.¹⁵ and Stamenkovic et al.:³

• (i) the initial stage of mass transfer control was negligible;

• (ii) the proportion of free fatty acid was negligible, and then the free fatty acid neutralization was not significant;

• (iii) the saponification reaction was insignificant;

• (iv) the values of the DG and MG molar conversions are represented by the sum of the 1,3-diglyceride (1,3-DG) and 1,2-diglyceride (1,2-DG) conversions; and 2-monoglyceride (2-MG) and 1-monoglyceride (1-MG), respectively.

The simulations were performed with a 0D dimension, considering a modeling relative tolerance error of 1 × 10⁻⁷. The algebraic equations solved by the backward differentiation formula (BDF), which is a time dependent solver, are the result of a discretization applied to the original differential equations that were solved by the finite element method (FEM) using the COMSOL Multiphysics^® interface.

Adjustment of the reaction rate constants

Reaction rate constants (k_n) for transesterification reactions, consists of three consecutive reversible reactions, where k_nf are the forward rate constants and k_nr are the reverse rate constants. The values of k_n were adjusted by the super modified simplex optimization method¹⁰ associated to the desirability functions.²⁰ These values have been assessed using the COMSOL Multiphysics^® software.¹⁸

Computer processing and program

The transesterification reaction was simulated using the COMSOL Multiphysics^® version 5.2 (COMSOL, Inc., Burlington, MA, USA).¹⁸ All results from the experiments were processed using a computer with Intel^® Core™ i7-4790 CPU^© 3.60GHz, 32GB RAM and 250GB HDD.

Results and Discussion

The molar conversion rate of the transesterification reaction was determined by equation 8, taking into account the conversion of TG into FAME molecules. This equation was adapted from the work published by Nieva-Echevarríaet al.,²¹ where they proposed different equations to quantify the molar percentage of mono-, di- and triglycerides in a complex lipid mixture by ¹H NMR. These equations were also adapted to be used on the determination of 1-MG (equation 9), 2-MG (equation 10), 1,2-DG (equation 11), TG (equation 12) and 1,3-DG (equation 13) present in the reaction medium during the transesterification reaction.

(8) $C_{\mathrm{FAME}}=0.67\times \frac{I_{\mathrm{FAME}}-\left(2\times \left(I_{1-\mathrm{MG}}+I_{1,2-\mathrm{DG}}\right)\right)}{I_{\mathrm{Acy}}}\times 100\%$

Equation 8 applies the area on the ¹H NMR spectrum from δ = 3.55 to δ = 3.78 ppm (I_FAME) to determine the molar percentage (conversion rate) of FAME in the reaction medium. In this region there is an overlapping of methoxylic protons of FAME molecules (3 protons, δ = 3.67 ppm), glyceryl group in 1,2-DG (I_1,2-DG, 2 protons at δ = 3.73 ppm) and half of the peak from the glyceryl group in 1-MG (I_1-MG, 1 proton at δ = 3.65 to δ = 3.73 ppm) and glyceryl group in 1-MG (I_1-MG, 2 protons at δ = 3.65 ppm) were subtracted from the I_FAME region.

The area of acyl groups (I_Acy, 2 protons at δ = 2.24 to δ = 2.40 ppm), present originally in TG molecules at beginning of the transesterification reaction, is used to relate the molar percentage of each compound to the initial amount of TG in the reaction medium, once these kind of protons are also present in MG, DG and FAME species.

Equations 9-11 were used for determination of molar percentage of 1-MG (C_1-MG), 2-MG (C_2-MG) and 1,2-DG (C_1,2-DG) in the reaction medium. For these types of protons there are no overlapping of any other protons in the considered region of integration. Thus, the equations take into account only the area of each specific protons and the protons of acyl groups to normalize the results, as described above. I_1-MG in equation 9 stands for the area of one proton from glyceryl group in 1-MG (δ = 3.94 ppm), I_2-MG in equation 10 for the area of one proton from glyceryl group in 2-MG (δ = 4.93 ppm) and I_1,2-DG in equation 11 for the area of one proton from glyceryl group in 1,2-DG (δ = 5.08 ppm).

(9) $C_{1-\mathrm{MG}}=6\times \left(I_{1-\mathrm{MG}}/I_{\mathrm{Acy}}\right)\times 100\%$

(10) $C_{2-\mathrm{MG}}=6\times \left(I_{2-\mathrm{MG}}/I_{\mathrm{Acy}}\right)\times 100\%$

(11) $C_{1,2-\mathrm{DG}}=6\times \left(I_{1,2-\mathrm{DG}}/I_{\mathrm{Acy}}\right)\times 100\%$

For determination of TG and 1,3-DG equations 12 and 13 were adapted from Nieva-Echevarría et al.²¹ Due to the overlapping of TG, and 1,2-DG protons in the integration region of the spectrum to determine C_TG, the I_1,2-DG protons signals were subtracted from the I_TG region. In equation 12, I_TG stands for the region between δ = 4.22 and δ = 4.40 ppm, where half of the signal (2 of 4 protons) from the glyceryl group present in TG is considered and half to the signal (1 of 2 protons) from the glyceryl group present in 1,2-DG. However, for C_1,3-DG determination all glyceryl protons from 1,3-DG must be considered and their overlapping with TG, 1-MG and 1,2-DG. In this way, the integration region (I_1,3-DG) goes from δ = 4.00 to δ = 4.40 ppm, and I_TG, I_1-MG and I_1,2-DG must be subtracted.

(12) $C_{\mathrm{TG}}=3\times \left(\left(I_{\mathrm{TG}}-I_{1,2-\mathrm{DG}}\right)/I_{\mathrm{Acy}}\right)\times 100\%$

(13) $\begin{array}{c}{C}_{1,3-\mathrm{DG}}=1.2\times \left[I_{1,3-\mathrm{DG}}\right.-\left(2\times I_{1-\mathrm{MG}}\right)-\left(2\times I_{1,2-\mathrm{DG}}\right)-\hfill \\ \left.\left(2\times \left(I_{\mathrm{TG}}-I_{1,2-\mathrm{DG}}\right)\right)\right]/I_{\mathrm{Acy}}\times 100\%\hfill \end{array}$

With those five equations all major and minor constituents of the reaction medium could be quantified, except methanol and glycerol whom were washed out previously to ¹H NMR analysis.

The possibility of determining the molar conversion of the different species over time allowed the monitoring of the kinetics of the transesterification reaction. With the molar conversion values of each species at the evaluated times it was possible to determine and optimize the values of the formation and reaction rate constants of the reaction by the super modified simplex method. Since it is a recursive and constrained method, values of upper and lower limits are specified for each parameter (Table 1). The choice of the rate constant ranges was established based on preliminary tests.

Table 1
Upper and lower limit of rate constants used in simplex optimization

Independent variables
Rate constants / (m3 mol-1 s-1)	Lower limit	Upper limit
k1f	0.20 × 10-6	0.33 × 10-6
k2f	0.90 × 10-6	1.20 × 10-6
k3f	2.80 × 10-6	3.3 × 10-6
k1r	0.06 × 10-6	0.26 × 10-6
k2r	0.50 × 10-6	1.20 × 10-6
k3r	0.15 × 10-6	0.39 × 10-6

Simplex optimization begins with the delimitation of lower and upper limits for each factor that will be controlled. Simplex optimization was performed in three fundamental steps: initial simplex; coordinate of reflection and coordinate of the motion vector, according to the methodology described by Bona et al.¹⁰ and Bordin et al.⁵ In the initial simplex, the algorithm suggests combinations of values of rate constants, that were used in the simulations of the reaction in geometry 0D. After establishing the initial simplex, the molar conversion values for each simplex were obtained with the FEM application. Then, new values of rate constants were suggested by the algorithm that was inserted into the reaction modeling interface on the COMSOL Multiphysics^® software.¹⁸ Thus, all values of the constants were automatically oriented towards the optimal response.¹⁰

In the optimization, the values of the formation and reverse rate constants presented stabilization from the simplex number 15 for the six variables evaluated (Figures 1a-b). Although the desirability function (D) has stabilized around the simplex 15, the simplex number 20 presented higher value of D (Figure 1c). This happens because is a joint optimization, that is, a maximization of the objective function, so that many combinations of the variables can cause oscillations in the values of the analyzed variables. The optimization is reached when the variables have stabilized, in this way, it is known that the simplex has reached its optimum and is no longer modified.⁵

Figure 1
Stabilization of forward and reverse rate constants during the optimization process.

Table 2 shows the values of the formation and reverse rate constants of the reactions (equations 2-4) optimized by the simplex method. The values were adjusted based on the experimental data of the molar conversion rates of triglycerides, monoglycerides, diglycerides, methanol, glycerol and FAME according to the proposed reaction mechanism.

Table 2
The kn values of formation and reverse reaction optimized by the super-modified simplex method

Rate constants / (10-6 m3 mol-1 s-1)
k1f	k2f	k3f	k1r	k2r	k3r
0.2498	1.1365	3.1336	0.2021	0.8842	0.2185

According to Table 2, the rate constant for the forward reactions at low temperatures normally has a magnitude higher than the corresponding reverse reaction.^3,17,22 For the rate constants of formation, the values decrease in the order of k_3f > k_2f > k_1f. This behavior is not observed for the reverse reactions. Considering only formation reactions, the reaction from TG to DG (k_1f) is slower than DG for MG (k_2f) and MG for GL (k_3f). This explains why there is accumulation of TG in the initial minutes of the reaction.

Low k_1f values for the first formation reaction were also reported by Vicente et al.¹⁵ and Vicente et al.,¹⁶ meaning that at low temperatures the conversion from triglyceride to diglyceride was the slowest reaction and was the one that controlled the process.

The low values of reverse rate constants indicate that these reactions are not favorable. According to Vicenteet al.,¹⁵ the values for the third reverse reaction (k_3r), where the reaction of the glycerin with the methyl ester to give monoglyceride and methanol was not favored, is considered irreversible. This fact is due to the immiscibility of methyl esters and glycerol, which involved a great mass transfer resistance in that direction.

The kinetic experimental profiles and modeling curves generated using the equations proposed for the determination of the conversion rates of the different compounds are shown in Figure 2. It should be noted that each value is the average of three observations for reactions performed under the same conditions.

Figure 2
Kinetic modeling curves (- simulated) and mean experimental points (● experimental) obtained by HF-1H NMR for the species present in the reaction.

The Figure 2 shows the rate of consumption of TG and formation of FAME and GL as well as the intermediate compounds. The intermediaries MG and DG, increased in the first few minutes of the reaction achieving a maximum, then decreased and finally stayed nearly constant.

According to Figure 2, the initial stages of the reaction, is slow, possibly by the temperature that the reactions were conducted and the immiscibility between TG and methanol and stirring intensity. From the 7 min forward is possible to observe a kinetically controlled region namely fast. Finally, a slow region, approximately 50 min, which represents the equilibrium state of the reaction.⁴

The experimental results showed that the second order kinetic model adequately describes the reaction conditions. In the research of Stamenkovic et al.,³ a similar kinetic profile was obtained using a 6:1 molar ratio of methanol to sunflower oil with 1% m/m KOH at a temperature of 20 °C at 200 rpm when proposed a simple kinetic model, for simulation of the TG conversion and the FAME formation in the latter regime: the fast irreversible second-order reaction was followed by the slow reversible second-order reaction close to the completion of the methanolysis reaction.

The reaction yields with the experimental data was 85% after 62 min reaction time, which had ± 2.9% deviations of simulated results (Table 3). The low deviation between the experimental and simulated data indicates that the FEM can be used to solve the kinetic model.

Table 3
The deviation between experimental data and simulated results

	Experimental	Simulated	Deviation
Conversion / %	85.43	87.94	2.94

Conclusions

The kinetic empirical model proposed for the transesterification reaction can be simulated by the finite element method under realistic conditions. The simulation of the reaction by FEM associated to the simplex optimization for the values of the rate constants were compared with the experimental results, validating its application. The estimated and optimized rate constants presented behavior similar to those reported by Stamenkovic et al.³ under similar conditions. This new approach allows a more rigorous control of the process as well as modulates the reaction time. Presenting encouraging prospects for kinetic studies in all types of environments.

Acknowledgments

The State University of Londrina (UEL) for the technical support and CAPES for granting the scholarship. The authors also thank the Spectroscopy Laboratory (LABSPEC) - UEL.

References

• ¹
  Chendynski, L. T.; Mantovani, A. C. G.; Savada, F. Y.; Messias, G. B.; Santana, V. T.; Salviato, A.; Di Mauro, E.; Borsato, D.; Fuel 2019, 242, 316.
• ²
  Oliveira, F. C.; Coelho, S. T.; Renewable Sustainable Energy Rev 2017, 75, 168.
• ³
  Stamenković, O. S.; Todorović, Z. B.; Lazić, M. L.; Veljković, V. B.; Skala, D. U.; Bioresour. Technol 2008, 99, 1131.
• ⁴
  Noureddini, H.; Zhu, D.; J. Am. Oil Chem. Soc 1997, 74, 1457.
• ⁵
  Bordin, M. S. P.; Borsato, D.; Cremasco, H.; Galvan, D.; Silva, L. R. C.; Romagnoli, É. S.; Angilelli, K. G.; Food Chem 2019, 273, 99.
• ⁶
  Qian, W.; Wu, J.; Yang, L.; Lin, X.; Chen, Y.; Chen, X.; Xiong, J.; Bai, J.; Ying, H.; Chem. Eng. J. 2012, 197, 424.
• ⁷
  Bona, E.; da Silva, R. S. S. F.; Borsato, D.; Silva, L. H. M.; de Souza Fidelis, D. A.; J. Food Eng 2007, 79, 771.
• ⁸
  Borsato, D.; Moreira, M. B.; Moreira, I.; Pina, M. V. R.; Silva, R. S. S. F.; Bona, E.; Food Sci. Technol 2012, 32, 281.
• ⁹
  Spendley, W.; Hext, G. R.; Himsworth, F. R.; Technometrics 1962, 4, 441.
• ¹⁰
  Bona, E.; Borsato, D.; Silva, R.; Herrerra, R.; Acta Sci 2000, 22, 1201.
• ¹¹
  Bezerra, M. A.; dos Santos, Q. O.; Santos, A. G.; Novaes, C. G.; Ferreira, S. L. C.; de Souza, V. S.; Microchem. J. 2016, 124, 45.
• ¹²
  Nelder, J. A.; Mead, R.; Comput. J. 1965, 7, 308.
• ¹³
  Neto, B. B.; Scarminio, I. S.; Bruns, R. E.; Planejamento e Otimização de Experimentos, 1^st ed.; Editora da UNICAMP: Campinas, São Paulo, Brazil, 1996.
• ¹⁴
  MATLAB; The MathWorks, Inc.: Natick, MA, USA, 2016.
• ¹⁵
  Vicente, G.; Martínez, M.; Aracil, J.; Esteban, A.; Ind. Eng. Chem. Res 2005, 44, 5447.
• ¹⁶
  Vicente, G.; Martinez, M.; Aracil, J.; Energy Fuels 2006, 20, 1722.
• ¹⁷
  Bambase, M. E.; Nakamura, N.; Tanaka, J.; Matsumura, M.; J. Chem. Technol. Biotechnol 2007, 82, 273.
• ¹⁸
  COMSOL Multiphysics^®, version 5.2; COMSOL, Inc., Burlington, MA, USA, 2016.
• ¹⁹
  Darnoko, D.; Cheryan, M.; J. Am. Oil Chem. Soc 2000, 77, 1263.
• ²⁰
  Derringer, G.; Suich, R.; J. Qual. Technol 1980, 12, 214.
• ²¹
  Nieva-Echevarría, B.; Goicoechea, E.; Manzanos, M. J.; Guillén, M. D.; Food Res. Int 2014, 66, 379.
• ²²
  Silva, N. L.; Rivera, E. C.; Batistella, C. B.; Maciel Filho, R.; Maciel, M. R. W.; Comput.-Aided Chem. Eng 2008, 25, 1101.

Publication Dates

• Publication in this collection
  20 Jan 2020
• Date of issue
  Feb 2020

History

• Received
  07 May 2019
• Accepted
  31 July 2019