Brazilian Journal of Chemical Engineering MODELING AND SIMULATION OF THE PROCESS OF DEHYDRATION OF BIOETHANOL TO ETHYLENE

The use of carbon-based waste biomass in the production of plastics can partially meet the growing demand for plastics in the near future. An interest in the production of ethylene from bioethanol has been renewing, motivated mainly by environmental appeal and economics. The main objective of this work is the development of a mathematical model for simulation and optimization of the production of ethylene by the dehydration of ethanol, improving the performance of the process. The phenomenological model proposed is based on mass, momentum and energy balances for the process. The results obtained are satisfactory in comparison with theoretical results and experimental data found in the literature.


INTRODUCTION
There are many renewable raw materials for the production of plastics, but only a limited number of petrochemical products could be produced from biomass through economically viable commercial technologies (Ferreira et al., 2009).Regarding the market of plastic and processing of biomass, the use of bioethanol is a reasonable goal in the medium and long term in order to meet the needs of these industries.Moreover, ethylene is the raw material most often used by the plastics industry and can be produced by available technologies, such as the dehydration of ethanol (Morschbacker, 2009).
Bioethanol-based polyethylene received a certificate as a renewable source by attesting the nature of renewable "green" plastic.In fact, all carbon atoms of ethylene, which end up forming polyethylene, come from CO2 in the atmosphere, which was fixed by sugar cane during its growth in the field.Therefore, the carbon isotopic composition of polyethylene is identical to the one found in the atmosphere: the bioethanolbased polyethylene contains a small amount -approximately 1.2 parts per trillion -of unstable isotope 14 C, which is continuously formed in the atmosphere by the action of cosmic rays.In comparison, the polyethylene obtained from oil or gas -petrochemical polyethylene -practically does not contain the isotope 14 C, since the raw materials were stored under the ground long enough for all unstable isotopes to decay (Carmo et al., 2012).This property allows the determination of the bio-carbon (biobased content) of the bioplastic using standardized methods (ASTM, 2011).
In addition to the environmental benefits, the polymer obtained from green ethylene has the same properties as that obtained from ethylene of fossil origin, allowing the processing industries to leverage their structure to process resin from a renewable source.

Brazilian Journal of Chemical Engineering
The production of ethylene from ethanol dehydration occurs via acid catalysis, using zeolites (HZSM-5), aluminum silicates and aluminum silicates modified with transition metals.The operating temperatures are between 180 and 500 °C with conversions and selectivities depending on the operating conditions.A predominant use of alumina in industrial processes is reported in the literature, probably due to its greater stability (Zhang et al., 2008).The studies in the literature on the transformation of bioethanol into olefins give emphasis to the selective production of ethylene by dehydration, using mainly zeolites HZSM-5, which are effective at temperatures below 300 °C, and catalyst modifications to moderate the strength of the acid sites to avoid secondary reactions of conversion of ethylene and to mitigate the formation of coke.The production of propylene occurs by the conversion of ethylene through a mechanism of oligomerization-cracking that requires temperatures above 350 °C, favoring also the reactions of coke formation and the consequently catalyst deactivation (Sheng et al., 2013).The intramolecular ethanol dehydration, which produces ethylene, is an endothermic and reversible reaction, and the intermolecular dehydration, which produces diethyl ether, is an exothermic and reversible reaction.In addition, ethylene formation by diethyl ether dehydration is an endothermic and reversible reaction.It is widely reported in the literature that, at low temperatures, the production of diethyl ether is predominant and, at higher temperatures, the formation of ethylene is favored.The evaluation of the main and secondary reactions occurring in the process of ethanol dehydration and their operating conditions are fundamental to the understanding of this process and to the proposal of improvements in the production of ethylene (Ramesh et al., 2012).
The process of dehydration of ethanol can operate in an isothermal mode (using a fluid heating) or in an adiabatic mode (using steam dilution).The effluent from the reactor is cooled to remove the water in a condensation column.The crude ethylene that exits at the top of this column is washed to remove acids and other water-soluble components and then passes through a drying bed, forming ethylene of high purity.The removal of the remaining impurities in the ethylene stream is done by distillation columns, producing a polymer-grade ethylene, which is sent to the polymerization plants (Morschbacker, 2009).
In September of 2010, a plant of Braskem S.A. located at Triunfo (RS, Brazil) was started up, being the first factory of green ethylene in the world, which uses 100% renewable raw materials, processing etha-nol from sugar cane.This plant has a production capacity of 200 thousand tons of ethylene per year, consuming about 500 million liters of bioethanol per year, and is responsible for the production of a wide range of polyethylene, to meet the growing demand for increasingly sustainable products.
Nevertheless, there are few works in the literature about the mathematical modeling of reactors for the dehydration of bioethanol.Golay et al. (1999) proposed a pseudo-homogeneous non-isothermal model for a bench-scale reactor.The model with five reactions and four chemical species was only qualitatively able to describe the experimental data.Kagyrmanova et al. (2011) presented a heterogeneous non-isothermal model for a pilot-scale reactor.The two-dimensional model with mass and thermal dispersion included five irreversible reactions and seven chemical species.However, weak assumptions such as constant axial velocity, constant pressure, constant density, and average values of thermodynamic and fluid dynamic properties, impaired the results compared with experimental data.
In this paper, a phenomenological mathematical model of a fixed bed catalytic reactor for the dehydration of bioethanol was proposed and implemented in an efficient computational code, which is able to simulate industrial or pilot plants for ethylene production.The mathematical model developed is based on mass, momentum and energy balances, and reaction rate expressions from the literature.The implementation of the model equations was made in the computational software Mathematica 10.

PROCESS MODELING
The main simplifying assumptions for the construction of the model are: pseudo-homogeneous system; one-dimensional system with axial mass and thermal dispersions; effective mass diffusion coefficient and thermal conductivity based on semi-empirical correlations; constant coolant temperature; constant catalyst activity; negligible heat dissipation by viscous forces and by mass diffusion; Newtonian fluid; uniform porosity; and dynamic viscosity as a function of the mixture composition.
An additional hypothesis considered for the reactor simulation is the quasi steady-state assumption (QSSA) for overall mass balance and for the momentum balance.Thus, it is assumed that the specific mass and the axial velocity have instantaneous response when the system is subjected to disturbances.
Brazilian Journal of Chemical Engineering Vol. 33, No. 03, pp. 479 -490, July -September, 2016 Based on the above simplifying assumptions, the principles of conservation of mass, energy and momentum when applied to the dynamic system in study, are presented in the dimensionless form, using the following definitions: resulting in the system of differential-algebraic equations:  Overall mass balance The initial conditions of these equations obey the following relations: The boundary conditions for concentration and temperature, at the inlet and outlet of the reactor, are based on the work of Langmuir (1908), and more intuitively by Danckwerts (1953), whose rigorous deduction was performed by Bischoff (1961) and are given by: The boundary conditions for the average axial velocity and pressure are specified at the inlet stream and are given by: The dimensionless parameters of the model are presented in Table 1.
Some definitions used in the dimensionless numbers of Table 1 deserve attention, such as the reference coefficient of effective mass dispersion ref M (D ) ,which is defined by: The reference effective coefficient of thermal conductivity The reference molar mass of the mixture   where iref Y is the reference mole fraction of the i-th component.Eq. of Blake-Kozeny BK The viscosity of the mixture   M μ for gases is calculated using a mixing rule based on the work of Wilke (1950), given by an average weighted by molar fraction, as follows:  are the heats of reaction at temperature T and the reference temperature T , respectively, and   are constants characteristic of the i-th component (Reid et al., 1987).
The equivalent diameter of particles   p D is the diameter of a perfect sphere with the same volume as the catalyst particles, therefore: where Cat D is the catalyst diameter and Cat L is the catalyst length.The bed void fraction   ε is obtained via the empirical equation of Haugley and Beveridge (Froment and Bischoff, 1990), expressed as: in which t D is the bed diameter.

KINETIC MODELING
The reactions considered in the process of dehydration of ethanol were based on the kinetic model presented in the work of Kagyrmanova et al. (2011).The experimental studies in their work, covering a temperature range of 350 -450 °C; total pressure of 1 atm; and residence times between 0.005 s and 2.2 s, revealed that the main products of the dehydration of ethanol on the surface of the catalyst based on alumina are ethylene, diethyl ether, acetaldehyde, hydrogen and unsaturated butenes in accordance with the following scheme of reaction: Kagyrmanova et al. (2011) did not adopt the reversibility of the reactions mentioned above; however, these phenomena are considered in the present work and the results show the importance of this consideration.For the sake of notation, the following reference is used: The kinetic model for this system is based on the law of mass action, which means that these reactions are considered elementary and their orders are established by their stoichiometric coefficients, generating power-law type models.Therefore, the model obtained for each of the five reactions can be generalized as follows: where r is the dimensionless reaction rate (j 1, ,5), , j k is the specific reaction rate with subscripts D and R indicating the direct and reverse reactions, respectively, C is the standard molar concentration   3 mol m  of ideal gas at the reference pressure, P (1 bar) and reference temperature, T (298.15K), i y is the dimensionless molar concentration of each component (i 1, ,7)   , i, j ν represents the stoichiometric coefficient of the i-th component in the j-th reaction, and j ν is the overall reaction order.The chemical equilibrium constant, j K , is defined by: T T j, T j, T j, T j P , j P , j ΔC  is the variation of the heat capacity in each reaction at constant pressure, given by: The specific reaction rate, j,D k , follows the Arrhenius law, and can be expressed by: where j  is the dimensionless activation energy of the j-th reaction,  is the dimensionless temperature, cat Ĉ is the mass concentration of the catalyst in 3 cat kg m  and 0 j k is the pre-exponential factor, whose units depend on the reaction.The mass concentration of catalyst is defined in terms of the mass of catalyst and the useful volume of the bed, i.e.: where m is the mass of catalyst, T D is the diameter of the bed and L is the axial length of the bed.

RESULTS AND DISCUSSION
Reactor data and operating conditions were based on the work of Kagyrmanova et al. (2011), and the information of the components was obtained from the thermodynamic database in Reid et al. (1987).

Estimation of Parameters
The experimental data required for parameter estimation are derived from a pilot plant described in the experimental work of Kagyrmanova et al. (2011).These data are the temperatures at different axial positions at the steady state of the plant and were used in the parameter estimation problem.
The estimated parameters were the pre-exponential factor of the reactions   with the software Mathematica 10 using the method of variable metric Levenberg-Marquardt with a 95% confidence level.
There was no statistical representation of the experimental data, because there were no replicates.Table 2 shows the values obtained for the parameters and their standard deviations.
The final value for the least-squares objective function was 4.76.The variance of prediction   2  was 6 3.1014x10  , statistically equal to the sample variance, which is zero (no replicates).Therefore, the model can be considered to be satisfactory and there is no apparent reason to be discarded (nor the possibility of a perfect model).Thus, the model can represent the experimental data satisfactorily.This means that the errors of prediction are not significantly higher than the experimental errors considered and that the model does not take better results than the data used to generate it.Therefore, it is not necessary to exert extra effort to refine the model and, conceptually, the experimental errors are not underestimated or overestimated.The determination coefficient   2 R obtained in the estimation process was 0.999998, which represents a great fit of the model to the experimental data.

Process Simulation
In this work, the discretization of the reactor in the axial direction was done using a fourth order central finite differences formula.The problem domain is divided into ν equal segments, k for a generic dependent variable ψ , the following approximations were applied: At the reactor boundaries ( x 0  and x 1  ), the boundary conditions are valid and the following approximations were applied: The magnitude of the error of fourth order is markedly lower than the error of second order.The approximation of the first derivative of fourth order shows an error of   order (Fornberg, 1998).The number of discretization points   ν used in this work was equal to 100, as a result of a mesh convergence analysis.The time step of integration was calculated from a variable-step BDF method of the SUNDIALS package, with absolute tolerance of 6 10  and relative 9 10  (Hindmarsh and Taylor, 1999).
The steady state was obtained through the use of the dynamic model for a time equal to five times the residence time ref (t 5L / v )  .This time used was necessary to ensure that all state variables had their temporal derivative, in any axial position, less than    It can be observed in Figure 1 that the ethanol is not fully consumed, reaching 98.8% of the value of its mole fraction in the equilibrium, even with the temperature at a level which is favorable for the ethanol dehydration reaction, which can be verified in Figure 2.This evidence suggests that a reactor with a larger size could be used to increase the conversion of ethanol or to reach the concentration of chemical equilibrium of ethanol, not taking in account the technical and economic aspects, like the pressure drop and the operating costs.In Figure 1, it is shown that diethyl ether is quickly formed in large quantity through intramolecular dehydration of ethanol, but suffers dehydration to ethylene, reaching negligible values of mole fraction at the end of the process.The unsaturated butenes only present a considerable amount when the concentration of ethylene is significant.Acetaldehyde and hydrogen have the same molar composition, which was expected, because they are formed in the same proportion by the same reaction and have null initial concentrations in the reactor.
By means of the analysis of Figure 3, it is possible to observe that the magnitudes of reactions 2 R and 4 R far outweigh the reactions 1 R , 3 R and 5 R , which is confirmed by observing Figure 1, which is in agreement with the literature.At the beginning of the reactor bed, the formation of ethyl ether is greater than that of ethylene, which is explained on the basis of the comparison between the values of the Damköhler numbers of the first and second reaction, in which 2 Da is greater than 1 Da in almost the whole bed.However, the formation of the main product becomes greater than ethyl ether after the first 1% of the length of the bed, an observation justified by several factors, among them: by the decomposition of this byproduct into ethylene, reaction 4 R that is virtually irreversible, and the intramolecular dehydration of bioethanol, given by the reaction rate 1 R .The operating conditions of the simulation also contribute to this effect, due to the high temperature, low ethanol concentration and ordinary total pressure, characteristics that favor the direct, endothermic and first order main reaction.
A peak temperature at the beginning of the reactor bed is observed, which is due to the high reaction rate 2 R , highly exothermic for the direct reaction.Ana- lyzing the endothermic nature of reactions 1 R , 3 R and 4 R and the exothermic nature of reactions 2 R and 5 R Brazilian Journal of Chemical Engineering Vol. 33, No. 03, pp. 479 -490, July -September, 2016 and knowing the decrease in temperature at the beginning of the reactor, which can be seen in Figure 2, there is an overall character of an endothermic reaction system.However, after this, the ethanol concentration remains lower, decreasing the reaction rates.In addition, the temperature mainly increases due to the effect of thermal exchange with the heating fluid and the formation of butenes.The activation energy of reaction 1 R is greater than the activation energy of reaction 2 R , so that a higher temperature will favor the reduction of the effects of this difference, favoring reaction 1 R .The reduction of the total pressure shows the volumetric expansion of the reaction mixture during the process.This expansion can also be checked by analyzing Figure 2, in which the specific mass decreases.The process presents overall molar variation resulting from the system of chemical reactions and, verifying that the axial velocity increases, it is possible to affirm that there is a significant increase in the number of moles.In Figure 2, it is also possible to note an almost linear profile for the pressure, which is common in problems with relatively low velocity.Looking at Figure 2, the pressure drop is found to be moderate, around 0.5%.This fact indicates that the pressure gradient in the energy balance has little influence on the temperature profile.
In Figure 2, the profiles of specific mass of the mixture and axial velocity have inverse behavior, which is expected because the product between these two variables must be constant throughout the bed, showing that there is not mass accumulation in the bed, in accordance with the model assumption.The specific mass of the mixture reduces by 43%, while the axial velocity increases 76%, which corroborate the justification of these variables not being considered to be constant in the process.
Figure 3 shows that only reactions 2 R and 3 R present considerable reversibility, which is an important observation, because in these reactions ethanol is converted into less interesting by-products.In addition, the reversibility of reactions 1 R and 4 R , where ethylene is formed, is negligible.
The Peclet thermal number, presented in Table 3, was calculated with a magnitude around 60, justifying the consideration of the thermal conductivity effect in the system analysis.However, the value found for the Peclet mass number (around 300) presents evidence that the effects of mass diffusive nature could be disregarded in this system.This means that the advection effects may be predominant in mass transport.However, the conversion of ethanol at the inlet of the reactor is greater than 20%, indicating the importance of mass diffusion effects and the relevance of considering this term in the model.The dimensionless heat exchange coefficient, illustrated in Table 3, shows an increase over the reactor due to the reduction of specific heat capacity of the reactive mixture.This indicates that a large amount of energy is supplied by the reactor's jacket along the reactor axial direction.Also illustrated in Table 3 is the increase in value of the reactor dimensionless group  , which suggests an increase in the contribution of mechanical energy in relation to the thermal.
Comparison with the Literature Kagyrmanova et al. (2011) carried out experimental work on a pilot plant in order to compare with the theoretical computational results obtained in the simulation of their proposed model.The experimental data presented in this reference are the axial profile of temperature at the steady state, in addition to the conversion of ethanol and the selectivities of certain components.For the temperature profile there were 14 experimental points, one of them being the measure of the inlet temperature.
Figure 4 illustrates the comparison of the axial temperature profile between the model proposed in this work, the model proposed by Kagyrmanova et al. (2011) and the experimental data.
Examining Figure 4, it is noted that the proposed model gives a better fit to the experimental data, especially near the end of the reactor.In order to quantify these differences, Table 4 presents the values of the residues (the difference between calculated and measured values) for each experimental point, and also the sum of square residuals for both models.
The absolute values of the residues found with the proposed model are at all points lower than those obtained with the model used by Kagyrmanova et al. (2011), with one exception for the temperature at position 0.5 m, demonstrating the better quality of the model proposed in the present work.The values for the conversion of ethanol and the selectivity of the products are shown in Table 5.Once again, the model proposed exhibited better adherence to the experimental data than the model proposed in the work of Kagyrmanova et al. (2011), possibly due to the different simplifying assumptions.

Brazilian Journal of Chemical Engineering
Therefore, it seems that the reversibility of the chemical reactions is effective in the kinetic scheme.Also, the variation of total pressure, axial velocity and specific mass of the mixture significantly influence the profiles, and the functional dependencies of the parameters with the state variables are relevant.Finally, it is possible to note that the radial dispersion of mass and heat is negligible due to the small diameter of the reactor and the resistance to intraparticle mass and heat transfer and interfacial tension can be neglected.
In this way, the theoretical system presented by Kagyrmanova et al. (2011) shows some complexity of a heterogeneous two-dimensional model without resulting in a better compliance with the experiments.Moreover, the simplifying assumptions adopted by these authors are strongly dissonant with the real system.
The simulation took about 8 seconds to solve the dynamic problem with the end time of simulation equal to 5 times the dimensionless residence time (19.49s).The computer used in this work had the followings specifications: Intel Core i7-4770 3.4 GHz, 8GB RAM, OS W8 64-bit.

CONCLUSION
A mathematical model of a fixed-bed catalytic reactor was developed for the process of dehydration of ethanol to ethylene, represented by a dynamic pseudo-homogeneous one-dimensional model.Effects were considered of the axial dispersion of mass and heat, as well as the reversibility of chemical reactions.
where i μ is the viscosity of the i-th component   Pa.s .By specifying the feed pressure   specific mass of the mixture in the feed   Mf ρ is also specified.Finally, the feed velocity v can be obtained by setting the mass flow rate   m  , considered constant in this work, i.e.: where A is the cross section of the tubular reactor.The thermodynamic properties, such as specific heat capacity of the mixture   P,M Ĉ and the standard enthalpy of each j-th reaction   j ΔH , are obtained by the following expressions:

15 10 
characterizing numerically the achievement of steady state.At this time instant, the state variables in the whole reactor were evaluated   0 z L   .Figures 1 to 3 show the spatial profiles at the steady state of the main dependent variables and some dimensionless parameters characteristic of the system.

Figure 1 :
Figure 1: Stationary axial profile of mole fractions of the components.

Figure 2 :
Figure 2: Stationary axial profile of temperature, pressure, axial velocity and specific mass of the mixture.

Figure 3 :
Figure 3: Stationary axial profile for the dimensionless parameters (Damköhler number and chemical equilibrium constant).

Figure 4 :
Figure 4: Comparison between the proposed model and the experimental and theoretical results of the literature for the axial temperature profile at the steady state.