# Steady-state modeling of reactive distillation columns/Modelagem de colunas de destilacao reativa em regime permanente.

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.

Material and methods

Modeling

The mass and energy balances was undertaken during the development of the model. It was presumed that outlet streams in each stage ([L.sub.j] and [V.sub.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, [xi] 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.

[FIGURE 1 OMITTED]

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

in which:

[R.sub.j] = [[N.sub.r].summation over (t=1)][[N.sub.c].summation over (t=1)] [v.sub.li] x [[xi].sub.lj] (3)

where:

F is the feed stream, U is the liquid side stream, W is the vapor side stream, [r.sub.B] is the reboil ratio, [r.sub.D] is the reflux ratio, q is the feed quality stream, v is the stoichiometric coefficient, [N.sub.C] is the number of components, [N.sub.P] is the number of stages and [N.sub.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:

[V.sup.(0).sub.2] = ([r.sub.D] + 1) x [D.sup.(0)] - [F.sub.1] (4)

[V.sup.(0).sub.j] = [V.sup.(0).sub.j-1] + [W.sub.j-1] + ([q.sub.j-1] - 1) x [F.sub.j-1] (5)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

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:

[L.sup.(0).sub.1] = [r.sub.D] x [D.sup.(0)] (7)

[L.sup.(0).sub.j] = [V.sup.(0).sub.j+1] = [V.sub.1] - D + [j.summation over (m=1)] [F.sub.m] + [j.summation over (m=2)] ([R.sub.m] - [U.sub.m] - [W.sub.m]), j = 2, 3, ..., [N.sub.p.sup.-1] (8)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

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.

[d.sub.1] - [X.sub.i1] + [u.sub.1] [x.sub.12] = [b.sub.1] (10)

[l.sub.j] - [X.sub.ij-1] + [d.sub.j] - [x.sub.ij] + [u.sub.j] - [x.sub.ij+1] = [b.sub.j], j = 2, 3, ..., [N.sub.p.sup.-1] (11)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

in which:

[l.sub.j] = [L.sub.j-1] j = 2, 3, ..., [N.sub.p] (13)

[d.sub.1] = - [[V.sub.1] x [K.sub.i1] + [L.sub.1] + D] (14)

[d.sub.j] = -[([U.sub.j] + [L.sub.j]) + ([W.sub.j] + [V.sub.j]) x [K.sub.ij]], j = 2, 3, ..., [N.sub.p.sup.-1] (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[u.sub.j] = [V.sub.j+1] x [K.sub.ij+1], j = 1, 2, ..., [N.sub.p.sup.-1] (17)

[b.sub.j] = - [F.sub.j] x [z.sub.ij] [[N.sub.r].summation over (l=1)] [v.sub.li] x [[xi].sub.lj], j = 1, 2, ..., [N.sub.p] (18)

The phase equilibrium relation is given as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

Equations (10)-(12) make up a tri-diagonal linear system with [N.sub.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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

Calculation of temperatures is undertaken by the re-arrangement of the phase equilibrium equation using the secant method to solve the equations:

F([T.sub.j]) = [[N.sub.c].summation over (i=1)][K.sub.ij] x [x.sub.ij] = 0, j = 1, 2, ... [N.sub.p] (21)

[T.sup.k+1.sub.j] = [T.sup.k.sub.j] - F([T.sup.k.sub.j]) x [[T.sup.k.sub.j] - [T.sup.k-1.sub.j]/ F([T.sup.k.sub.j]) - F([T.sup.k-1.sub.j]] (22)

Mole fractions of the vapor phase are calculated directly by the phase equilibrium equation:

[y.sub.ij] = [K.sub.ij] x [x.sub.ij], i = 1, 2, ..., [N.sub.c]; j = 1, 2, ..., [N.sub.p] (23)

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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

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:

[Q.sub.C] = [F.sub.1] x [h.sub.Feed 1] + [V.sub.2] x [H.sub.2] - [V.sub.1] x [H.sub.1] - ([L.sub.1] + D) x [h.sub.1] (26)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

The combination of molar balances, definition of reboil ratio and energy balances, with modifications, provide the following set of linear equations:

[[alpha].sub.j] x [V.sub.j] + [[beta].sub.j] x [V.sub.j+1] = [[gamma].sub.j], j = 2, 3, ..., [N.sub.p.sup.-1] (28)

in which:

[[alpha].sub.j] = [h.sub.j-1] - [H.sub.j] (29)

[[beta].sub.j] = [H.sub.j+1-] - [h.sub.j] (30)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

in which:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

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.

D = [V.sub.2] + [F.sub.1]/[r.sub.D] + 1 (36)

Equations of the molar balances may be rearranged so that the liquid phase streams may be calculated from the following equations:

[L.sub.1] = [r.sub.D] x D (37)

[L.sub.j] = [V.sup.j+1] - [V.sub.1] - D + [j.summation over (m=1)] [F.sub.m] + [j.summation over (m=2)] ([R.sub.m] - [U.sub.m] - [W.sub.m]), j = 2, 3, ..., [N.sub.p.sup.-1] (38)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

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.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

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 equilibrium-limited and occurs on stages 2 to 13. The following chemical equation describes the 2-pentene metathesis:

2[C.sub.5][H.sub.10]][] [C.sub.4][H.sub.8] + [C.sub.6][H.sub.12] (41)

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).

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

A code in FORTRAN was developed to solve the case. The program converged with 236 iterations, with 1 x [10.sup.-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 ?so-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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (42)

The equation for the reaction rate was given by Venimadhavan et al. (1994):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

Peng-Robinson equation of state was employed to represent the non-ideal behavior of the vapor phase, whereas Wilson model was used for the non-ideal 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.sup.-10] tolerance, using Intel Core2 QUAD Q6600 2,40GHz processor.

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.

[FIGURE 4 OMITTED]

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.

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 5 OMITTED]

Conclusion

An algorithm was provided in this work for the solution of a reactive distillation column that need solution for non-linear equations system.

Equations in the model were divided into sets and each set was solved separately. A methodology was also provided to produce initial guess which constitutes a critical step in the solution of nonlinear equations system.

The suggested methodology for the production of initial estimates was efficient with values close to those calculated by simulation. This fact accelerated and increased the convergence warranty.

The algorithm used was also efficient since the results obtained with the code developed were very close to those in the literature. There were no convergence problems in the solution of the equation system.

Nomenclature

B Bottoms;

D Distillate;

[D.sub.a] Damkohler number;

[F.sub.j] Feed rate of stage j;

h Enthalpy of liquid phase;

H Enthalpy of vapor phase;

[L.sub.j] Liquid rate of stage j;

[N.sub.c] Total number of components;

[N.sub.p] Total number of stages;

[q.sub.j] Feed quality of stage j;

[Q.sub.j] Heat load from stage j;

[r.sub.B] Reboil ratio;

[r.sub.D] Reflux ratio;

[T.sub.j] Temperature of stage j;

[U.sub.j] Liquid side rate of stage j;

v Volume;

[V.sub.j] Vapor rate of stage j;

[W.sub.j] Lateral vapor outlet of stage j;

[x.sub.ij] Liquid mole fraction of component i on stage j;

[[gamma].sub.ij] Vapor mole fraction of component i on stage j;

[z.sub.ij] Mole fraction of component i on stage j.

Doi: 10.4025/actascitechnol.v34i1.9535

References

ALFRADIQUE, M. F.; CASTIER, M. Modeling and simulation of reactive distillation columns using computer algebra. Computers and Chemical Engineering, v. 29, n. 4, p. 1875-1884, 2005.

BAUR, R.; HIGLER, A. P.; TAYLOR, R.; KRISHNA, R. Comparison of equilibrium stage and nonequilibrium stage models for reactive distillation. Chemical Engineering Journal, v. 76, n. 1, p. 33-47, 2000.

BROYDEN, C. G. A class of methods for solving nonlinear simultaneous equations. Mathematics of Computation, v. 19, n. 92, p. 577-593, 1965.

CHEN, F.; HUSS, R. S.; MALONE, M. F.; DOHERTY, M. F. Simulation of kinetic effects in reactive distillation. Computers and Chemical Engineering, v. 24, n. 11, p. 2457-2472, 2000.

CHENG, Y. C.; YU, C. C. Effects of feed tray locations to the design of distillation and its implication to control. Chemical Engineering Science, v. 60, n. 17, p. 4661-4677, 2005.

CORAZZA, F. C.; OLIVEIRA, J. V.; CORAZZA, M. L. Application of a subdivision algorithm for solving nonlinear algebraic systems. Acta Scientiarum. Technology, v. 30, n. 1, p. 27-38, 2008.

FERNHOLZ, G.; ENGELL, S.; KREUL, L.-U.; GORAK, A. Optimal operation of a semi-batch reactive distillation column. Computers and Chemical Engineering, v. 24, n. 2-7, p. 1569-1575, 2000.

GANGADWALA, J.; KIENLE, A. MINLP optimization of butyl acetate synthesis. Chemical Engineering and Processing, v. 46, n. 2, p. 107-118, 2007.

KHALEDI, R.; BISHNOI, P. R. A method for modeling two- and three-phase reactive distillation columns. Industrial and Engineering Chemical Research, v. 45, n. 17, p. 6007-6020, 2006.

KROUMOV, A. D.; MODENES, A. N.; WENZEL, B. M. Desenvolvimento de um modelo da cinetica enzimatica da transesterificacao de oleos vegetais para producao de biodiesel. Acta Scientiarum. Technology, v. 29, n. 1, p. 9-16, 2007.

MARQUINI, M. F.; MARIANI, D. C.; MEIRELLES, A. J. A.; SANTOS, O. A. A.; JORGE, L. M. M. Simulacao e analise de um sistema industrial de colunas de destilacao de etanol. Acta Scientiarum. Technology, v. 29, n. 1, p. 23-28, 2007.

POLING, B. E.; PRAUSNITZ, J. M.; O'CONNELL, J. P. The properties of gases and liquids. 5th ed. Digital Engineering Library: McGraw-Hill, 2004.

QI, Z.; KIENLE, A.; STEIN, E.; MOHL, K.-D.; TUCHLENSKI, A.; SUNDMACHER, K. MTBE decomposition in a reactive distillation column. Chemical Engineering Research and Design, v. 82, n. 2, p. 185-191, 2004.

SEADER, J. D.; HENLEY, E. J. Separations process principles. New York: John Wiley and Sons, 1998.

SINGH, B. P.; SINGH, R.; KUMAR, M. V. P.; KAISTHA, N. Steady-state analyses for reactive distillation control: An MTBE case study. Journal of Loss Prevention in the Process Industries, v. 18, n. 4-6, p. 283-292, 2005.

TAYLOR, R.; KRISHNA, R. Modelling reactive distillation. Chemical Engineering Science, v. 55, n. 22, p. 5183-5229, 2000.

VENIMADHAVAN, G.; BUZAD, G.; DOHERTY, M. F.; MALONE, M. F. Effect of kinetics on residue curve maps for reactive distillation. American Institute of Chemical Engineers Journal, v. 40, n. 11, p. 1814-1824, 1994.

UNG, S.; DOHERTY, M. F. Vapor-liquid phase equilibrium in system with multiple chemical reactions. Chemical Engineering Science, v. 50, n. 1, p. 23-48, 1995.

YAWS, C. L. Yaws' Handbook of thermodynamics and physical properties of chemical compounds. New York: Knovel, 2003.

Received on March 4, 2010.

Accepted on October 26, 2010.

License information: This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Vilmar Steffen * and Edson Antonio da Silva

Curso de Engenharia Quimica, Universidade Estadual do Oeste do Parana, R. da Faculdade, 645, 85903-000, Toledo, Parana, Brazil. * Author for correspondence. E-mail: eq.vilmar@hotmail.com

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.

Material and methods

Modeling

The mass and energy balances was undertaken during the development of the model. It was presumed that outlet streams in each stage ([L.sub.j] and [V.sub.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, [xi] 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.

[FIGURE 1 OMITTED]

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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (2)

in which:

[R.sub.j] = [[N.sub.r].summation over (t=1)][[N.sub.c].summation over (t=1)] [v.sub.li] x [[xi].sub.lj] (3)

where:

F is the feed stream, U is the liquid side stream, W is the vapor side stream, [r.sub.B] is the reboil ratio, [r.sub.D] is the reflux ratio, q is the feed quality stream, v is the stoichiometric coefficient, [N.sub.C] is the number of components, [N.sub.P] is the number of stages and [N.sub.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:

[V.sup.(0).sub.2] = ([r.sub.D] + 1) x [D.sup.(0)] - [F.sub.1] (4)

[V.sup.(0).sub.j] = [V.sup.(0).sub.j-1] + [W.sub.j-1] + ([q.sub.j-1] - 1) x [F.sub.j-1] (5)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (6)

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:

[L.sup.(0).sub.1] = [r.sub.D] x [D.sup.(0)] (7)

[L.sup.(0).sub.j] = [V.sup.(0).sub.j+1] = [V.sub.1] - D + [j.summation over (m=1)] [F.sub.m] + [j.summation over (m=2)] ([R.sub.m] - [U.sub.m] - [W.sub.m]), j = 2, 3, ..., [N.sub.p.sup.-1] (8)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

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.

[d.sub.1] - [X.sub.i1] + [u.sub.1] [x.sub.12] = [b.sub.1] (10)

[l.sub.j] - [X.sub.ij-1] + [d.sub.j] - [x.sub.ij] + [u.sub.j] - [x.sub.ij+1] = [b.sub.j], j = 2, 3, ..., [N.sub.p.sup.-1] (11)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

in which:

[l.sub.j] = [L.sub.j-1] j = 2, 3, ..., [N.sub.p] (13)

[d.sub.1] = - [[V.sub.1] x [K.sub.i1] + [L.sub.1] + D] (14)

[d.sub.j] = -[([U.sub.j] + [L.sub.j]) + ([W.sub.j] + [V.sub.j]) x [K.sub.ij]], j = 2, 3, ..., [N.sub.p.sup.-1] (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[u.sub.j] = [V.sub.j+1] x [K.sub.ij+1], j = 1, 2, ..., [N.sub.p.sup.-1] (17)

[b.sub.j] = - [F.sub.j] x [z.sub.ij] [[N.sub.r].summation over (l=1)] [v.sub.li] x [[xi].sub.lj], j = 1, 2, ..., [N.sub.p] (18)

The phase equilibrium relation is given as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

Equations (10)-(12) make up a tri-diagonal linear system with [N.sub.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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

Calculation of temperatures is undertaken by the re-arrangement of the phase equilibrium equation using the secant method to solve the equations:

F([T.sub.j]) = [[N.sub.c].summation over (i=1)][K.sub.ij] x [x.sub.ij] = 0, j = 1, 2, ... [N.sub.p] (21)

[T.sup.k+1.sub.j] = [T.sup.k.sub.j] - F([T.sup.k.sub.j]) x [[T.sup.k.sub.j] - [T.sup.k-1.sub.j]/ F([T.sup.k.sub.j]) - F([T.sup.k-1.sub.j]] (22)

Mole fractions of the vapor phase are calculated directly by the phase equilibrium equation:

[y.sub.ij] = [K.sub.ij] x [x.sub.ij], i = 1, 2, ..., [N.sub.c]; j = 1, 2, ..., [N.sub.p] (23)

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):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (24)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (25)

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:

[Q.sub.C] = [F.sub.1] x [h.sub.Feed 1] + [V.sub.2] x [H.sub.2] - [V.sub.1] x [H.sub.1] - ([L.sub.1] + D) x [h.sub.1] (26)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

The combination of molar balances, definition of reboil ratio and energy balances, with modifications, provide the following set of linear equations:

[[alpha].sub.j] x [V.sub.j] + [[beta].sub.j] x [V.sub.j+1] = [[gamma].sub.j], j = 2, 3, ..., [N.sub.p.sup.-1] (28)

in which:

[[alpha].sub.j] = [h.sub.j-1] - [H.sub.j] (29)

[[beta].sub.j] = [H.sub.j+1-] - [h.sub.j] (30)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (31)

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (32)

in which:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (33)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (34)

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.

D = [V.sub.2] + [F.sub.1]/[r.sub.D] + 1 (36)

Equations of the molar balances may be rearranged so that the liquid phase streams may be calculated from the following equations:

[L.sub.1] = [r.sub.D] x D (37)

[L.sub.j] = [V.sup.j+1] - [V.sub.1] - D + [j.summation over (m=1)] [F.sub.m] + [j.summation over (m=2)] ([R.sub.m] - [U.sub.m] - [W.sub.m]), j = 2, 3, ..., [N.sub.p.sup.-1] (38)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (39)

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.

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (40)

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 equilibrium-limited and occurs on stages 2 to 13. The following chemical equation describes the 2-pentene metathesis:

2[C.sub.5][H.sub.10]][] [C.sub.4][H.sub.8] + [C.sub.6][H.sub.12] (41)

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).

[FIGURE 2 OMITTED]

[FIGURE 3 OMITTED]

A code in FORTRAN was developed to solve the case. The program converged with 236 iterations, with 1 x [10.sup.-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 ?so-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:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (42)

The equation for the reaction rate was given by Venimadhavan et al. (1994):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (43)

Peng-Robinson equation of state was employed to represent the non-ideal behavior of the vapor phase, whereas Wilson model was used for the non-ideal 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.sup.-10] tolerance, using Intel Core2 QUAD Q6600 2,40GHz processor.

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.

[FIGURE 4 OMITTED]

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.

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 5 OMITTED]

Conclusion

An algorithm was provided in this work for the solution of a reactive distillation column that need solution for non-linear equations system.

Equations in the model were divided into sets and each set was solved separately. A methodology was also provided to produce initial guess which constitutes a critical step in the solution of nonlinear equations system.

The suggested methodology for the production of initial estimates was efficient with values close to those calculated by simulation. This fact accelerated and increased the convergence warranty.

The algorithm used was also efficient since the results obtained with the code developed were very close to those in the literature. There were no convergence problems in the solution of the equation system.

Nomenclature

B Bottoms;

D Distillate;

[D.sub.a] Damkohler number;

[F.sub.j] Feed rate of stage j;

h Enthalpy of liquid phase;

H Enthalpy of vapor phase;

[L.sub.j] Liquid rate of stage j;

[N.sub.c] Total number of components;

[N.sub.p] Total number of stages;

[q.sub.j] Feed quality of stage j;

[Q.sub.j] Heat load from stage j;

[r.sub.B] Reboil ratio;

[r.sub.D] Reflux ratio;

[T.sub.j] Temperature of stage j;

[U.sub.j] Liquid side rate of stage j;

v Volume;

[V.sub.j] Vapor rate of stage j;

[W.sub.j] Lateral vapor outlet of stage j;

[x.sub.ij] Liquid mole fraction of component i on stage j;

[[gamma].sub.ij] Vapor mole fraction of component i on stage j;

[z.sub.ij] Mole fraction of component i on stage j.

Doi: 10.4025/actascitechnol.v34i1.9535

References

ALFRADIQUE, M. F.; CASTIER, M. Modeling and simulation of reactive distillation columns using computer algebra. Computers and Chemical Engineering, v. 29, n. 4, p. 1875-1884, 2005.

BAUR, R.; HIGLER, A. P.; TAYLOR, R.; KRISHNA, R. Comparison of equilibrium stage and nonequilibrium stage models for reactive distillation. Chemical Engineering Journal, v. 76, n. 1, p. 33-47, 2000.

BROYDEN, C. G. A class of methods for solving nonlinear simultaneous equations. Mathematics of Computation, v. 19, n. 92, p. 577-593, 1965.

CHEN, F.; HUSS, R. S.; MALONE, M. F.; DOHERTY, M. F. Simulation of kinetic effects in reactive distillation. Computers and Chemical Engineering, v. 24, n. 11, p. 2457-2472, 2000.

CHENG, Y. C.; YU, C. C. Effects of feed tray locations to the design of distillation and its implication to control. Chemical Engineering Science, v. 60, n. 17, p. 4661-4677, 2005.

CORAZZA, F. C.; OLIVEIRA, J. V.; CORAZZA, M. L. Application of a subdivision algorithm for solving nonlinear algebraic systems. Acta Scientiarum. Technology, v. 30, n. 1, p. 27-38, 2008.

FERNHOLZ, G.; ENGELL, S.; KREUL, L.-U.; GORAK, A. Optimal operation of a semi-batch reactive distillation column. Computers and Chemical Engineering, v. 24, n. 2-7, p. 1569-1575, 2000.

GANGADWALA, J.; KIENLE, A. MINLP optimization of butyl acetate synthesis. Chemical Engineering and Processing, v. 46, n. 2, p. 107-118, 2007.

KHALEDI, R.; BISHNOI, P. R. A method for modeling two- and three-phase reactive distillation columns. Industrial and Engineering Chemical Research, v. 45, n. 17, p. 6007-6020, 2006.

KROUMOV, A. D.; MODENES, A. N.; WENZEL, B. M. Desenvolvimento de um modelo da cinetica enzimatica da transesterificacao de oleos vegetais para producao de biodiesel. Acta Scientiarum. Technology, v. 29, n. 1, p. 9-16, 2007.

MARQUINI, M. F.; MARIANI, D. C.; MEIRELLES, A. J. A.; SANTOS, O. A. A.; JORGE, L. M. M. Simulacao e analise de um sistema industrial de colunas de destilacao de etanol. Acta Scientiarum. Technology, v. 29, n. 1, p. 23-28, 2007.

POLING, B. E.; PRAUSNITZ, J. M.; O'CONNELL, J. P. The properties of gases and liquids. 5th ed. Digital Engineering Library: McGraw-Hill, 2004.

QI, Z.; KIENLE, A.; STEIN, E.; MOHL, K.-D.; TUCHLENSKI, A.; SUNDMACHER, K. MTBE decomposition in a reactive distillation column. Chemical Engineering Research and Design, v. 82, n. 2, p. 185-191, 2004.

SEADER, J. D.; HENLEY, E. J. Separations process principles. New York: John Wiley and Sons, 1998.

SINGH, B. P.; SINGH, R.; KUMAR, M. V. P.; KAISTHA, N. Steady-state analyses for reactive distillation control: An MTBE case study. Journal of Loss Prevention in the Process Industries, v. 18, n. 4-6, p. 283-292, 2005.

TAYLOR, R.; KRISHNA, R. Modelling reactive distillation. Chemical Engineering Science, v. 55, n. 22, p. 5183-5229, 2000.

VENIMADHAVAN, G.; BUZAD, G.; DOHERTY, M. F.; MALONE, M. F. Effect of kinetics on residue curve maps for reactive distillation. American Institute of Chemical Engineers Journal, v. 40, n. 11, p. 1814-1824, 1994.

UNG, S.; DOHERTY, M. F. Vapor-liquid phase equilibrium in system with multiple chemical reactions. Chemical Engineering Science, v. 50, n. 1, p. 23-48, 1995.

YAWS, C. L. Yaws' Handbook of thermodynamics and physical properties of chemical compounds. New York: Knovel, 2003.

Received on March 4, 2010.

Accepted on October 26, 2010.

License information: This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Vilmar Steffen * and Edson Antonio da Silva

Curso de Engenharia Quimica, Universidade Estadual do Oeste do Parana, R. da Faculdade, 645, 85903-000, Toledo, Parana, Brazil. * Author for correspondence. E-mail: eq.vilmar@hotmail.com

Table 1. Operational conditions of reactive distillation column in Case 1. Variables Specifications Pressure All stages 1 atm Ratio Reflux 4.0 Reboil 4.151 Feed Saturated liquid Stage 7 Flow rate 100 kmol [h.sup.-1] Temperature 310.1 K Pressure 1 atm Mole fraction 2-pentene (1.0) Table 2. Comparison between results obtained in this work and those by Chen et al. (2000) for Case 1. Variables Specifications Chen et This work al. (2000) Rate (kmol Distillate 50.0 50.0 [min..sup.-1]) Bottoms 50.0 50.0 Mole fraction (D) 2-pentene 0.0199 0.0220 2-butene 0.9801 0.9779 3-hexene 7.08x 6.30x [10.sup.-5] [10.sup.-5] Mole fraction (B) 2-pentene 0.0199 0.0220 2-butene 6.50x 6.29x [10.sup.-5] [10.sup.-5] 3-hexene 0.9801 0.9779 Table 3. Operational conditions of distillation column in Case 2. Variables Specifications Pressure All stages 11 atm Ratio Reflux 7 Reboil 10,44 * Feed 1 Liquid Stage 10 Rate 711.30 kmol [h.sup.-1] Temperature 320 K Pressure 11 atm Mole fractions MeOH (1.00) Feed 2 Vapor Stage 11 Rate 1965.18 kmol [h.sup.-1] Temperature 350 K Pressure 11 atm Mole fractions i-B (0.36) n-B (0.64) Molar holdup ** Reactive stages (j = 4-11) 68529.3 kmol * Adjusted for B = 640.8 kmol [h.sup.-1]; ** Da = 100, reaction equilibrium column (CHEN et al., 2000). Table 4. Comparison between results obtained in this work and those by Singh et al. (2005) for Case 2. Variables Specifications Singh et This work al. (2005) Rate (kmol Distillate 2035.6 1394.6 [h.sup.-1]) Bottoms 640.8 640.8 Mole fractions (D) i-B 0.0491 0.0477 MeOH 0.0518 0.0505 MTBE 0.0003 0.0000 n-B 0.8988 0.9018 Mole fractions (B) i-B 0.0001 0.0002 MeOH 0.0001 0.0000 MTBE 0.9986 0.9968 n-B 0.0012 0.0030