# Numerical Analysis for the Synthesis of Biodiesel Using Spectral Relaxation Method

Authors: Z. G. Makukula; S. S. Motsa; S. Shateyi
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Biodiesel is an alternative diesel fuel chemically defined as the
mono-alkyl esters of long chain fatty acids derived from vegetable
oils or animal fat. It is becoming more attractive as an alternative
fuel due to the depleting fossil fuel resources. A mathematical
model for the synthesis of biodiesel from vegetable oils and animal
fats is presented in this study. Numerical solutions of the model
are found using a spectral relaxation method. The method, originally
developed for boundary value problems, is an iterative scheme based
on the Chebyshev spectral collocation method developed by decoupling
systems of equations using Gauss-Seidel type of techniques. The
effects of the reaction rate constants and initial concentrations of
the reactants on the amount of the final product are being
investigated. The accuracy of the numerical results is validated by
comparison with known analytical results and numerical results
obtained usingode45, an efficient explicit 4th and
5th order Runge-Kutta method used to integrate both linear and
nonlinear differential equations.


---

## Body

## 1. Introduction

Due to the continuous uncertainty and increasing environmental impact associated with the utilization of petroleum-based diesel fuel, the demands for alternative fuels have increased considerably in recent years [1]. Methyl and ethyl esters derived from vegetable oil or animal fat, known as biodiesel, have good potential as alternative diesel fuel [2]. Biodiesel is synthesized through a chemical process called transesterification. During this chemical process, triglycerides (TG) in vegetable oils or animal fats react with short chain alcohols such as methanol (MeOH) and ethanol in the presence of homogeneous basic catalysts such as sodium hydroxide (NaOH) and potassium hydroxide (KOH) to form fatty acid methyl esters (biodiesel) and glycerin [3, 4].Biodiesel has many advantages over petroleum-based diesel fuels. It is biodegradable and  nontoxic and produce less particles, smoke, and carbon monoxide [2]. It is also renewable and usable in a variety of diesel engines with minimum or no modification necessary [3, 5]. It burns clean, and do not form engine deposits or generate harmful pollutants which results in a significant reduction of the types of pollutants that contribute to pollution and global warming and releases up to 85% fewer cancer-causing agents [6–10].Experiments have been carried out to investigate different aspects in the production of biodiesel [3, 4, 11–15]. In this work, we carry out a numerical study of the transesterification process for the synthesis of biodiesel from vegetable oils. We develop a mathematical model from the experimental works of Santos et al. [2], Burnham et al. [6], and Gunvachai et al. [7], where the reaction equations and reaction rates are given. In this study, the chemical reaction equations and rates are reduced to a set of nonlinear first order differential equations that are solved using a spectral relaxation method. The method decouples nonlinear systems of equations using ideas similar to those of the Gauss-Seidel iterative scheme to give rise to a linear system which is solved sequentially using the Chebyshev spectral method. The method's applicability to initial value problems is also explored in this work. The effects of the reaction rate constants and initial concentrations of the reactants on the amount of the final product are being investigated. The accuracy of the numerical results is validated by comparison with known analytical results and numerical results obtained using ode45, an efficient explicit 4th and 5th order Runge-Kutta method used to integrate both linear and nonlinear differential equations.


## 2. Mathematical Formulation

Biodiesel (BD) is produced by the transesterification of large, branched triglyceride (TG) molecules into smaller, straight chain molecules with methanol (MeOH) in the presence of an alkali or acid as a (nonreactive) catalyst. The process occurs in three consecutive reactions which in turn produce three by-products, di-glyceride (DG), mono-glyceride (MG), and glycerol (GL) [2, 4, 6, 7, 15, 16]. The thermal transesterification reactions proposed in [12, 13] may be characterized as
(1)TG+MeOH⟶k1DG+R1COOH3,(2)DG+MeOH⟶k2MG+R2COOH3,(3)MG+MeOH⟶k3GL+R3COOH3,
where k1, k2, and k3 are reaction rate constants. The alternative mathematical model proposed by Kusdiana and Saka [14] ignores the two intermediate reactions and assumes the reaction
(4)TG+MeOH⟶kGL+RCOOH3.In this study we consider the three-stage model defined by (1)–(3). Following [4, 6], we assume that there is only one biodiesel product produced in all three reactions so that R1=R2=R3=R. In [6], an adaptive Runge-Kutta integrator was used to simulate the set of ODEs that represent the reaction network defined by (1)–(3). In this study, we derive a system of differential equations that model the change in the concentrations of the reactants and the product species with time from (1)–(3) by using the notation,
(5)x1=[TG],x2=[MeOH],x3=[DG],x4=[BD],x5=[MG],x6=[GL],
where the square bracket denotes the concentration of the given chemical at time t. These equations are
(6)x1′(t)=-k1x1x2,x2′(t)=-k2x1x2-k2x2x3-k3x2x5,x3′(t)=k1x1x2-k2x2x3,x4′(t)=k1x1x2+k2x2x3+k3x2x5,x5′(t)=k2x2x3-k3x2x5,x6′(t)=k3x2x5,
where the prime denotes the time derivative, subject to the initial conditions
(7)x1(0)=x1,0,x2(0)=x2,0,x3(0)=x4(0)=x5(0)=x6(0)=0.In (6), the positive terms in the differential equations are those that contribute to an increase in the xith variable (i=1,…6), that is to say those terms that contribute positively to the concentration of the different products. The negative terms are those contributing to the decay in the concentration of the different species. Effectively, we start with definite known concentrations of the reacting species with no by-products. As t→∞, the concentrations of the reacting species are expected to approach zero and those of the by-products to approach some limiting values, xi,∞  (i=3,4,5,6).


## 3. Method of Solution

In this section, (6) are solved using the spectral relaxation method (SRM). Principally, the method has been developed for similarity boundary layer problems with at least one the essential profiles such as velocity, temperature, or concentration decaying exponentially. The governing systems of equations are decoupled in a manner parallel to the Gauss-Seidel relaxation method normally used for solving linear algebraic system of equations. Worth noting also is that the method bears some similarities with the Jacobi and Gauss-Seidel waveform relaxation methods earlier developed by Lelarasmee et al. [17] for solutions of initial value problems (IVPs) in very large-scale electrical networks. However, these methods may result into nonlinear differential equations which require further linearization using appropriate linearization schemes while the SRM gives rise into a sequence of linear differential equations which can be directly solved using standard methods. Spectral methods [18, 19] are then used in this work to solve the resulting scheme hence, the name spectral relaxation method. Spectral methods have been chosen in this work because of their well-documented advantages including high accuracy. In developing the sequence of linear equations using the SRM, the following guidelines are used.(1)
The equations have to be arranged in such a way that the equation with the least number of unknowns is placed first on the list.(2)
Variables sayZ1,Z2,Z3,…,Zn can be used to assign the ordered equations in [1], such that Zi, (i=1,2,3…) is an unknown function in the ith equation. It is chosen to be the unknown function associated with the highest order derivative of the ith equation.(3)
To generate the iterative scheme, in each equation, only terms linear inZi are to be evaluated at the current iteration level, r+1. All other terms, linear and nonlinear, in the other functions are to be evaluated from the previous iteration, r, with an exception from the second equation, where updated solutions for Zi-1 obtained from the previous i-1 equation will be used. Nonlinear terms in Zi are also evaluated from the previous iteration. Terms made up of a product of Zi and its derivatives, the derivative terms are to be evaluated from r.In view of the guidelines above, we develop the iterative scheme for (6) as follows:
(8)x˙1,r+1+k1x2,rx1,r+1=0,x1,r+1(0)=x1,0,x˙3,r+1+k2x2,rx3,r+1=k1x1,r+1x2,r,x3,r+1(0)=0,x˙2,r+1+(k2x1,r+1+k2x3,r+1+k3x5,r+1)x2,r+1=0,x2,r+1(0)=x2,0,x˙5,r+1+k3x2,r+1x5,r+1=k2x2,r+1x3,r+1,x5,r+1(0)=0,x˙4,r+1=(k1x1,r+1+k2x3,r+1+k3x5,r+1)x2,r+1,x5,r+1(0)=0,x˙6,r+1=k3x2,r+1x5,r+1,x6,r+1(0)=0.To solve the iteration schemes (8) we use the Chebyshev spectral collocation method. For brevity, we omit the details of the spectral methods and refer interested readers to [18, 19]. Before applying the spectral method, it is convenient to transform the domain on which the governing equation is defined to the interval [-1,1] on which the spectral method can be implemented. We use the transformation t=tF(τ+1)/2 to map the interval [0,tF] to [-1,1], where tF is a finite time. The basic idea behind the spectral collocation method is the introduction of a differentiation matrix D which is used to approximate the derivatives of the unknown variables x at the collocation points as the matrix vector product
(9)dxdt=∑k=0N‍Djkx(τk)=DX,j=0,1,…,N,
where N+1 is the number of grid points known as collocation points, D=2D/tF, and X=[x(τ0),x(τ1),…,x(τN)]T is the vector function at the collocation points τj. Applying the Chebyshev spectral collocation method in (8) gives
(10)A1x1,r+1=B1,x1,r+1(ξN)=x1,0,A3x3,r+1=B3,x3,r+1(ξN)=0,A2x2,r+1=B2,x2,r+1(ξN)=x2,0,A5x5,r+1=B5,x5,r+1(ξN)=0,A4x4,r+1=B4,x4,r+1(ξN)=0,A6x6,r+1=B6,x6,r+1(ξN)=0,
where
(11)A1=D+diag⁡(k2x2,r),B1=0,A3=D+diag⁡(k1x2,r),B3=k1x1,r+1x2,r,A2=D+diag⁡(k2x1,r+1+k2x3,r+1+k3x5,r),B2=0,A5=D+diag⁡(k3x2,r+1),B5=k2x2,r+1x3,r+1,A4=D,B4=(k1x1,r+1+k2x3,r+1+k3x5,r+1)x2,r+1,A6=D,B6=k3x2,r+1x5,r+1.From the above equations,diag⁡ is a diagonal matrix and xj,  j=1,…6, are the values of the functions xj when evaluated at the collocation points. To start the SRM iterative scheme (20), we use the initial conditions (7) as initial solutions. The scheme is repeated until the solution converges to the numerical solution obtained using the ode45.


## 4. Results and Discussion

In this section we present SRM results of the governing equations (6) subject to the initial conditions (7). In experiments, the contents of TG, DG, MG, and GL in samples of reaction product are analyzed periodically using, for example, gas chromatography and the mass fractions of each component calculated based on the concentration of some internal standard [4, 20]. However in this work we use concentration-time curves to determine the concentration of each reacting component and products at any given time. We also present results showing the effect of the reaction rate constants on the main product, biodiesel, and the effect of the concentration of the main reactant, methanol, on the concentration of biodiesel produced. The accuracy of the SRM in this work was validated using ode45, a fourth and fifth order Runge-Kutta scheme whose tolerance level was set to be seven decimal digits (10-7). A further comparison with similar results in the literature was made. Figures 1, 2, 3, and 4 show the conversion of the reaction and yield of biodiesel and the other by-products with reaction time. In all figures, the solid line represents the ode45 solution while the filled circles, squares, and diamonds represent the SRM solutions. Convergence of the SRM was reached after the 15th iteration when N=80  in all simulations. Figure 1 shows the variation of [BD(t)] with reaction rate constant k1 at different concentrations of MeOH. The concentration of TG was taken to be 0.25 and k2=0.2 while k3=0.13. The concentration of BD produced shows to be directly proportional to k1 as it increases with increase in k1. The concentration of BD produced doubles when [MeOH] is also doubled. The SRM solutions agree with the ode45 results.Concentration-time curves showing the effect ofk1 on the concentration of biodiesel produced for different concentrations of MeOH. The solid lines represent ode45 solutions while the figures represent SRM solutions for k2=0.20, k3=0.13, [TG]=0.25.
(a)
(b)The effect ofk2 on the concentration-time curves for biodiesel produced at different concentrations of MeOH. The solid lines represent ode45 solutions while the figures represent SRM solutions for  k1=0.10, k3=0.13, [TG]=0.25.
(a)
(b)Concentration-time curves forBD at different concentrations of MeOH and that of  TG and MeOH when k1=0.10, k2=0.20, k3=0.13. The solid lines represent ode45 solutions while the figures represent SRM solutions.
(a)
(b)Figure 4
The variation of[DG], [MG], and [GL] with time when [TG]=[MeOH]=0.25, k1=0.10, k2=0.20, k3=0.13. The solid lines denote ode45 solutions while the figures represent SRM solutions.In Figure2, the effect of k2 is shown to decrease [BD(t)]. In the figure, k1=0.10, k3=0.13, and [TG]=0.25. The concentration of biodiesel produced in this case is shown to be equal to that of the methanol used. A good agreement between the two solutions is observed.The variation of[BD(t)] with the concentration of the main reactant, MeOH, is shown in Figure 3. Increasing [MeOH] increases that of  BD produced. In the same figure, it is shown that when k1=0.10,k2=0.20, and k3=0.13, about 40% of [TG] gets used up in the reaction, while methanol gets exhausted after about 120 days. The results here are consistent with the analytical results in our earlier study [21]. The two sets of solutions continue to show a very good agreement.Concentration-time curves for the three by-products are shown in Figure4. Di-glyceride is produced in stronger concentrations followed by ML and lastly DG.


## 5. Conclusion

A spectral relaxation method has been used successfully in this work to solve a system of differential equations governing the thermal transesterification reactions in the production of biodiesel from vegetable oils. The method proved useful in solving initial value problems of the type in this study. The results were consistent with numerical results obtained usingode45and those in the works of Sibanda et al. [21]. In the study we investigated the effect of the reaction rate constants and that of the concentration of methanol as the main reactant on the final concentration of biodiesel produced. For optimal results, the study showed that the first reaction rate constant k1 should be kept large while the value of k2 should be small. The results showed that that is increasing k1 increases biodiesel production while increasing k2 decreases biodiesel production. The amount of biodiesel produced increases with increase in the initial concentration of methanol used. Methanol is used up in the reaction after sometime. Experimental research is costly compared to theoretical studies. It is however essential to use methods that give accurate results for trustworthy findings. The spectral relaxation method promises to be a viable tool for theoretical studies in the biosciences.



---

*Source: 601374-2014-01-23.xml*