International Journal of Chemical Engineering
Volume 2016 (2016), Article ID 7843081, 8 pages
http://dx.doi.org/10.1155/2016/7843081

Research Article

Biodiesel Production by Reactive Flash: A Numerical Simulation

Alejandro Regalado-Méndez,1,2,3 Sigurd Skogestad,² Reyna Natividad,¹ and Rubí Romero¹

¹Centro Conjunto de Investigación en Química Sustentable UAEMex-UNAM, Universidad Autónoma del Estado de México, Carretera Toluca-Atlacomulco Km 14.5, Unidad San Cayetano, 50200 Toluca, MEX, Mexico
²Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU), 7034 Trondheim, Norway
³Universidad del Mar, Ciudad Universitaria SN, Puerto Ángel, 70902 San Pedro Pochutla, OAX, Mexico

Received 26 December 2015; Revised 11 March 2016; Accepted 4 April 2016

Academic Editor: Deepak Kunzru

Copyright © 2016 Alejandro Regalado-Méndez et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Reactive flash (RF) in biodiesel production has been studied in order to investigate steady-state multiplicities, singularities, and effect of biodiesel quality when the RF system approaches to bubble point. The RF was modeled by an index-2 system of differential algebraic equations, the vapor split () was computed by modified Rachford-Rice equation and modified Raoult’s law computed bubble point, and the continuation analysis was tracked on MATCONT. Results of this study show the existence of turning points, leading to a unique bubble point manifold, , which is a globally stable flashing operation. Also, the results of the simulation in MATLAB® of the dynamic behavior of the RF show that conversion of triglycerides reaches 97% for a residence time of 5.8 minutes and a methanol to triglyceride molar flow ratio of 5 : 1.

1. Introduction

In the last decade, biofuels production has been worldwide motivated by the need of reducing greenhouse gases in order to slow down climate change. In addition, biofuels production has also been encouraged by the unstable oil prices, the reduction of petroleum reserves, and environmental penalties in stationary and mobile sources of pollution. In this context, biodiesel emerges as a viable renewable energy alternative to petroleum diesel. This clean renewable fuel is superior to diesel oil in terms of sulfur and aromatic content. Also, it is environmentally safe, nontoxic, and biodegradable. In general, biodiesel is derived from a transesterification reaction of triglycerides in vegetable oils or animal fats with alcohol (i.e., methanol and ethanol) under the presence of catalyst [1]. Biodiesel manufacturing at large scale demands the development of innovative and efficient processes. In this sense, reactive distillation (RD) seems to offer interesting and desirable features. RD is a hybrid process where a chemical transformation and a separation take place in the same vessel. The product is removed at the same time that it is formed. This characteristic makes it possible to overcome the thermodynamic equilibrium of the reaction as well as the vapor-liquid equilibrium (VLE) [2]. RD have many advantages, which include reduced capital cost, complete conversion in equilibrium-limited reactions, heat integration, and reduced waste generation [3, 4]. However, this makes of RD a rather complex hybrid process because of the interactions between reaction and phase equilibrium. This process is highly nonlinear and, at time, postures challenges in terms of operation and control. Multiple steady states exist in conventional chemical reactors operated in continuous mode (i.e., continuous stirred tank reactors) and in nonreactive and reactive distillation columns. Investigations of the physical phenomena that cause multiplicities of steady states and Hopf bifurcations in RD have been reported by many researchers around the word [5–17]. These studies show that there are various reasons for the occurrence of output and input multiplicity, which include multiple steady states caused by heat effects, or kinetic instabilities [6, 18, 19], or presence of singularities due to nonlinear transformation and decoupling of energy balance [20]. Existence of azeotropes can also lead to multiplicity in both nonreactive and reactive distillation columns [21] as well as the presence of reaction hysteresis due to the interaction between nonreactive and reactive sections [12]. Moreover, according to Waschler et al. [22], existence of multiple steady states has been identified for light-boiling reactant, sufficiently large difference in boiling points, and reaction orders lower than the physical parameter α (the α parameter is a measure of the phase-equilibrium-driven self-inhibition of the reaction mechanism), and Purohit et al. [17] propose a new relatively simple method to identify multiplicity in reactive distillation due to the interaction of reaction and distillation. Their method compared to the bifurcation technique makes the analysis simpler and provides more insight into the influences of different factors on multiplicity behavior. Despite the many investigations of these important and interesting phenomena in RD, only a few efforts to understand the basic causes of steady states multiplicity have been reported in the literature. Few investigators studied the multiplicity of steady states behavior of the reactive flash in order to understand these phenomena in RD systems. Rodríguez et al. [23] demonstrated that input and output multiplicity lead to an isobaric, adiabatic reactive flash for a binary mixture due to the presence of vapor-liquid equilibrium. Lakerveld et al. [24] studied the steady-state behavior of a reactive flash by singularity theory. The reactive system was characterized by an exothermic isomerization reaction with first-order kinetics and a light-boiling reactant. By using the Damköhler number as continuation parameter, twenty-five bifurcation diagrams were found, exhibiting three steady states and five feasibility boundaries, when the heat of reaction, the activation energy, or the relative volatility is increased. Rodríguez et al. [7] investigated an isobaric reactive flash with controlled kinetics and a constant split fraction. They found that in this situation the coupling of the energy balance does not drive multiplicities with the material balances and phase equilibrium equations. Their solutions revealed multiplicity when the heat of reaction was small compared to the heat of vaporization. This provided some insight into causes of multiplicity in RD columns when the constant molal overflow approximation works well. Furthermore, Ruiz et al. [25] and Ruiz et al. [8] demonstrate the existence of Hopf bifurcation, limit points, and isolas with intersecting branches for an equilibrium and nonequilibrium reactive separation process when this is modeled by a system of ordinary differential equations and a set of differential algebraic equations.

Jaime-Leal et al. [9] introduced a new approach and conditions to identify input multiplicity in reactive flash, based on the application of reaction-invariant composition variables. Finally, Harney et al. [26] demonstrate that dynamic behavior of reactive flash and reactive distillation represented by index-2 system of differential algebraic equations (DAEs) can be reduced to an ordinary differential equations system (ODE) by single differentiation; the resulting Jacobian matrix will have a null eigenvalue at every steady state with multiplicity of at least the dimension . These null eigenvalues must be accounted for when determining the stability of a steady state. Moreover, Alvarez-Ramirez [27] found singular dynamics in a simple reactive flash model and described different steady-state operating scenarios with at least one of them belonging to the one-phase operation of the system, including a globally stable flashing to unfeasible operation leading to emptiness of liquid phase. As noted earlier, steady-states multiplicity can occur in RD columns due to high nonlinearity [4, 22]. These multiplicities have impact on desired high conversion at steady state, quality of products, dynamic behavior, and performance and control system design for RD systems. Recent studies on control in RD systems have referred to the need to avoid controlling outputs exhibiting multiplicity [28, 29]. A short review of studies of steady-state multiplicity focused on reactive flash and distillation shows that few investigators found steady-state multiplicity when the temperature of separation process or temperature of reactive system leads to boiling point, and the split fraction is constant. In the present work, authors search for dynamic singularity by dimensionless heat of reaction rate when the system temperature achieves the bubble point temperature. This issue is accomplished by studying an isobaric, reactive flash without a constant split fraction and a biodiesel production from triglycerides with alkali catalyst. Specifically, reactive flash model is developed to analyze the operating conditions and is numerically solved to illustrate possible steady-state multiplicity. The dynamic behavior is represented by an index-2 system of differential algebraic equations.

2. Simulation Methodology

The dynamic simulations were numerically carried out in MATLAB by solving the model represented by (1)–(4) using Petzold’s method [30], where the consistent initial conditions [31] are given for and bubble point temperature. In the steady state, the system of DAEs forms a set of nonlinear algebraic equations. To solve them in MATLAB, the Newton-Raphson method was applied. Besides, the identification of multiple steady states was tracked using MATCONT from MATLAB toolbox [32]. Also, the vapor split fraction () and the bubble point temperature were computed using the modified Rachford-Rice equation [25] and modified Raoult’s law by Newton-Raphson method, respectively. For simulation purposes, the fluid package was set as Wilson equation to compute the activity coefficients as recommended by Carlson [33] and Suthar and Joshipura [34]. Also, extended Antoine equation to compute the vapor pressure was used. Finally, the Antoine constants and binary coefficients for Wilson equation were obtained from ASPEN Properties PLUS® software. An algorithm flowchart to solve the system of DAEs is shown in Figure 1.

Figure 1: Algorithm flowchart reactive flash solutions.

2.1. Model

The reason for modeling the reactive flash as a DAE system is that the equations for the vapor-liquid equilibrium calculations are implicit. Hence, the reactive flash model is obtained from mass and heat balances in dynamic conditions. Figure 2 illustrates the reactive flash process. The assumptions for the employed model are as follows: (i) the chemical reaction is accomplished in the homogeneous liquid phase, (ii) vapor and liquid phases are well mixed, and (iii) the vapor holdup is negligible compared to the liquid holdup (). The reactive flash model is represented by the set of differential algebraic equations given in (1)–(4).

Figure 2: Diagram of reactive flash adapted from [23, 24].

The dynamic representation model of a reactive flash with reactions, components, nonconstant molar holdup , liquid molar feed rate , and nonideal vapor-liquid equilibrium behavior can be written as follows.

As the liquid mixture achieves its bubble point (), evaporation rises spontaneously. The total mass balance is given by

The component mass balance is described as follows:

By considering heating properties, the energy balance is given by

The restrictions of the system are represented by the algebraic equations (4), which are the thermodynamic considerations (i.e., equilibrium equation and the -value prediction):

The saturation pressure () is given by the extended Antoine’s equation and activity coefficients () are computed by Wilson model as recommended by Carlson [33] and Suthar and Joshipura [34]. The vapor split fraction () was computed by using the modified Rachford-Rice equation [25]:

2.2. Kinetic Model

The particular case considered here is an adaptation of the transesterification of triglycerides (oils/fats) by reaction with alcohol in the presence of NaOH as catalyst to produce fatty acid alkyl esters and glycerol. The reaction proceeds in three steps as shown in the following reactions [36]:

The reaction rates (, , and ) of each reaction are given as follows:

where , , , , , and represent the concentration of triglyceride, diglyceride, monoglyceride, glycerol, alkyl esters (biodiesel), and alcohol, respectively. Also, , and are forward and backward rate constants. The temperature dependency of the rate constant () is expressed by Arrhenius’ law [37]:

The kinetic parameters are shown in Table 1. The heat of the reaction is  J/mol for each reaction, calculated from the formation heat data [36].

Table 1: Kinetic parameters [35].

Component
(L/mol min)
(J/mol)

Triglyceride
1.469 × 10⁸
58782.6720

Alcohol
105100
44962.0452

Diglyceride
1.19 × 10¹⁰
67193.9532

Alkyl ester
1.725 × 10⁸
58225.8276

Monoglyceride
24940
30031.9164

Glycerol
1.469 × 10⁸
46042.2396

2.3. Global Regularity Condition on Manifold Point

Consider an index-2 system of differential algebraic equations given by

where (), , and and are both in an open neighborhood, , of in .

Definition 1. Let and be the liquid mol fraction of component th and the system pressure, respectively. Given a thermodynamic model based on activity coefficients, the bubble point temperature is given, in general as implicit equation of the form:

is a regular manifold point of dimension , if

The structure of depends, of course, on the parameter . Even for a simple system of DAEs, it may not be satisfied for some values of .

The manifold point is the state space for the dynamic system of DAEs defined by (9), which induces a vector field on .

The vector field may not be well defined at all points of . At any point , we have and is not singular; the is uniquely defined by

If = is singular at point , then the vector field is not well defined at that point. Typically, such singular point lies in codimension 1 submanifolds of [24]. The next theorem is written in adapted form from previous works such as Rheinboldt [38], Sudarsan and Keerthi [39], and Okay et al. [40].

Theorem 2. Suppose is regular manifolds point for all near , and that at point , . Then, is said to be causal.

Lemma 3. A singular point of the systems of the DAEs in steady state at given parameter is also an equilibrium point.

Proof. Suppose is a singular point of the decoupled DAEs at the parameter value such that

Therefore, a singular point at the parameter is an equilibrium point.

At , the evaporation at the bubble point spontaneously occurs when the excessive heat in the liquid phase is captured by transferring it into vapor phase. It is worth pointing out that (14) is rather important in the modeling of reactive flash, since the reactive separation process can operate in one or two phases. In consequence, the point will be called bubble point manifold:

3. Results

The set of DAEs system is conformed by 8 differential equations and 12 algebraic constraints, which come from mass and energy balances and thermodynamic considerations (equilibrium equation and -value prediction), respectively. The main variables are residence time , liquid and vapor mole fractions, and temperature. Constants are activity coefficients that come from Wilson model equation. Also, vapor split fraction can be computed from modified Rachford-Rice equation. Finally, is selected as bifurcation parameter, since it can induce large steady-state changes in composition and temperature, as shown in Figure 4.

The index-2 system of DAEs was solved by using ode15i MATLAB toolbox when  K mol⁻¹. A perfect holdup is assumed ( is constant). Hence, implies that . Also, the molar flow ratio 1 : 5 of triglyceride/alcohol was chosen to perform the dynamic behavior and the phase maps because this relationship leads to the maximum of biodiesel molar fraction [36]. Figure 3 depicts the dynamic behavior of the molar fraction for the most important components. Also, the settling time, , of the reactive flash and the steady-state molar fraction for biodiesel were determined to be 5.8 minutes and 0.46, respectively. In addition, the molar fraction of triglyceride at minutes is around 0.004. This implies that a 97% conversion was reached.

Figure 3: Dynamic behavior of the reactive flash process.

Figure 4: Bifurcation diagram for as a function of the parameter .

The two-phase operating mode does not exhibit steady-state multiplicity. The overall steady-state multiplicity is introduced by the multiplicity of the one-phase operating mode. Figure 4 shows the continuation path for the biodiesel molar fraction and steady-state stability properties evaluated at large values of the dimensionless reaction enthalpy , where a region of steady-state multiplicity exists. The dashed line indicates the transition from unstable to stable steady states, which is delineated by turning points, respectively. In addition, note that the two-phase branch emerges in a nondifferential way from the one-phase branch. In other words, the bifurcation diagram is not differentiable at = 0. This fact is a consequence of the discontinuous nature of the reactive flash process because the system operation goes from one- to two-phase mode. Also, the blue line leads to low molar fraction biodiesel; meanwhile, the red dashed line leads to high molar fraction biodiesel corresponding to liquid phase operation and black line reaches a unique equilibrium point corresponding to a globally stable flashing operation.

The singular trajectories of the system of DAEs are reflected in approaches to the bubble point manifold. Figure 5 illustrates the phenomena by showing the behavior of the system trajectories as the bubble point manifold is approached. It can be observed that the dynamics trajectories converge to a stable equilibrium (0.46, 478.41 K). The dynamics trajectories display a sliding behavior as they achieve the bubble point manifold to converge to a stable equilibrium. The singular nature of these dynamics as the bubble point manifold is attained by interaction between liquid-vapor separation and chemical transformation.

Figure 5: Behavior of the system trajectories as the manifold is attained.

The reactive flash drawing has an asymptotic dynamic behavior such as that exhibited by continuous stirred tank chemical reactors.

The topological form of the bifurcation shapes of this work is of type “” such as typically displayed in continuous stirred tank chemical reactor and reactive flash. This is supported by other studies (Rodríguez et al. [7], Ruiz et al. [8], Alvarez-Ramirez [27], Jaime-Leal et al. [9], and Harney et al. [26]).

The shapes of the trajectories of the bifurcation map displayed in Figure 5 are similar to spiral points (stable node) [41, 42] as is sketched in autonomous differential equation system, for example, in stirred tank chemical reactors.

Finally, modeling and analyzing a reactive flash provide important insights for understanding the design, operation, and control of higher order process. For example, in this case study, the feasible flashing region is not too sensitive to dimensionless reaction enthalpy at large values, because it implies relatively small changes of composition and temperature. Therefore, implementing composition control and temperature control in reactive stages should be avoided on reactive distillation.

4. Conclusions

The singularities in reactive flash when the system temperature approaches to the bubble point were investigated. It is shown that the differential algebraic equations display steady-state multiplicities, revealing the existence of turning points and leading to the bubble point manifold, this unique equilibrium point corresponding to a globally stable flashing operation. Results indicate that reactive distillation stages can display discontinuous dynamics by the large changes in the heat of reaction rate.

Nomenclature

: | Liquid feed rate ()
: | Liquid molar holdup (mol)
: | Liquid exit flow ()
: | Specific liquid molar exit rate (=)
: | Vapor exit flow ()
: | Temperature (K)
: | Liquid mole fraction of component
: | Feed mole fraction of component
: | Vapor mole fraction of component
: | Total pressure
: | Vapor pressure of component (Pa)
: | ⁡ -value prediction
: | Rate of reaction for component in reaction ()
: | Total rate of reaction
: | Heat capacity ()
: | External heat input to flash, divided by feed rate ()
: | Total heat input to flash (=) ().

Greek

: | Residence time (=)
: | Vapor split fraction =
: | Stoichiometric coefficient of component in reaction
: | Reaction volume
: | Enthalpy of reaction (J)
: | Enthalpy of vaporization (J)
: | Modified enthalpy of reaction (=) (K)
: | Modified enthalpy of vaporization (=) (K).

Competing Interests

The authors declare that they have no competing interests regarding the publication of this paper.

Acknowledgments

The authors are grateful to PRODEP (before PROMEP) for financial support through PROMEP/103.5/10/5552.

References

A. C. Dimian, C. S. Bildea, F. Omota, and A. A. Kiss, “Innovative process for fatty acid esters by dual reactive distillation,” Computers & Chemical Engineering, vol. 33, no. 3, pp. 743–750, 2009.

R. A. Pai, M. F. Doherty, and M. F. Malone, “Design of reactive extraction systems for bioproduct recovery,” AIChE Journal, vol. 48, no. 3, pp. 514–526, 2002.

K. Alejski and F. Duprat, “Dynamic simulation of the multicomponent reactive distillation,” Chemical Engineering Science, vol. 51, no. 18, pp. 4237–4252, 1996.

R. Taylor and R. Krishna, “Modelling reactive distillation,” Chemical Engineering Science, vol. 55, no. 22, pp. 5183–5229, 2000.

R. Baur, R. Taylor, and R. Krishna, “Bifurcation analysis for TAME synthesis in a reactive distillation column: comparison of pseudo-homogeneous and heterogeneous reaction kinetics models,” Chemical Engineering and Processing: Process Intensification, vol. 42, no. 3, pp. 211–221, 2003.

F. Chen, R. S. Huss, M. F. Doherty, and M. F. Malone, “Multiple steady states in reactive distillation: kinetic effects,” Computers & Chemical Engineering, vol. 26, no. 1, pp. 81–93, 2002.

I. E. Rodríguez, A. Zheng, and M. F. Malone, “Parametric dependence of solution multiplicity in reactive flashes,” Chemical Engineering Science, vol. 59, no. 7, pp. 1589–1600, 2004.

G. Ruiz, M. Diaz, and L. N. Sridhar, “Singularities in reactive separation processes,” Industrial & Engineering Chemistry Research, vol. 47, no. 8, pp. 2808–2816, 2008.

J. E. Jaime-Leal, A. Bonilla-Petriciolet, J. G. Segovia-Hernández, S. Hernández, and H. Hernández-Escoto, “Analysis and prediction of input multiplicity for the reactive flash separation using reaction-invariant composition variables,” Chemical Engineering Research and Design, vol. 90, no. 11, pp. 1856–1870, 2012.

R. Monroy-Loperena and J. Alvarez-Ramirez, “On the steady-state multiplicities for an ethylene glycol reactive distillation column,” Industrial & Engineering Chemistry Research, vol. 38, no. 2, pp. 451–455, 1999.

A. Kumar and P. Daoutidis, “Modeling, analysis and control of ethylene glycol reactive distillation column,” AIChE Journal, vol. 45, no. 1, pp. 51–68, 1999.

M. G. Sneesby, M. O. Tadé, and T. N. Smith, “Reaction hysteresis: a new cause of output multiplicity in reactive distillation,” Developments in Chemical Engineering and Mineral Processing, vol. 7, no. 1-2, pp. 41–56, 1999.

M. A. Al-Arfaj and W. L. Luyben, “Comparative control study of ideal and methyl acetate reactive distillation,” Chemical Engineering Science, vol. 57, no. 24, pp. 5039–5050, 2002.

A. P. Higler, R. Taylor, and R. Krishna, “Nonequilibrium modelling of reactive distillation: multiple steady states in MTBE synthesis,” Chemical Engineering Science, vol. 54, no. 10, pp. 1389–1395, 1999.

T. E. Güttinger and M. Morari, “Predicting multiple steady states in equilibrium reactive distillation. 2. Analysis of hybrid systems,” Industrial & Engineering Chemistry Research, vol. 38, no. 4, pp. 1649–1665, 1999.

T. E. Güttinger and M. Morari, “Predicting multiple steady states in equilibrium reactive distillation. 1. Analysis of nonhybrid systems,” Industrial & Engineering Chemistry Research, vol. 38, no. 4, pp. 1633–1648, 1999.

J. L. Purohit, S. M. Mahajani, and S. C. Patwardhan, “Analysis of steady-state multiplicity in reactive distillation columns,” Industrial & Engineering Chemistry Research, vol. 52, no. 14, pp. 5191–5206, 2013.

A. Uppal, W. H. Ray, and A. B. Poore, “On the dynamic behavior of continuous stirred tank reactors,” Chemical Engineering Science, vol. 29, no. 4, pp. 967–985, 1974.

L. F. Razón and R. A. Schmitz, “Multiplicities and instabilities in chemically reacting systems—a review,” Chemical Engineering Science, vol. 42, no. 5, pp. 1005–1047, 1987.

E. W. Jacobsen and S. Skogestad, “Multiple steady states in ideal two-product distillation,” AIChE Journal, vol. 37, no. 4, pp. 499–511, 1991.

N. Bekiaris, G. A. Meski, C. M. Radu, and M. Morari, “Multiple steady states in homogeneous azeotropic distillation,” Industrial and Engineering Chemistry Research, vol. 32, no. 9, pp. 2023–2038, 1993.

R. Waschler, S. Pushpavanam, and A. Kienle, “Multiple steady states in two-phase reactors under boiling conditions,” Chemical Engineering Science, vol. 58, no. 11, pp. 2203–2214, 2003.

I. E. Rodríguez, A. Zheng, and M. F. Malone, “The stability of a reactive flash,” Chemical Engineering Science, vol. 56, no. 16, pp. 4737–4745, 2001.

R. Lakerveld, C. S. Bildea, and C. P. Almeida-Rivera, “Exothermic isomerization reaction in a reactive flash: steady-state behavior,” Industrial and Engineering Chemistry Research, vol. 44, no. 10, pp. 3815–3822, 2005.

G. Ruiz, L. N. Sridhar, and R. Rengaswamy, “Isothermal isobaric reactive flash problem,” Industrial & Engineering Chemistry Research, vol. 45, no. 19, pp. 6548–6554, 2006.

D. A. Harney, T. K. Mills, and N. L. Book, “Numerical evaluation of the stability of stationary points of index-2 differential-algebraic equations: applications to reactive flash and reactive distillation systems,” Computers and Chemical Engineering, vol. 49, pp. 61–69, 2013.

J. Alvarez-Ramirez, “Singular reactive flash dynamics,” Chemical Engineering and Processing: Process Intensification, vol. 69, pp. 119–125, 2013.

S.-B. Hung, M.-J. Lee, Y.-T. Tang et al., “Control of different reactive distillation configurations,” AIChE Journal, vol. 52, no. 4, pp. 1423–1440, 2006.

M. V. Pavan Kumar and N. Kaistha, “Role of multiplicity in reactive distillation control system design,” Journal of Process Control, vol. 18, no. 7-8, pp. 692–706, 2008.

K. E. Brenan and L. R. Petzold, “The numerical solution of higher index differential/algebraic equations by implicit methods,” SIAM Journal on Numerical Analysis, vol. 26, no. 4, pp. 976–996, 1989.

S. L. Campbell, “Consistent initial conditions for singular nonlinear systems,” Circuits, Systems, and Signal Processing, vol. 2, no. 1, pp. 45–55, 1983.

A. Dhooge, W. Govaerts, Y. A. Kuznetsov, H. G. Meijer, and B. Sautois, “New features of the software MatCont for bifurcation analysis of dynamical systems,” Mathematical and Computer Modelling of Dynamical Systems: Methods, Tools and Applications in Engineering and Related Sciences, vol. 14, no. 2, pp. 147–175, 2008.

E. C. Carlson, “Don't gamble with physical properties for simulations,” Chemical Engineering Progress, vol. 92, no. 10, pp. 35–46, 1996.

K. Suthar and M. Joshipura, “A comparative study on predictions of vapor liquid equilibrium of biodiesel systems,” in Proceedings of the 2nd International Conference on Current Trends in Technology (NUiCONE '11), pp. 1–6, Institute of Technology, Nirma University, Ahmedabad, India, December 2011.

H. Noureddini and D. Zhu, “Kinetics of transesterification of soybean oil,” Journal of the American Oil Chemists' Society, vol. 74, no. 11, pp. 1457–1463, 1997.

M. Agarwal, K. Singh, and S. P. Chaurasia, “Simulation and sensitivity analysis for biodiesel production in a reactive distillation column,” Polish Journal of Chemical Technology, vol. 14, no. 3, pp. 59–65, 2012.

N. de Lima da Silva, E. Ccopa Rivera, C. B. Batistella, D. Ribeiro de Lima, R. Maciel Filhob, and M. R. Wolf Maciel, “Biodiesel production from vegetable oils: operational strategies for large scale systems,” in 18th European Symposium on Computer Aided Process Engineering—ESCAPE 18, B. Braunschweig and X. Joulia, Eds., pp. 1001–1007, Elsevier, Lyon, France, 2008.

W. C. Rheinboldt, “Solving algebraically explicit DAEs with the MANPAK-manifold-algorithms,” Computers and Mathematics with Applications, vol. 33, no. 3, pp. 31–43, 1997.

R. Sudarsan and S. S. Keerthi, “Numerical approaches for solution of differential equations on manifolds,” Applied Mathematics and Computation, vol. 92, no. 2-3, pp. 153–193, 1998.

I. Okay, S. L. Campbell, and P. Kunkel, “Completions of implicitly defined linear time varying vector fields,” Linear Algebra and its Applications, vol. 431, no. 9, pp. 1422–1438, 2009.

S. E. LeBlanc and D. R. Coughanowr, Process Systems Analysis and Control, McGraw-Hill, New York, NY, USA, 3rd edition, 2009.

Z. Bilicki, C. Dafermos, J. Kestin, G. Majda, and D. L. Zeng, “Trajectories and singular points in steady-state models of two-phase flows,” International Journal of Multiphase Flow, vol. 13, no. 4, pp. 511–533, 1987.