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}).
Bioethanolbased 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 CO_{2} 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 biocarbon (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.
The production of ethylene from ethanol dehydration occurs via acid catalysis, using zeolites (HZSM5), 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 HZSM5, 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 oligomerizationcracking 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 watersoluble 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 polymergrade 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 ethanol 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 pseudohomogeneous nonisothermal model for a benchscale 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 nonisothermal model for a pilotscale reactor. The twodimensional 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: pseudohomogeneous system; onedimensional system with axial mass and thermal dispersions; effective mass diffusion coefficient and thermal conductivity based on semiempirical 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 steadystate 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.
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 differentialalgebraic equations:
1. Overall mass balance
2. Component mass balance
3. Energy balance
4. Momentum balance
5. Equation of state
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 (D_{Mref}),which is defined by:
The reference effective coefficient of thermal conductivity (k_{Href}) is defined by:
The reference molar mass of the mixture (M_{Mref}) is defined by:
where (Y_{iref}) is the reference mole fraction of the ith component.
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:
where μ_{i} is the viscosity of the ith component (Pa.s).
By specifying the feed pressure (P_{f}), temperature (T_{f}) and composition (y_{if}), the specific mass of the mixture in the feed (ρ_{Mf}) is also specified. Finally, the feed velocity can be obtained by setting the mass flow rate (v_{f}), considered constant in this work, i.e.:
where _{AT} 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 jth reaction (ΔH_{j}), are obtained by the following expressions:
where and , _{TO} are the heats of reaction at temperature T and the reference temperature T_{O}, respectively, and α_{i}(J mol^{1} K^{1}), β_{i}(J mol^{1} K^{2}), γ_{i}(J mol^{1} K^{3}) and δ_{i}(J mol^{1} K^{4}) are constants characteristic of the ith component (^{Reid et al., 1987}).
The equivalent diameter of particles (D_{p}) is the diameter of a perfect sphere with the same volume as the catalyst particles, therefore:
where D_{Cat} is the catalyst diameter and L_{Cat} 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 D_{t} 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 powerlaw type models. Therefore, the model obtained for each of the five reactions can be generalized as follows:
where r_{j} is the dimensionless reaction rate (j=1,...,5), K_{j}C^{ºvj} = k_{j,D} / k_{j,R}, k_{j}, is the specific reaction rate with subscripts D and R indicating the direct and reverse reactions, respectively, Cº is the standard molar concentration (mol m^{3}) of ideal gas at the reference pressure, Pº (1 bar) and reference temperature, Tº (298.15 K), y_{i} is the dimensionless molar concentration of each component (i=1,...,7), ν_{i, j} represents the stoichiometric coefficient of the ith component in the jth reaction, and ν_{j} is the overall reaction order. The chemical equilibrium constant, K_{j}, is defined by:
where , _{Tº} and , _{Tº} are the standard enthalpy of formation and the standard Gibbs energy at the reference temperature, respectively, and is the variation of the heat capacity in each reaction at constant pressure, given by:
The specific reaction rate, K_{j,D}, follows the Arrhenius law, and can be expressed by:
where γ_{j} is the dimensionless activation energy of the jth reaction, θ is the dimensionless temperature, is the mass concentration of the catalyst in kg_{cat} m^{3} and k_{0j} is the preexponential 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_{cat} is the mass of catalyst, D_{T} is the diameter of the bed and 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 preexponential factor of the reactions (K_{0j}) , the effective factor of mass diffusivity (D_{0M}) and the effective factor of thermal conductivity (k_{0H}). The estimation was carried out with the software Mathematica 10 using the method of variable metric LevenbergMarquardt 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.
Parameters  Estimate  Standard deviation 

k_{01}  4.54x10^{2}  4.34x10^{2} 
k_{02}  6.43x10^{1}  6.85x10^{4} 
k_{03}  2.39x10^{3}  1.64x10^{2} 
k_{04}  2.19x10^{6}  3.85x10^{2} 
k_{05}  6.85x10^{3}  4.98x10^{2} 
D_{0M}  7.56 x10^{3}  8.83x10^{4} 
k_{0H}  4.25x10^{1}  6.47x10^{4} 
U_{T}  2.36x10^{1}  2.61x10^{3} 
The final value for the leastsquares objective function was 4.76. The variance of prediction was 3.1014x10^{6}, 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 (R_{2}) 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, x_{k} = κ/ν, η = ΔX_{κ} = X_{κ+1}  X_{κ}, for κ=0,1,...,ν. In each inner point κ, for a generic dependent variable Ψ, the following approximations were applied:
for κ = 1,...,ν1.
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 versus obtained by the approximation of second order; the approximation of the second derivative by fourth order shows an error of versus obtained by approximation of second 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 variablestep BDF method of the SUNDIALS package, with absolute tolerance of 10_{6} and relative 10_{9} (^{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 (t = 5L/v_{ref}). This time used was necessary to ensure that all state variables had their temporal derivative, in any axial position, less than 10_{5} 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.
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 R_{2} and R_{4} far outweigh the reactions R_{1}, R_{3} and R_{5}, 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 Da_{2} is greater than Da_{1} 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 R_{4} that is virtually irreversible, and the intramolecular dehydration of bioethanol, given by the reaction rate R_{1}. 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 R_{2}, highly exothermic for the direct reaction. Analyzing the endothermic nature of reactions R_{1}, R_{3} and R_{4} and the exothermic nature of reactions R_{2} and R_{5} 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 R_{1} is greater than the activation energy of reaction R_{2} so that a higher temperature will favor the reduction of the effects of this difference, favoring reaction R_{1}.
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 R_{2} and R_{1} present considerable reversibility, which is an important observation, because in these reactions ethanol is converted into less interesting byproducts. In addition, the reversibility of reactions R_{1} and R_{4}, 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.
Model (Kagyrmanova et al.) 
Data (Kagyrmanova et al.) 
Proposed Model 


χ_{C2H5OH} (%)  98.5  98.9  98.7 
S_{1/2,4,5,7} (%)  98.8  98.6  98.5 
S_{2/2,4,5,7} (%)  0.21  0.11  0.21 
S_{3/2,4,5,7} (%)  0.26  0.17  0.18 
S_{4/2,4,5,7} (%)  0.90  1.10  1.10 
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.
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 twodimensional 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.49 s). The computer used in this work had the followings specifications: Intel Core i74770 3.4 GHz, 8GB RAM, OS W8 64bit.
CONCLUSION
A mathematical model of a fixedbed catalytic reactor was developed for the process of dehydration of ethanol to ethylene, represented by a dynamic pseudohomogeneous onedimensional model. Effects were considered of the axial dispersion of mass and heat, as well as the reversibility of chemical reactions.
The discretization of the axial variable by a fourthorder centered finite differences formula presented a maximum error of when using 100 points to achieve the required accuracy. The integration package SUNDIALS was appropriate and provided effective results.
The proposed model exhibited better adherence to the experimental data of a pilot plant when compared with the model available in the literature. The model can also be used to describe and optimize operations in industrialscale plants.
NOMENCLATURE
A  Area (m^{2}) 
C  Molar concentration (mol m^{3}) 
C  Heat capacity (J kg^{1}K^{1}) 
D  Diameter (m) 
Ea  Activation energy (J mol^{1}) 
g  Acceleration of gravity (m s^{2}) 
G  Gibb's molar energy (J mol^{1}) 
H  Molar enthalpy (J mol^{1}) 
K  Chemical equilibrium constant () 
L  Length (m) 
m  Mass (kg) 
M  Molar mass (kg mol^{1}) 
P  Total pressure (Pa) 
R  Universal gas constant (J mol^{1}K^{1}) 
r  Dimensionless rate of reaction () 
t  Time (s) 
T  Temperature (K) 
U  Heat exchange coefficient (W m^{2}K^{1}) 
v  Velocity (m s^{1}) 
x  Dimensionless axial coordinate () 
y  Dimensionless molar concentration () 
z  Axial coordinate (m) 
Greek Symbols  
ε  Bed porosity () 
η  Interval discretization () 
Φ  Dimensionless specific mass () 
κ  Point of discretization () 
μ  Dynamic viscosity (Pa s) 
ν  Stoichiometric number () 
ω  Dimensionless velocity () 
Π  Dimensionless Pressure () 
Ψ  Generic variable ( ) 
ρ  Specific mass (kg m^{3}) 
σ  Standard deviation ( ) 
τ  Dimensionless time () 
θ  Dimensionless temperature () 