Steady-state modeling of reactive distillation columns

An algorithm for the solution of the mathematical model featuring reactive distillation process in steady-state columns is analyzed. It has been presumed that each stage’s outlet streams in the model were in phase equilibrium conditions and that chemical kinetics was described by a reaction rate model. Within the context of the developed algorithm a procedure to solve a set of equations in a sequential form and a methodology to produce the initial estimates was defined. Broyden’s method was employed to solve the equations that model the chemical reactions. Algorithm was evaluated by study cases of 2-pentene metathesis and MTBE synthesis. The simulation results were close to results available in the literature and the proposed algorithm proved to be efficient since in both cases the convergence towards a solution was found.


Introduction
The analysis of a plant, by simulation, within the development of new processes, may frequently show beforehand whether it is technically and economically viable.The simulation process in already operating plants may optimize the operational conditions for better quality products, decrease energy consumption and other losses in the process (MARQUINI et al., 2007).
Several experimental situations may be simulated through modeling and final results are then tested with real experiments (KROUMOV et al., 2007).
Reactive distillation is a process in which a chemical reaction and separation by distillation continuously occur within a single operation (KHALEDI; BISHNOI, 2006).Reactive distillation may be an alternative to conventional processes.Distillation and reaction simultaneously are only possible if the conditions of both operations are combined, or rather, the reaction must have a high conversion in temperature and pressure levels compatible to the distillation column's operational conditions.
Reactive distillation shows characteristics that makes the process highly attractive: (i) process simplification decreases costs with equipments; (ii) greater conversions may be obtained for chemical equilibrium-limited reactions with cost reduction through the recycling of excess reagents; (iii) higher selectivity is obtained; (iv) the employment of smaller amount of catalysts for the same conversion degree; (v) reaction conditions may avoid possible azeotropes of the reaction's product mixture and thus avoid the use of auxiliary solvents; (vi) benefits of heat integration are obtained by the heat produced in the chemical reaction for vaporization and, consequently, hot spots are avoided (TAYLOR; KRISHNA, 2000).
The complex and non-linear nature of several issues in Chemical Engineering is frequently described by phenomenological or empirical models represented by non-linear algebraic equations (CORAZZA et al., 2008).The mathematical models of steady-state reactive distillation involve the solution of a system of non-linear equations.The equations' nonlinearity is mainly produced from equations used in the modeling of phase equilibrium and from equations used for the calculation of the extent of reaction.Newton-Raphson method and its adaptations are rather common to solve these problems, although the methods may fail when good initial guess are lacking or more than one solution is possible.Current research provides a methodology for initial guess.The equations that model the chemical reactions will be solved by modified Newton-Raphson method as suggested by Broyden (1965).
In spite of their great importance in the simulation of reactive distillation, initial guess are not well represented in the literature.Research by Baur et al. (2000), Fernholz et al. (2000), Qi et al. (2004) and Cheng and Yu (2005) showed simulations of reactive distillation columns but failed to forward a methodology to obtain the initial estimates.Alfradique and Castier (2005) simulated a reactive distillation column and solved all the equations of the model simultaneously by the Newton-Raphson method.Initial guess used by the above authors comprised linear temperature profiles and results available in the literature.Khaledi and Bishnoi (2006) developed a repeatable algorithm for simulations, in a stead state, for reactive distillation processes, limited by chemical and phase equilibrium, while presuming as initial guess the linear profile of temperature and pressure, considered a nonreactive system for molar streams flow rate.In the catalyst distillation of butyl acetate synthesis, Gangadwala and Kienle (2007) forwarded a study using improved MINLP.As initial guess, the authors used a solution of a model which took into consideration the equimolar flow between the phases and solved the problem by employing results taken from the literature.
Scarcity of investigations on a detailed methodology to produce initial guess demonstrates the need for the development of research on the matter.

Modeling
The mass and energy balances was undertaken during the development of the model.It was presumed that outlet streams in each stage (L j and V j ) were in phase equilibrium conditions.Reaction rate equations were used to calculate the extent of reaction.Molar holdup and pressure were presumed constant.
Figure 1a.shows the schematic representation of all the stages of a distillation column.The algorithm given by Seader and Henley (1998) for conventional distillation columns was modified for reactive distillation columns and incorporated the equations to calculate the extent of reaction.Figure 1b shows the adapted algorithm, where T is the temperature, ξ is the extent of reaction, B is the bottoms; D is the distillated product; V is the vapor flow; L is the liquid flow; x is the mole fraction of the liquid phase, y is the mole fraction of the vapor phase, h is the enthalpy of the liquid phase, H is the enthalpy of the vapor phase, Q is the stage's heat flow to the neighborhood; and indexes i, j and k represent the component, stage and reaction, respectively.
The linear profile between the reboiler (highest temperature) and the condenser (lowest temperature) is used in the initial guess of temperature at each stage.The highest and the lowest saturation temperatures of pure components are employed in the column's operation pressure, calculated by equations that model the vapor pressure of pure components.
In the case of the extent of reaction, a guess is used that presents a maximum rate within a determined stage and a quadratic rate fall of the same up to the extreme reactive stages (first and last reaction stages).The highest extent of reaction is obtained when feed streams reach the chemical equilibrium in the temperature of feed.The most common feed characteristics of a reactive distillation column should be taken into account to choose the stage which will receive the highest rate of the extent of reaction.The two most common feed forms in a column of reactive distillation are: (1) all feed occurs in a single stage and (2) feed occurs in two stages; in this case, the most volatile reagent is feed close to the reboiler and the less volatile reagent close to the condenser.In the first case, the stage chosen is the feed stage; in the second case, the stage is that in which feed occurs closer to the reboiler; if the feed does not occur in a reactive stage, the closest reactive stage, stage should be chosen.For the initial guess of distillate and bottoms, it is presumed that for each mol of liquid that evaporates another mol of vapor condenses, and vice-versa.Liquid and vapor streams are only affected by feed factors (flow rate and quality), reaction and side streams.When the above is applied to molar global and for the liquid phase balances, with rearrangement of expressions, the following equations are given: Guess B and D , Eqs. ( 1) and ( 2) Guess j L , Eqs. ( 7)-( 9) x , Eqs. ( 10)-( 12) error ≤ tol yes Print results no Normalize i j x , Eq. ( 20) Calculate j T , Eq. ( 22) Calculate i j y , Eq. ( 23) Calculate j h and j H , Eqs. ( 24) and ( 25) 26) and ( 27) Calculate j V , Eqs. ( 28) and ( 32) Calculate B and D , Eqs. ( 35) and (36) Calculate j L , Eqs. ( 37)-( 39) in which: where: F is the feed stream, U is the liquid side stream, W is the vapor side stream, r B is the reboil ratio, r D is the reflux ratio, q is the feed quality stream, v is the stoichiometric coefficient, N C is the number of components, N P is the number of stages and N r is the number of reactions.
When the equimolar flow between the phases is once more taken into account, with molar balances in each stage for the vapor phase, the estimate of the vapor streams are given by the following equations: ( ) Through global molar balances, comprising the stages between the first and each of the intermediate ones, the following equations for guess of liquid streams are obtained: Once calculated the set of guess values, it may be started the algorithm iterative calculations.
Equations of component mass balances in each stage and the condition of phase equilibrium give the following set of linear equations for mole fractions of the liquid phase.
in which: 1 , 2,3 ,..., ( ) The phase equilibrium relation is given as Equations ( 10) -( 12) make up a tri-diagonal linear system with N p unknown values that should be determined for each component.
Since the initial guess imprecision and the mole fractions of each component are calculated separately, especially in the first iterations, they may have values without any physical meaning (such as, negative, a sum different from unity and others).So that the process of algorithm convergence may be accelerated by normalization, undertaken by the following equation: Calculation of temperatures is undertaken by the re-arrangement of the phase equilibrium equation using the secant method to solve the equations: ( ) Mole fractions of the vapor phase are calculated directly by the phase equilibrium equation: , 1, 2,..., ; 1, 2,..., Enthalpy of the liquid and vapor phases are calculated by the elemental reference state (standard enthalpy of formation of pure components in ideal gas state): The following equations, used for the calculation of heat duties in the condenser and reboiler were obtained from the first stage and global energy balances equations respectively: ( ) ( ) The combination of molar balances, definition of reboil ratio and energy balances, with modifications, provide the following set of linear equations: in which: and in which: For the calculation of distillate and bottoms rates, only the reflux and reboil ratios definition and molar balances in the reboiler and condenser are used.
Equations of the molar balances may be rearranged so that the liquid phase streams may be calculated from the following equations: When the reaction is kinetically-limited, the extents of reaction are calculated by the reaction rate equation.When the reaction is very fast, tha is, when the reaction is equilibrium-limited, the extents of reaction may be calculated by the equilibrium equation.This is the most complex step in the algorithm and the step which limits the convergence of the simulation.The use of a method for the solution of a set of non-linear equations is required since most methods demand the calculation of the derivatives of the functions.The derivatives are numerically evaluated due to the difficulty in obtaining algebraic derivatives for this set of equations.Broydon's method was employed for the solution of this algorithm step.
The end of the iterative process occurs when the relative variation of some main variables of the model are very low.The smallness of the variation depends on tolerance (tol) chosen.

Results and discussion
Two examples taken from the literature about reactive distillation were used to evaluate the proposed methodology: 2-pentene metathesis, investigated by Chen et al. (2000), and MTBE synthesis, studied by Singh et al. (2005).Results from the two case studies will be shown in the following section.

Case 1: 2-pentene metathesis
The aim of the column in this particular case is that the metathesis of 2-pentene takes place to form 2-butene and 3-hexene, coupled to the separation of the products.The column's operational conditions were taken from Chen et al. (2000).A 14-stage column was used, was supposed that reaction is Table 1 shows the column's operational conditions.Reboil ratio was adapted so that molar flow rates would be equal for the distillate and bottoms.
Whereas the Peng-Robinson equation of state was employed to represent the behavior of the vapor phase, the groups contribution method UNIFAC was used to represent excess Gibbs's free energy, with parameters from Poling et al. (2004).Physical, critical and thermodynamic properties were taken from Yaws (2003).A code in FORTRAN was developed to solve the case.The program converged with 236 iterations, with 1 x 10 -10 tolerance, using Intel Core2 QUAD Q6600 2,40GHz processor.
Figures 2a and b show comparisons between initial guess and simulation results respectively for temperature and extent of reaction.Initial guess profiles were extremely close to results.
Figures 3a and b show comparisons between initial guess and simulation results respectively for molar flow rates of the liquid and vapor phases.Initial guess of rates of both phases were extremely close to the values calculated by simulation .
Table 2 shows relevant variables in the reactive distillation process obtained in this work, which results is in agreement with reference.Case 2: MTBE synthesis The aim of the column in Case 2 is the etherification of iso-butene (i-B) and methanol (MeOH) in the presence of an inert n-butane (n-B) and strong acid catalyst to produce methyl tert-butyl ether (MTBE) (SINGH et al., 2005;VENIMADHAVAN et al., 1994) and the separation of the product from inert.Data about the column were taken from Singh et al. (2005).A 17-stage column was used and it was supposed that reaction occurs at stages 4 to 11.The MTBE synthesis occurs according to the reaction: The equation for the reaction rate was given by Venimadhavan et al. ( 1994): 3187.0 4464 exp 8.33 10 exp 6820.0 Peng-Robinson equation of state was employed to represent the non-ideal behavior of the vapor phase, whereas Wilson model was used for the nonideal behavior of the liquid phase to represent Gibbs's excess free energy.The parameters of the binary iteration for Wilson model was taken from Ung and Doherty (1995).Physical, critical and thermodynamic properties were taken from Yaws (2003).
Table 3 shows the column operational conditions.
A code in FORTRAN was developed to solve the case.The program converged with 471 iterations, with 1 x 10 -10 tolerance, using Intel Core2 QUAD Q6600 2,40GHz processor.CHEN et al., 2000).
Figure 4a shows the comparison of temperature profiles obtained from initial guess and simulation results.Initial guess were very close to simulation results.Figure 4b shows the comparison between initial guess and simulation results for the extent of reaction.Initial guess in this case were not so close to simulation results, although the profile was similar.Figures 5a and b show comparisons between initial guess and simulation results of rate profiles for liquid and vapor phases.Initial estimates of moles discharge in the two phases are very close to simulation-calculated rates.

Temperature (K)
Table 4 shows relevant variables in the reactive distillation process obtained in this work.
The results obtained in this work for Case 2 are reasonably close to those by Singh et al. (2005).Certain differences may be justified due to some differences in thermodynamic models.The rate of distillate stream values are different, however, the rate provided by Singh et al. (2005) is not correct, because due to the reaction, represented by Equation ( 42), there is a decrease in the number of moles.

Figure 1 .
Figure 1.(a) -Schematic representation of a reactive distillation column; (b) -Algorithm for the simulation of reactive distillation columns operate in a stead state.

Figure 2 .Figure 3 .
Figure 2. Comparison between the simulation results and initial guess in Case 1: (a) -temperature; (b) -extent of reaction.

Figure 4 .
Figure 4. Comparison between simulation results and initial guess in Case 2: (a) -temperature; (b) -extent of reaction.

Figure 5 .
Figure 5.Comparison between simulation results and initial guess in Case 2: (a) -Liquid rate; (a) -Vapor rate.

Table 1 .
Operational conditions of reactive distillation column in Case 1.

Table 2 .
Chen et al. (2000) results obtained in this work and those byChen et al. (2000)for Case 1.

Table 3 .
Operational conditions of distillation column in Case 2.

Table 4 .
Comparison between results obtained in this work and those by Singh et al. (2005) for Case 2.