Biodiesel Production by Reactive Flash: A Numerical Simulation

1Centro 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 2Department of Chemical Engineering, Norwegian University of Science and Technology (NTNU), 7034 Trondheim, Norway 3Universidad del Mar, Ciudad Universitaria SN, Puerto Ángel, 70902 San Pedro Pochutla, OAX, Mexico


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][6][7][8][9][10][11][12][13][14][15][16][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 2 International Journal of Chemical Engineering 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 y.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 steadystate 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.

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  CIC ∈ [0, 1] 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 PLUS5 software.An algorithm flowchart to solve the system of DAEs is shown in Figure 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) x CIC Read parameters: F, T 0 , P, z i , a, C p , q ext , k i , 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).
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 ( bbp ), 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 ( sat  ) 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 ( 1 ,  2 , and  3 ) of each reaction are given as follows: where  TG ,  DG ,  MG ,  G ,  AE , and  ROH represent the concentration of triglyceride, diglyceride, monoglyceride, glycerol, alkyl esters (biodiesel), and alcohol, respectively.Also,  1 ,  2 , . .., and  6 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 Δ  = −5.07× 10 3 J/mol for each reaction, calculated from the formation heat data [36].
Definition 1.Let   ∈ [0, 1] 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  bbp 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  bbp .
Even for a simple system of DAEs, it may not be satisfied for some values of  bbp .
The manifold point  is the state space for the dynamic system of DAEs defined by ( 9), which induces a vector field on .
Therefore, a singular point ( ss ,  ss ) ∈ ( ss ) at the parameter  = [      ]  is an equilibrium point.
At  =  bbp , 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:

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,  rxn 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  = −100 K mol −1 .A perfect holdup is assumed ( is constant).Hence, / = 0 implies that   = 1 − .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  = 5.8 minutes is around 0.004.This implies that a 97% conversion was reached.
The two-phase operating mode does not exhibit steadystate 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  rxn ,  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.
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  Figure 5: Behavior of the system trajectories as the manifold is attained.[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  rxn 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.

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.

Figure 3 :
Figure 3: Dynamic behavior of the reactive flash process.

Figure 4 :
Figure 4: Bifurcation diagram for  = −100 K mol −1 as a function of the parameter  rxn .