Investigation of Numerical Solution Approaches to Multicomponent

Feb 15, 1997 - Higher CPU time associated with higher order of the polynomial necessitates use of ... tributed in time and at least the axial directio...
0 downloads 0 Views 180KB Size
Download PDF
1738

Ind. Eng. Chem. Res. 1997, 36, 1738-1746

Investigation of Numerical Solution Approaches to Multicomponent Batch Distillation in Packed Beds R. M. Wajge, J. M. Wilson,† J. F. Pekny, and G. V. Reklaitis* School of Chemical Engineering, Purdue University, West Lafayette, Indiana 47906

Finite difference and orthogonal collocation techniques are used to obtain numerical solution of the differential contactor model that represents a packed column. In the orthogonal collocation approach, polynomial approximation to the model equations results in a sparse system of equations that is much smaller in dimension than with the other methods. Exploiting this sparsity and choosing the proper order for approximating polynomial are the critical issues. Higher CPU time associated with higher order of the polynomial necessitates use of finite elements in the collocation techniques. We study the accuracy and efficiency issues underlying all these approaches with the help of hydrocarbon mixture and ethanol esterification examples. The later case study is considered with both with and without the presence of reaction. 1. Introduction The ancient chemical engineering operation of batch distillation has regained the attention of chemical engineers in manufacturing environments, where demand does not justify a dedicated continuous separation facility (Macchietto and Mujtaba, 1995). Specifically, batch distillation in packed columns has been recognized as an effective method for carrying out separation of fine chemicals that is both more flexible and often cheaper than a staged column (Srivastava and Joseph, 1984). The dynamic nature of batch distillation introduces several key simulation issues which do not arise in the continuous case. First, dynamic models must be distributed in time and at least the axial direction (Aly et al., 1990a). Secondly, an efficient solution technique must be employed for the large systems of equations which typically arise. Conventional simulation packages are unable to solve the differential contactor model equations which represent batch distillation simultaneously for all variables over the time span of interest without using an equation-tearing algorithm to mitigate numerical difficulties (Hitch and Rousseau, 1988). This results in an iterative approach that requires the repetitive solution of the system of model equations until convergence is obtained. Efficiency in the execution of models also becomes critical when such models are used to determine optimal operation policies in real time. While in the domain of staged models, short-cut methods are fairly established for reducing computational time (Diwekar and Madhavan, 1991; Sundaram and Evans, 1993), the corresponding continuous models have received little attention in the literature. This work addresses key numerical issues involved in the development of an efficient simulation model and solution method for packed batch columns. 2. Packed Column Simulation 2.1. Model Formulation and Past Work. The models used to represent packed column distillations differ from their staged counterparts in that the relationships between vapor and liquid stream compositions are no longer solely determined by equilibrium relations * Author to whom correspondence is to be addressed. † Present address: Monsanto Crop Protection Division, Muscatine, IA 52761. S0888-5885(96)00296-5 CCC: $14.00

Figure 1. Differential contactor model. Table 1. Model Variables number

description

symbol

NC NC NC NC 1 1 1

liquid composition gas composition liquid interfacial composition gas interfacial composition liquid flow rate gas flow rate interfacial temperature

xj yj x*j y*j L′ G′ T*

but also by mass transfer effects. One such model is the differential contactor model which is depicted in Figure 1 and in eqs 1-7. The model variables are given in the Table 1 (Hitch and Rousseau, 1988). Equation 1 represents the unsteady state individual component material balances.

∂bLj ∂bGj ∂(L′xj) ∂(G′yj) + ) ∂t ∂t ∂z ∂z

(1)

The mass transfer relationships for both the liquid and vapor phases for each component are given by eqs 2 and 3, respectively. © 1997 American Chemical Society

Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 1739

∂bLj ∂(L′xj) ) + k′xja(x* j - xj) ∂t ∂z

(2)

∂bGj ∂(G′yj) )- k′yja(yj - y* j) ∂t ∂z

(3)

The relationship between the interfacial vapor and liquid mole fractions is given by eq 4, while the compositions in both phases sum to 1.0 as given by eq 5.

y* j ) mjx* j NC

(4)

NC

y* ∑ j ) ∑mjx* j ) 1.0 j)1 j)1

(5)

Equations 6 and 7 represent the unsteady state adiabatic energy balance and the unsteady state total material balance, respectively.

∂(BLHL) ∂(L′HL) ∂(G′HG) ) ∂t ∂z ∂z

(6)

∂BL ∂BG ∂L′ ∂G′ + ) ∂t ∂t ∂z ∂z

(7)

The bulk liquid- and gas-phase temperature are assumed to be same as the interfacial temperature. The details of the model equations and the underlying assumptions can be found in Hitch and Rousseau (1988). By integrating these equations, the behavior of the packed column can be predicted. Unfortunately, the exact solution of the differential contactor model is difficult to obtain due to the presence of algebraic variables and equations. One numerical method that has received attention in the literature is the use of finite difference techniques to convert the derivative expressions into algebraic expressions. The resulting algebraic system can be readily solved using a linear or nonlinear equation-solving package. Von Rosenberg and Hadi (1980) used the finite difference method on a simplified version of the differential contactor model to determine the steady state behavior of a packed column. The application they considered involved the separation of a mixture of toluene, p-xylene, m-xylene, and methylcyclohexane. The results obtained indicated that the predicted system behavior using 20 points agreed with that predicted using 50 points over a 20 ft column. No results were reported for the actual computational times required to obtain the steady state results. The finite difference method has also been applied by Hitch and Rousseau (1988) to the differential contactor model. The system they chose to investigate contained propane, n-butane, and n-hexane. In contrast to Von Rosenberg and Hadi (1980), Hitch and Rousseau used the simulation to study the dynamic column behavior. In their formulation, they incorporated detailed physical property calculations rather than treating the physical properties as being constant throughout the column. By varying boil-up rate, column height, and reflux ratio, Hitch and Rousseau found that the differential contactor model gave results that matched qualitative expectations. However, the use of finite difference methods requires the solution of large systems of equations. Even though these systems are highly structured and that structure can be exploited by block tridiagonal

solution method, a formulation that requires fewer grid points can yield significant computational savings. This suggests the use of potentially more efficient methods of solving the differential contactor model, viz., the method of orthogonal collocation. Srivastava and Joseph (1984) applied orthogonal collocation to a simplified form of the differential contactor model to determine the steady state and dynamic behavior of a packed distillation column. Using the simplifying assumptions of constant liquid and vapor flow rates, linear equilibrium relationships, and saturated liquid and vapor streams, they modeled the steady state behavior of the system consisting of water and ethanol in a 1.2 m column. Their work showed that a fourth-order polynomial could give the same accuracy as a finite difference technique with 25 grid points. However, as with the finite difference method, orthogonal collocation can become computationally intensive as the polynomial order is increased. For example, when two collocation points were added to increase the polynomial to sixth order, the solution time nearly doubled. One way of keeping the polynomial order low is to use finite elements combined with some form of collocation. Recently, Aly et al. (1990a,b) used the finite elements method together with the Galerkin method to solve the differential contactor model. The Galerkin method is more difficult to implement than the orthogonal collocation techniques as it involves solution of integral equations as opposed to solution of algebraic equations. The specific example system studied consisted of cyclohexane and toluene, and the results indicated that the finite element method gave the same solution as the finite difference approach. The finite difference method was found to be slower in convergence than the finite element method. In addition, they found that in going from four to six elements, the time to determine the solution doubled, while in going from four to eight, the solution times tripled. 2.2. Reactive Distillation. The incorporation of reactions into batch distillation processes is receiving increasing attention in the simulation literature. Cuille and Reklaitis (1986) discussed reactive batch distillation taking place in tray or staged columns. Reuter et al. (1989) studied simulation and control aspects of reactive distillation processes. Doherty and Buzad (1992) presented a survey of the available design techniques for reactive distillation and pointed out distinct advantages of this unit operation. The most important of these advantages is that reactive distillation can be used to overcome yield limitations imposed by chemical equilibrium. With the presence of a chemical reaction, the unsteady state individual component material balance for the differential contactor model, eq 1 transforms into eq 8. The inclusion of a reaction term (which is a

∂bLj ∂bGj ∂(L′xj) ∂(G′yj) + ) + Rj ∂t ∂t ∂z ∂z

(8)

function of composition of the species and temperature) adds to the numerical difficulties presented by the nonlinearities and the index of the differential algebraic equations. 2.3. Model Implementation. The system of eqs 1-7 constitutes 4NC + 3 equations. The independent model variables are given in Table 1. The underlying goal of the finite difference or collocation methods is to

1740 Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997

convert the 4NC + 3 equations into a system of ordinary differential equations which can be integrated simultaneously using a readily available ODE solver such as LSODE (Hindmarsh, 1980). On the basis of an inspection of the model equations, three numerical problems can be identified: (i) the differential equations contain algebraic variables; (ii) the differential equations contain nonlinear terms; (iii) eqs 4 and 5 represent purely nonlinear algebraic equations. Specifically, the differential contactor model equations yield three types of algebraic variables, T, x* j, and y* j. In general, systems of differential algebraic equations cannot be treated as ordinary differential equations due to problems associated with index (Unger et al., 1995). Finding the index of the system of model equations is highly complex due to the presence of numerous correlations for the estimation of physical properties and holdup. However, it can be assumed that the system will have index greater than one, in general. Such systems are known to present numerical difficulties in terms of initialization and convergence of the solution algorithm (Unger et al., 1995). In addition, the nonlinearity of eqs 4 and 5 necessitates an iterative approach. Therefore, an equation-tearing or -pairing algorithm must be imcorporated that not only addresses the problems posed by index but also handles the nonlinear terms in the model equations. To linearize the model equations, the dependent variable functionality must be eliminated. This can be accomplished by treating the values of dependent variables such as surface tension and viscosity as constants based on previous values for the independent variables. This paired approach was proposed by Hitch and Rousseau (1988). It essentially requires that the solution for a specific variable be determined by the equation that has the greatest impact on the value of that variable. Since distillation is normally performed on a set of components with a narrow boiling point spread, changes in enthalpy will represent changes in total flow rates but will have little effect on interfacial temperatures. The interfacial temperatures are best determined by the composition profiles throughout the column. This implies that in eqs 1-4 the values for L′, G′, BL, and BG are treated as constants determined from the model solution at the previous time step. Likewise, for eqs 6 and 7, the values for HL, HG, BL, and BG are treated as constants and an iterative approach is required to ensure that the solution converges (column holdup may play an important role in the rate of convergence, in general, and a higher column holdup may necessitate use of smaller step size). Once the model equations have been linearized, they take on the following form: Equation 1 or the unsteady state individual component balance becomes eq 9.

∂bLj ∂bGj L′ ∂xj G′ ∂yj + ) ∂t ∂t ∂z ∂z

(9)

The liquid and vapor phase mass transfer relationships, eqs 2 and 3, become eqs 10 and 11, respectively.

∂bLj L′ ∂xj ) + k′xja(x* j - xj) ∂t ∂z

(10)

∂bGj G′ ∂yj )- k′yja(yj - y* j ) ∂t ∂z

(11)

The interfacial vapor and liquid relationship, eq 4, along with the summation relationship, eq 5, remain unchanged during the linearization. However, the unsteady state adiabatic energy balance and material balance, eqs 6 and 7, become eqs 12 and 13.

BL ∂HL HL ∂L′ HG ∂G′ ) ∂t ∂z ∂z 0)

∂L′ ∂G′ ∂z ∂z

(12) (13)

Using the equation-tearing method to bypass the problems associated with index as well as the nonlinear terms, eqs 9-11 as well as eq 4 are solved simultaneously to yield the new bulk and interfacial compositions xj, x*j, yj, and y* j. The new interfacial compositions are then used to determine the interfacial temperature, T, which is in turn used to determine new gas and liquid enthalpies, HG and HL. Finally, the gas and liquid enthalpies are used in conjunction with eqs 12 and 13 to determine the liquid and gas flow rates throughout the column, L′ and G′. These values, as well as all constants based on the newly calculated compositions, are then substituted back into the model eqs 9-11 and 4, and the process is repeated until the difference between consecutive solutions has met the convergence criterion. The complete details of this algorithm can be found in Hitch and Rousseau (1988). Though the algorithm cannot be applied directly to azeotropic systems, the extended use of coordinate transformation (Anderson and Doherty, 1984; Diwekar, 1991) could be considered. 3. Simulation Results and Discussion 3.1. Finite Difference Method. Using the above equation-tearing algorithm, the finite difference technique was employed to convert the model equations into a system of linear algebraic equations. We first considered the light hydrocarbon system (propane, butane, and hexane) studied by Hitch and Rousseau (1988). Although our results match the results of Hitch and Rousseau (1988) qualitatively (Figure 2) in the sense that the components are removed in order of decreasing volatility, the actual composition profiles do not agree since Hitch and Rousseau obtained much lower purities. The maximum propane and butane concentrations they obtained were 82% and 70%, respectively, as opposed to our values of 99% and 97%. Such high concentrations are also reported by Mujtaba and Macchietto (1988) with the system of propane, butane, pentane, and hexane. They observed that the systems composed of low molecular weight hydrocarbons behave as a series of binary batch distillations. Our results (Figure 2) obtained using a grid of size 0.1 (for both time and height) also demonstrate the same general sequential binary behavior as that observed by Mujtaba and Macchietto. We assume that the differences in our composition profiles and those generated by Hitch and Rousseau are attributable to different physical properties data. Because of the unavailability of such details in the literature, it is difficult to find out the exact physical property routine(s) causing the quantitative differences, and without actual experimental data to evaluate the results from the differential contactor

Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 1741

Figure 2. Comparison of finite difference results (our results are indicated by solid lines, while points present the results of Hitch and Rousseau).

Figure 3. Finite difference coefficient matrix.

model, a further assessment of the results is not possible (see Appendix for further comments). The finite difference solution approach is computationally time consuming. The results shown in Figure 2 required a Sparc 2 workstation to run for a period of nearly 48 h. This results in nearly 3.5 h of CPU time for every hour of simulated time. If the ultimate goal of simulating packed column distillation is to incorporate the simulation within a control or optimization framework, then the simulation must be much more efficient. The required speedup in solution times can be obtained either by reducing the number of equations to be solved or by solving the equations more efficiently. In the finite difference approach, the set of equations take on the block tridiagonal structure shown in Figure 3. Each block of the coefficient matrix is of size 4 (NC), where the total number of blocks is given by 3(NGP 2) + 4. Even though efficient block tridiagonal solvers exist, the computational overhead required to evaluate all physical and chemical properties at each grid point becomes excessive as the number of grid points is increased. Therefore, a numerical technique that can reduce the number of equations to be solved through a reduction in the number of grid points within the column should yield computational savings. One such method is orthogonal collocation. 3.2. Orthogonal Collocation Method. 3.2.1. Global Collocation. The goal of collocation or polynomial approximation techniques is to reduce the size of the systems to be solved while still yielding highly accurate

solutions. The difference between the approximation and the actual solution, or residual, decreases with increasing polynomial order. However, by increasing the order of the polynomial, the amount of work required to solve the system of equations increases on the order of n3 when dense matrix techniques are used. Hence, the choice of polynomial order depends on the task and accuracy required. Specifically, as long as the function F(z) to be approximated by a nth degree polynomial is a polynomial of degree 2n + 1 or less, the approximation will be exact at the basis points (Carnahan et al., 1969). In those cases where F(z) is of degree greater than 2n + 1, even though the residual is forced to vanish at the basis points, the approximation need not be exact. With respect to distillation simulation, when the primary interest is in obtaining overhead compositions rather than detailed composition profiles throughout the column, the lowest order polynomial subject to allowable approximation error should be selected. However, if the goal is to determine accurate composition profiles, the polynomial should be chosen to yield as low a residual as possible, and its order n should be such that the function degree is less than or equal to 2n + 1. With respect to the size of the system of equations to be solved per iteration, global collocation required the solution of a 96 × 96 system (using a seventh-order Legendre polynomial for the same above hydrocarbon example) that is 7.9% dense, while the finite difference technique results in the block tridiagonal structure with blocks that are 8.3% dense. The size of the system of equations to be solved is described by 4(NC)(NCP). If the number of collocation points (NCP) is less than the number of grid points used in the finite difference method (NGP), then the computational overhead associated with the calculation of chemical and physical properties will be reduced. In addition, by eliminating the holdup within the column or assuming pseudo steady state, the similarity between no-holdup and lowholdup systems allows for a further reduction in work required per simulation. By eliminating column holdup with the assumption that the dominant holdup is located in the reboiler, the model equations are no longer integrated through time; rather, they are solved at a point in time depending on the integration of the reboiler balance equations (due to the temperature dependency of the reaction term in eq 8, paired approach cannot be applied unless the column holdup is neglected). This approach is readily implemented using the integration package LSODE (Hindmarsh, 1980). The flow chart of this simulation (henceforth referred to as FE simulation) can be seen in Figure 4. To study the effect of order of the polynomial on the composition profiles and the performance of the FE simulation, we chose the system consisting of ethyl acetate (EtOAc), ethyl alcohol (EtOH), acetic acid (AcOH), and water. This is a reactive system, but we assume no chemical reaction among the components for the initial analysis. The input parameters and computational results for collocation are given in Tables 2 and 3, respectively. The reflux profiles used are shown in Figure 5. Unlike the light hydrocarbon example studied earlier, this system has a very narrow boiling point spread between the two components with the highest volatility. As a result, the system cannot be expected to behave as a series of binary separations. These expectations are confirmed by the composition profiles generated by run 1 as shown in Figure 6. The results

1742 Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997

Figure 4. Finite element simulation flow chart, FE. Table 2. Simulation Parameters for Collocation Runs run

reflux

boil-up (kmol/h)

initial feed (kmol)

1-3, 7-9 10-12 4-6

A

150

10-EtOAc, 10-EtOH, 10-AcOH, 10-water

B

150

10-EtOAc, 10-EtOH, 10-AcOH, 10-water

Table 3. Computational Results for Collocation Method run

height (ft)

element

order

sim. (h)

CPU (h)

1 2 3 4 5 6 7 7s 8 9 9s 10

10 10 10 10 10 10 10 10 10 10 10 20

1 1 4 1 1 4 1 1 2 4 4 4

3 4 3 3 4 3 3 3 3 3 3 3

3:50 3:50 3:50 2:50 2:50 2:50 3:00 3:00 3:00 3:00 3:00 3:00

1:04 1:46 10:43 0:33 0:53 8:53 1:04 0:16 3:43 7:34 6:30 12:19

using global collocation (runs 2, 4, 5, 7, and 7s) show that collocation techniques using low-order polynomials can produce solutions in reasonable CPU times. However, as the size of the system increases (either due to an increase in the height of the distillation column or in the order of the approximating polynomial), the reductions in CPU times obtained via global collocation can be lost. Therefore, we must use a method that employs the low-order polynomials while yielding the desired accuracy. This can be accomplished by using collocation on finite elements. 3.2.2. Collocation on Finite Elements. The application of the finite element approach results in a linear system as outlined in Figure 7. The size of the

Figure 5. Reflux profiles A (left) and B (right).

individual blocks along the diagonal is controlled by the same relationship that was used for the global collocation, viz., 4(NC)(NCP). This system is mildly block tridiagonal and was solved using the same block tridiagonal solver developed for the finite difference approach whenever three or more elements were used. To study accuracy and efficiency issues between global collocation and collocation over finite elements, several runs were taken with the same simulation input parameters as used in the global collocation runs (Table 2). The effect of order of polynomial on the composition profiles (runs 1, 2, 4, and 5) can be seen from the Figures 8 and 9. The composition profiles generated by collocation over finite elements (runs 3 and 6) are shown in these figures as a measure of accuracy. We observe that the composition profiles generated by third-order poly-

Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 1743

Figure 6. Global collocation results (run 1).

Figure 7. Finite element coefficient matrix.

Figure 9. Collocation results: runs 4-6.

Figure 8. Collocation results: runs 1-3.

nomial match with those generated by a fourth-order polynomial or by collocation over finite elements within 5%. The same is true for the results of collocation over finite elements methods with a lower number of elements (runs 7, 8, and 9). Hence, if 5% deviation in the composition profiles is acceptable, then for the system under consideration the global collocation method with a third-order polynomial is adequate.

With respect to the computational results given for the finite difference simulation as well as the collocation techniques, the results are directly influenced by the use of dense matrix or sparse matrix solution techniques. Consider the use of dense matrix techniques. The system of equations that is developed through the use of finite difference techniques is block tridiagonal, where the individual block size is 4(NC). If the behavior of a 10 ft column is simulated using the finite difference method with a grid spacing of 0.1 feet, 100 grid points will be required. For a four-component system, the block size and overall system size will be 16 × 16 and 16002, respectively. The solution of the block tridiagonal system can then be determined in the order of 100(163) operations. For the same column modeled by a thirdorder Legendre polynomial approximation, the block size and system size will be 64 × 64, the solution of which will require order 643 operations. For a column simulated by two third-order polynomials, the solution time will be of the order of 2(643) elementary mathematical operations. These estimates suggest that the solution using the global polynomial approach will take roughly half the operations that the finite difference

1744 Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 Table 5. Input Parameters for Runs 1R, 2R, and 3R

Figure 10. Increase in nonzero matrix elements with incresaing elements. Table 4. Increase in Nonzero Elements with Order num. comp.

order

NZ elements

3 3 3

3 4 7

177 279 729

method requires, while the operations required for the finite difference method and finite element collocation approach are roughly identical. However, these comparisons do not take into account the computational overhead associated with the increased number of grid points used by the finite difference method when compared to collocation methods. While the finite difference approximation will require physical and chemical properties to be calculated at 100 grid points, the collocation methods will require only 4 or 8 points. Hence, the number of calculations is reduced by an order of magnitude for the collocation techniques compared to the finite difference methods. When using Gaussian elimination routines or block tridiagonal routines based on dense Gaussian elimination, the increase in nonzero terms has little effect on the solution time; rather, it is the overall system size, or block size, that dominates. However, sparse matrix routines that manipulate only nonzero elements can be used to improve solution times as evident from run 7s. The solution time of run 7s (using sparse matrix method) is less than that of run 7 by a factor of 4. Specifically, the addition of grid points to the finite difference approach or collocation over finite elements approach results in a linear increase in the nonzero elements within the coefficient matrix. For a threecomponent system, the increase in nonzero elements associated with an increase in the number of third- and fourth-order polynomial segments (or finite elements) within the column is linear as given in Figure 10, whereas the increase in nonzero coefficient matrix elements with increasing polynomial order is nonlinear (Table 4). As was noted earlier, the individual blocks in the block tridiagonal structure developed by the use of the finite difference approach were roughly 8% dense. Hence, the sparse matrix routines will be forced to manipulate on the order of 3600 elements for a 10 ft column with 3 components. From Figure 10, this is

parameter

value

height elements order boil-up (kmol/h) initial feed (kmol)

10 1 3 150 0-EtOAc, 30-EtOH, 30-AcOH, 0-water

roughly an order of magnitude higher than the number of elements that a two-segment collocation simulation will use. Finally, with respect to the computational results shown in Table 3, the CPU overhead associated with physical and chemical property convergence necessitated by the equation-tearing algorithm and discussed earlier is evident. Specifically, the coefficient matrices used in both runs 9 and 10 are identical with respect to nonzero element location. The difference lies in the fact that run 9 repreents a 10 ft column while run 10 represents a 20 ft column. The CPU time required for run 10 on a sparc 2 is nearly twice that required for run 9. This is solely due to an increase in the number of iterations for the convergence of the chemical and physical properties within the column. Hence, in simulating longer column height, sparse matrix techniques must be employed to reduce computation times. 4. FE Simulation Applied to Reactive Distillation The previous discussion validated the FE approach for packed columns. In those simulations involved with the testing of both global collocation and collocation on finite elements, no reaction was assumed to take place. In this section, we apply the FE approach to the same chemical system but with the presence of the esterification reaction as given in eq 14. The sparsity of the coefficient matrix and the implementation of the algorithm remain the same for the reactive distillation cases; however, the system stiffness maybe higher depending on the kinetics. This would impose the need for smaller step size to achieve the desired accuracy, thus increasing the computational time. For the system under consideration, since the rate of reactions are relatively low, system stiffness does not change significantly. Since the assumption was made earlier that the reboiler holdup dominates the holdup within the column, the reaction is considered to take place only in the reboiler and not within the packing.

EtOH + AcOH S EtOAc + water

(14)

Holland (1981) gives the forward and reverse rate constants for this reaction, which are incorporated into the reactor balance equations.

L ‚‚‚ (-7150 T ) (g mol)(min) L -7150 k ) 7380 exp( ‚‚‚ T ) (g mol)(min)

kf ) 29 000 exp

r

(15) (16)

Due to the close boiling points of product EtOAc and reactant EtOH (350.2 and 351.7 K, respectively), the reflux ratio will have a significant effect on the purity of the product. Three simulation runs were taken to observe this effect. Input parameters for all these runs are listed in Table 5. Number of moles distilled and CPU time on a Sparc 2 are given in Table 6. The reflux profiles used and the compositions changes are shown

Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997 1745 Table 6. Results for Reactive Distillation Runs run

kmol distilled

CPU (h)

1R 2R 3R

58.37 56.70 49.04

7:39 5:42 3:54

Figure 11. Run 1R.

Figure 12. Run 2R.

Figure 13. Run 3R.

in Figures 11-13. Run 1R uses constant reflux, whereas runs 2R and 3R use total reflux (which allows the system to attain the chemical equilibrium at which the product purity is highest) and product withdrawal alternately. The reflux profile shown is in terms of the

percentage of overhead vapor flow rate recycled to the column (internal reflux ratio). We observe that total reflux operation results in an ethyl acetate concentration of up to 51%, which agrees well with the experimental and simulation results reported by Savkovic-Stevanovic et al. (1992). Macchietto and Mujtaba (1995), however, give high value of purity (75%). We again assume that the differences in the results are mainly due to the difference in physical properties data. Table 6 shows that the CPU time required for simulation increases as number of moles processed increases. This is due to the fact that toward the end of the distillation composition changes are rapid, thus taking more time for the convergence. During the total reflux stage (no product collection), the system reaches chemical equilibrium, thereby showing high overhead product purity. The purity drops as soon as product withdrawal starts. Thus, the chemical equilibrium becomes an important factor (in addition to the phase equilibrium) in designing reflux policy for reactive distillation. Hence, the reflux policies used in runs 2R and 3R, which allow the system to reach chemical equilibrium periodically, yield the product of higher purity compared to the constant reflux policy (run 1R). However, for obtaining the best trade-off between the purity of the product and the production rate, rigorous optimization study is required. The objective can be to maximize amount of (high-purity) product in the product cut and to maximize purity of product in the off-cut (so that its reprocessing is easy). The duration of these cuts should be such that the reactant level does not drop beyond a certain level at which the reverse reaction becomes dominant. The offcuts produced could be used to recycle the reactants provided that a suitable reprocessing scheme is established. This leads to the consideration of the optimization of the operation of an entire campaign of batches (Macchietto and Mujtaba, 1995) rather than focusing on the processing of a single batch (Wajge and Reklaitis, 1995). 5. Conclusions The simulation algorithm based on differential contactor model and orthogonal collocation techniques simulates multicomponent (nonreactive and reactive) batch distillation satisfactorily. While the collocation techniques are computationally more efficient compared to the finite difference method, the order of the approximating polynomial needs to be chosen judiciously to obtain a good trade-off between accuracy and efficiency. In cases where global collocation requires a very high order of polynomial for the accuracy in the results, collocation over the finite elements in conjunction with a low-order polynomial proves to be more efficient. The simulation results indicate that as the size of the system of equations increases or as composition profiles become widely separated within the column, the computational time increases substantially. Here, the use of sparse matrix techniques to solve the sparse block diagonal structure of the coefficient matrix offers significant reduction in CPU time. This efficiency can well be utilized in designing on-line optimal operation policies. The study on reactive batch distillation emphasizes the need for rigorous optimization study for designing optimal reflux policy and optimal campaign structure.

1746 Ind. Eng. Chem. Res., Vol. 36, No. 5, 1997

Nomenclature a ) interfacial area bGj ) vapor holdup of component j bLj ) liquid holdup of component j BG ) gas holdup BL ) liquid holdup  ) energy transfer rate G′ ) superficial vapor flow rate Hl ) liquid enthalpy Hv ) gas enthalpy kf ) forward reaction rate constant kr ) reverse reaction rate constant k′xi ) liquid phase mass transfer coefficient for component i k′yi ) vapor phase mass transfer coefficient for component i lij ) component liquid flow rate L ) liquid flow rate L′ ) superficial liquid flow rate mi ) phase equilibrium constant for component i n ) number of cuts N ) mass transfer rate NC ) number of components NCP ) number of collocation points NGP ) number of grid points Ri ) rate of formation of component i t ) time T ) temperature vij ) component vapor flow rate V ) vapor flow rate xi ) liquid phase mole fraction of component i x* ) liquid phase interfacial mole fraction yi ) vapor phase mole fraction y* ) vapor phase interfacial mole fraction z ) height from the bottom of the column Abbreviations AcOH ) acetic acid EtOH ) ethanol EtOAc ) ethyl acetate

Appendix Due to the unavailability of the experimental results in the literature, it is not possible to verify the quantitative accuracy of the composition profiles obtained in this and other studies. Such experimental results could also be helpful in finding the exact physical properties routines that cause these discrepancies. The details on the physical properties data used in the different studies can be found in Wilson (1993), Mujtaba (1989), and Hitch (1986). Literature Cited Aly, S.; Pibouleau, L.; Domenech, S. Treatment of Batch, Packed Distillation by a Finite Element Method. I. Steady State Model with Axial Dispersion. Int. Chem. Eng. 1990a, 30 (3), 452463. Aly, S.; Pibouleau, L.; Domenech, S. Treatment of Batch, Packed Distillation by a Finite Element Method. II. Steady State Model with Axial and Radial Dispersion and a Dynamic Model. Int. Chem. Eng. 1990b, 30 (3), 464-478.

Anderson, N. J.; Doherty, M. F. An Approximate Model for Binary Azeotropic Distillation Design. Chem. Eng. Sci. 1984, 39 (11), 11-18. Carnahan, B.; Luther, H. A.; Wilkes, J. A. Applied Numerical Methods; John Wiley and Sons: New York, 1969. Cuille, P. E.; Reklaitis, G. V. Dynamic Simulation of Multicomponent Batch Rectification with Chemical Reactions. Comput. Chem. Eng. 1986, 10 (4), 389-398. Diwekar, U. M. An Efficient Design Method for Binary, Azeotropic, Batch Distillation Columns. AIChE J. 1991, 37 (10), 15711578. Diwekar, U. M.; Madhavan, K. P. Multicomponent Batch Distillation Column Design. Ind. Eng. Chem. Res. 1991, 30, 713721. Doherty, M. F.; Buzad, G. Reactive Distillation by Design. Chem. Eng. Res. Des. 1992, 70A, 448-458. Hindmarsh, A. C. LSODE and LSODI: Two New Initial Value Ordinary Differential Equation Solvers. ACM Signum Newsletter 15 (4), 1980. Hitch, D. M. Numerical Simulation of Continuous-Contact Separation Processes. Ph.D. Dissertation, North Carolina State University, Raleigh, NC, 1986. Hitch, D. M.; Rousseau, R. W. Simulation of Continuous Contact Separation Processes Multicomponent Batch Distillation. Ind. Eng. Chem. Res. 1988, 27, 1466-1473. Holland, C. D. Fundamentals of Multicomponent Distillation; McGraw-Hill Book Co.: New York, 1981. Macchietto, S.; Mujtaba, I. M. Design of Operation Policies for Batch Distillation. In Batch Processing Systems Engineering: Current Status and Future Directions; Reklaitis, G. V., Ed.; Springer-Verlag: Berlin, 1995. Mujtaba, I. M. Optimal Operational Policies in Batch Distillation. Ph.D. Dissertation, Imperial College, London, 1989. Mujtaba, I. M.; Macchietto, S. Optimal Control of Batch Distillation. Presented at the IMACS World Conference, Paris, 1988. Reuter, E.; Wozny, G.; Jeromin, L. Modeling of Multicomponent Batch Distillation Processes with Chemical Reaction and their Control Systems. Comput. Chem. Eng. 1989, 13 (4/5), 499570. Savkovic-Stevanovic, J.; Misic-Vukovic, M.; Bonciccaricic, G.; Trisovic, B.; Jezdic, S. Reactive Distillation with Ion Exchangers. Sep. Sci. Technol. 1992, 27 (5), 613-630. Srivastava, R. K.; Joseph, B. Simulation of Packed Bed Separation Processes using Orthogonal Collocation. Comput. Chem. Eng. 1984, 8 (1), 43-50. Sundaram, S.; Evans, L. B. Multicomponent Batch Distillation Column Design. Ind. Eng. Chem. Res. 1993, 32, 511-518. Unger, J.; Kroner, A.; Marquardt, W. Structural Analysis of Differential-Algebraic Equation SystemssTheory and Applications. Comput. Chem. Eng. 1995, 19 (8), 867-882. VonRosenberg, D. U.; Hadi, M. S. Numerical Simulation of Multicomponent Packed Tower Distillation Problems. Chem. Eng. Commun. 1980, 4, 313. Wajge, R. M.; Reklaitis, G. V. Campaign Optimization of Multicomponent Reactive Batch Distillation. Presented at the AIChE annual conference, Miami, FL, 1995. Wilson, J. M. The Numerical Simulation and Operating Policy Optimization of Reactive Packed Column Batch Distillation. M.S. Dissertation, Purdue University, West Lafayette, IN, 1993.

Received for review May 28, 1996 Revised manuscript received January 3, 1997 Accepted January 4, 1997X IE960296A

X Abstract published in Advance ACS Abstracts, February 15, 1997.