Ind. Eng. Chem. Res. 1990, 29, 1364-1382
1364
Multicomponent Three-phase Azeotropic Distillation. 2. Phase-Stability and Phase-Splitting Algorithms Brett P. Cairns and Ian A. Furzer* Department of Chemical Engineering, University of Sydney, Sydney, New South Wales, Australia 2006
A literature review of three-phase distillation calculation procedures highlights the large number of efforts that have been made to account for the appearance of a second liquid phase on a distillation tray. The numerical problems due to the introduction of a second liquid phase during an iteration scheme are summarized in a table containing the known 15 solution sets. Many of these solutions contain a liquid-phase-splitting algorithm. A new algorithm based on Michelsen’s phase-stability test and phase-splitting algorithm has been developed. New procedures in the iteration method and to identify liquid-phase instability have been produced and verified. These new procedures use analytical partial derivatives which can be determined from modern thermodynamic models to describe nonideal liquid-phase behavior. The predicted separation performance of a three-phase azeotropic distillation column using this new algorithm is in good agreement with the experimental data collected in part 1 of this series. Literature Review Extension of existing mathematical models describing standard two-phase distillation to cover three-phase distillation is relatively straightforward. Difficulties arise, however, when attempting to find a solution to the equation set. A standard three-phase stage, without side streams, is shown in Figure 1. At steady state, the following MESH equations can be written for stage n: component balance Vn+lYn+l,t (Vn.Yn,i
+ Ll,n-lXl,n-l,i
+ L2,n-1X2,n-l,1 + FnXF,n,r i = 1, 2, **., nc (1) =0
+ L l , n X l , n , i + &,nX2,n,i)
enthalpy balance vn+,Hn+l
Ll,n-lhl,n-l
+ &,n-lh,n-1
-I- FnhF,n
+
Qn
-
+ L I , n h l , n + ~ 5 2 , n h 2 , n )= 0 ( 2 )
(VnHn+,
equilibrium equations Yn,i
i = 1, 2,
= Kl,n,iXl,n,r
Yn,r = K 2 , n , i ~ 2 , n , i XZ,n,r
= K*n,rXl,n,r
nc
(3)
i = 1, 2 , ..., n,
(4)
i = 1, 2, * * * > nc
(5)
**.I
It is important to notice that only two of the equilibrium equations are independent, as division of eq 3 by eq 4 yields K*n,r = Kl,n,r/Kz,n,i (6) summation equations (7)
The total number of independent equations for a column consisting of N stages all within the three-phase region is then Nnc N 2Nn, 3N
N(3n, + 4)
(component balance) (enthalpy balance) (equilibrium equations) (summation equations)
(10)
0888-5885/90/2629-1364$02.50/0
It can be shown that, for a column operating with NP stages within the three-phase region, the number of equations describing the process is given by N(2nc + 3) + NP(nc+ 1)
which equals N(3n,
(11)
+ 4) when NP equals N.
Review of Previous Simulation Methods Method of Successive Flashes. To the best of our knowledge, the first three-phase distillation model was developed by Bril’ et al. (197413, 1975) as part of a systematic investigation of the heteroazeotropic distillation process (Bril’ et al., 1973a-d, 1974a,b, 1975, 1976, 1977). Bril’ and co-workersbegan their investigation by studying the predictions of activity coefficient models for VLE, VLLE (Bril’ et al., 1973a,b), and ternary LLE (Bril’ et al., 1973~).An algorithm was developed to solve the LLE calculations (Bril’ et al., 1973d), which was later expanded to cover multicomponent, multiphase equilibrium (Bril’ et al., 1974a). The VLLE algorithm was then used to perform a limited simulation of a total reflux three-phase column (Bril’ et al., 1974b). Bril’ et al. (1975) expanded this further for a column operating at partial reflux by modifying the relaxation method described by Hanson et ai. (1962) to include three phases. This method used a flash equation that lead to the naming of the technique as the method of successive flashes. Bril’ and co-workers used their algorithm to investigate the heteroazeotropic distillation of ethanol-water and ethyl acetate in one- and two-column complexes with a decanting vessel. 2-Propanol dehydration with a benzene entainer was also examined. Buzzi Ferraris and Morbidelli (1981) have also looked at using a multiflash method to solve three-phase distillation problems. Buzzi Ferraris and Morbidelli presented several possible flash sequences that could be used but noted, as did Bril’ et al. (1975),that, while this method was extremely stable in approaching the solution, it usually required many iterations. The flash techniques developed by Buzzi Ferraris and Morbidelli (1981) allowed the authors to check the validity of the solutions they obtained with other methods (Newton’s and boiling point). This was necessary because these methods, as applied by Buzzi Ferraris and Morbidelli, required a prior knowledge of the three-phase region. If this was underspecified, the solution obtained would be incorrect, as a single liquid phase would be located at an unstable stage. The use of the flash technique to check the validity of solutions represented 0 1990 American Chemical Society
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1365 ‘n
I, n-1
Yn,i
l,n-l,i
STAGE
n
X2,n,i
L 1,n X
1,n , i
Figure 1. Three-phase stage.
a way of analyzing the stability of the calculated liquidphase compositions. It is interesting to note that the algorithm by Bril’ et al. (1975) requires an initial estimation of equilibriumphase compositions, yet no mention of prior knowledge of the three-phase region is discussed by the authors. The reason for this may lie in the method used to calculate the LLE where Bril’ and co-workers (1973d) accounted for the possibility that the mixture may be homogeneous. Pucci et al. (1986) more recently combined an isenthalpic flash with a stage-to-stage algorithm where the MESH equations at the stage were solved simultaneously by a Newton-Raphson method. The idea appears to essentially automate the procedure discussed by Buzzi Ferraris and Morbidelli (1981). The isenthalpic flash basically acts as a stability test, indicating when the liquid phase is unstable. To save computation time, Pucci et al. also allowed specifications of a region in which phase splitting would be likely to occur. Outside this region, only a single liquid phase was considered. The procedure used by Pucci and co-workers involves the examination of the liquid at the temperature obtained from a standard two-phase isenthalpic flash. Rather than applying a direct stability test, Pucci et al. opted for solving the isoactivity criterion. The authors noted tha this method suffers from the strong attraction to the trivial solution and, therefore, proposed a technique based on infinite dilution activity coefficients to initialize the LLE calculation. If a set of phase compositions is found from this starting point, the mixture is considered to be three phase and is appropriately flashed to produce the corresponding vapor composition. If no LLE solution is found, the mixture is considered to be stable, and the previous two-phase isenthalpic flash results are used. This method was employed by Pucci et al. to simulate the dehydration of a butyl-beer mixture and was compared with some limited experimental data obtained from a 5-cm-diameter glass column with 60 sieve trays. Approximate Methods. Guo (1986) modified the multiflash approach to investigate systems containing water and hydrocarbons. The essential features of the algorithm are similar to that of Pucci et al. (1986). A flash calculation is used to determine if the liquid phase is unstable, and this is combined with a stage-to-stage calculation where stage equations are solved simultaneously by using a Newton-Raphson method. The stability check
used is based on the comparison of the vapor-phase fugacity of water resulting from the flash calculation and the liquid fugacity of the water phase. Guo’s approximation assumes the water phase is pure and hence concludes the mixture is unstable if f’,’ f k w (12) Although not cited by Guo (1986), this type of stability method was put forward by Buzzi Ferraris and Morbidelli (1982) as part of their approximate treatment of systems in which one of the two liquid phases could essentially be regarded as a pure component. Buzzi Ferraris and Morbidelli (1982) noted that, in many applications when phase separation occurs, each component is soluble almost only in one of the two immiscible liquid phases. This allows one of the phases to be considered pure, which simplifies the calculation procedure considerably. The authors used this technique to solve two examples: an 18-stage separator fed with acrylonitrile-acetonitrile-water and a 12-stage column fed with a mixture of propanol-butanol and water. Large savings in computer time are shown for the approximate calculation, which produces acceptable results when the solubility in the water phase is small. Other approximate methods proposed to cut down on computer time have been presented by Walraven and van Rompay (1987), who utilized the inside-out convergence method in a short-cut algorithm, and by Niedzwiecki et al. (1980),who introduced the mixed K-value approach to allow modification of an existing single-liquid-phase program that used the Newton-Raphson technique. The Niedzwiecki et al. approximation restricts the calculations to systems with distribution of only one or a few components over both liquid phases. Global Newton-Raphson and Mixed K-Value Method. The global application of the Newton-Raphson method to all the MESH equations provides the opportunity to correct simultaneously all the variables and does away with the problems encountered with tearing methods where the system of nonlinear equations is broken up into smaller subgroups. One advantage is that the NewtonRaphson method converges on the solution quadratically. The application of the Newton-Raphson method to the simulation of single-phase distillation problems has been well reviewed by Vickery and Taylor (1987). The paper most cited as introducing the simultaneous correction technique was written by Naphtali and Sandholm (1971), but as Vickery and Taylor point out, the method can be traced back to Whitehouse (1964). The scaling up of the single-liquid-phase case to accommodate the possibility of a liquid-phase split has been addressed by Buzzi Ferraris and Morbidelli (1981). They showed that the set of equations describing three-phase distillation can be grouped stage-by-stageto give the desired block tridiagonal structure similar to the two-phase case of Naphtali and Sandholm (1971). In their formulation of the problem, Buzzi Ferraris and Morbidelli (1981) used all the equations, including the liquid-liquid equilibrium set. This then requires a prior knowledge of the number of three-phase stages, as the LLE calculation is solved along with all the other equations and, therefore, offers no scope for the introduction of a stability test. Although not explicitly stated by the authors, reasonably good estimates of the LLE-phase compositions for the designated three-phase stages would also have to be provided to avoid the strong attraction of the trivial single-phase solution. The authors did note, however, that, because of the high nonideality of liquid mixtures where phase separation can occur, the derivatives of the equilibrium ratios with respect to the compositions and tem-
1366 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990
perature cannot be neglected if the convergence qualities of the Newton-Raphson method are to be preserved. To avoid the problem of needing advance knowledge of the three-phase region, Niedzwiecki et al. (1980) introduced the mixed or pseudo-K-value concept. This allowed the LLE calculation to be performed outside the main Newton-Raphson equations and, therefore, frees up the restrictions of having to designate which stages experience a liquid-phase split. By calculating the LLE separately, a stability test, or at the very least a search for a solution via the isoactivity equations, can be performed to decide whether the liquid phase is unstable. It is perhaps easy to understand the motivation behind developing the mixed K-value approach when one considers that by working with the averaged variables the equation set is reduced to the conventional single-phase case. The underlying concept of this approach is to treat the two-liquid-phase mixture as a single phase and to calculate thermodynamic properties defined by an appropriate averaging of the two liquid phases present. The average K value is given by
R; = y , / f ,
(13)
where the average or overall liquid mole fraction is given by Ri = SX1,' + (1 - S ) X 2 , , (14) For an ideal vapor phase,
Ri= PP
(15)
where (16) The paper by Niedzwiecki et al. (1980) is only very brief and provides no explicit equations. They do state, however, that, when this method is used in conjunction with the Newton-Raphson method, derivatives of the K values are required with respect to the total liquid-phase composition and temperature. Baden (1984)and Baden and Michelsen (1987) have also looked a t the mixed K value approach combined with the Newton-Raphson method. Baden and Michelsen (1987) worked with activities and used the following
pi = d , / f ,
(17)
which is an equivalent expression to (16). Baden and Michelsen (1987), like Niedzwiecki et al. (1980), noted that, in order to evaluate the temperature and compositionalderivatives, it was necessary to take into account that the equilibrium-phase compositions, as well as phase fraction, were affected by changes in temperature and overall composition. Baden and Michelsen posed their general Naphtali-Sandholm algorithm in terms of component liquid flow rates and presented a scheme that utilized the calculation of the isoactivity condition to arrive at the derivatives of the equilibrium activities with respect to temperature and total component flows, These exact derivatives are prerequisites for quadratic convergence of the Naphtali-Sandholm method. The overall scheme used by Bden and Michelsen appears in principle to be similar to that put forward by Niedzwiecki et al. (1980). That is, the general equations forming the framework of the standard Naphtali-Sandholm algorithm remain unchanged, and only modifications are required to the section calculating the liquid-
phase thermodynamic properties. To decide whether or not to base the K value, and its derivatives, on the mixed or standard equilibrium ratio requires an ability to predict whether the liquid phase is unstable and will form two liquid phases. Some authors (Bril' et al. 1973d; Pucci et al., 1986; Kinoshita et al., 1983b; Buzzi Ferraris and Morbidelli, 1981; Block and Hegner, 1976) used the isoactivity condition for this purpose. If a solution to this cannot be found, then the mixture is considered stable. Unfortunately, this method suffers from the strong attraction to the single-phase trivial solution. Another approach is to perform a test to decide whether or not the liquid is unstable before attempting to solve the isoactivity criterion. Baden and Michelsen performed a stability test, based on a tangent-phase stability analysis and the work of Michelsen (1982a,b), whenever a set of liquid-phase activity coefficients is required. Depending on the result of the test, the isoactivity condition can then be solved and the mixed K value and its derivatives calculated and returned to the Naphtali-Sandholm algorithm. Bubble Point Method. One of the first three-phase distillation algorithms presented in the open literature was by Block and Hegner (1976), who modified a class of tearing techniques known as the bubble point method. In general, the bubble point method uses the summation and equilirium equations to determine the stage temperature from a standard bubble point calculation. Vickery and Taylor (1987) noted that this technique was introduced by Amundson and Pontinen (1958) and has undergone significant improvements in the intervening years, the best known of which is due to Wang and Henke (1966). The Block and Hegner (1976) method of decoupling the three-phase distillation equations uses the average liquid-phase mole fractions, xi, as the iteration variables. The equations are then broken into groups, with the isoactivity condition being solved to give equilibrium-phase compositions. If no solution can be found, the mixture is considered stable. The authors make no mention of any technique used to generate starting guesses for the liquid-liquid calculations based on the compositions obtained from the previous set. Boston and Shah (1979) point out, however, that there is an inconsistency in the way Block and Hegner have formulated this approach. The LLE calculations are performed at a temperature different from the bubble point temperature, and hence, the liquid-phase compositions and the phase fractions will not be the same as if performed at the true bubble point. They will only match at the solution, and Boston and Shah (1979) argue that this would impair algorithm performance. The third and last group of equations solved give the vapor and liquid flow rates. Block and Hegner then use a Newton-Raphson correction to update the iteration variables, and the process is repeated until the mass balance residuals and average mole fraction corrections have converged. Block and Hegner (1976) solved two three-phase distillation examples. The first looks at the distillation of a l-propanol-water-butanol system where the butanol concentrates sufficiently on the lower trays to form an additional liquid phase. The effect of drawing off an aqueous phase side stream is also investigated. The second example looks at the removal of water from a saturated mixture of butanol and butyl acetate. Buzzi Ferraris and Morbidelli (1981) have also looked at the bubble point method. They proposed an algorithm similar to that of Block and Hegner's; however, they did not use the Newton-Raphson technique to update the average mole fraction variables. Instead, they optted for
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1367 the procedure proposed by Boston and Sullivan (1972), which the authors argue eliminates numerical instabilities due to differences between two numbers of the same order of magnitude. Buzzi Ferraris and Morbidelli also noted the importance of avoiding the trivial solutions when attempting to calculate the LLE values. To do this, the authors used a method based on the work of Buzzi Ferraris and Casapollo (1972), which restricts the unknown phase compositions to lie within a specified range. The obvious disadvantage of this approach is that it is system dependent. Moreover, it is necessary to know a priori the location of the three-phase trays. Kinoshita et al. (1983b) extended their two-phase modified bubble point algorithm (Kinoshita et al., 1983a) to account for a liquid-phase split. This approach is essentially the same as that of Block and Hegner (1976), although the deviation function used to generate the average liquid mole fraction corrections (via Newton-Raphson) is different. Kinoshita and co-workers reported great difficulty in avoiding the trival solution when solving the LLE subset of equations. Initially a Newton-Raphson method was used, but convergence to the equilibrium phases could not be achieved. The authors then used a minimization approach utilizing the Simplex method of Nelder and Mead (1965). This was performed repeatedly for different initial conditions, and if all results converged to the trivial solution, the mixture was considered stable. The method used to choose the initial conditions was not discussed. One further three-phase application of the bubble point method has been briefly discussed by Baumgartner et al. (1985). They used a modified Wang-Henke procedure to model ethanol dehydration with cyclohexane. Newton’s method is reported to have been used to assist in adjusting the vapor flow profile to enable satisfactory convergence. The authors also state that phase splitting was assessed by using a modified method based on the LLE algorithm of Prausnitz et al. (1980). Inside-Out Methods. Boston and Sullivan (1974) introduced a method for solving distillation problems in which the complicated equilibrium and enthalpy expressions were replaced by simple models and the primitive iteration variables, such as temperature and mole fractions, were replaced by variables within these models that were relatively free of any cross interations. This method has become known as the inside-out method. Boston and Shah (1979) extended the Boston and Sullivan (1974) method to incorporate two liquid phases. As with the original inside-out algorithm, the three-phase extension used the physical property parameters as iteration variables but approximated the average K values by using an expression developed by Boston and Britt (1979). Further, for the three-phase case, a new expression for the ratio of the liquid-phase activity coefficients was developed by using the concept of Boston and Britt (1979). Boston and Shah (1979) utilized Broyden’s method to update the vapor and liquid flow rates. A significant feature of the Boston and Shah (1979) method was the development of a phase-splitting algorithm. This algorithm bases the initial estimates of the phase split on what the authors have called the “maximum effective infinite dilution activity”. This is described further in the following sections. Boston and Shah solve several problems with the inside-out algorithm. The first example looks at a waterhydrocarbon system and the effect of withdrawing the water phase from the trays. To demonstrate the reliability and robustness of the algorithm and stability test, the
authors solve a second example, a 10-stage separator fed with acetone-chloroform and water, at temperatures up to 150 K outside the final solution. Ross (1979) and Ross and Seider (1981) also developed an algorithm based on the inside-out concept of Boston and Sullivan (1974) that included Murphree tray efficiencies. Prokopakis et al. (1981) modified this method slightly to investigate the dehydration of 2-propanol and ethanol with cyclohexane and benzene as the respective entrainers. Ross and Seider (1981) have made some modifications to the original inside-out method of Boston and Sullivan (1974) and the three-phase application of Boston and Shah (1979). As Ross (1979) points out, the flow rates are calculated by using the bounded Wegstein method (Myers and Seider, 1976) instead of the quasi-Newton method of Broyden. Ross (1979) cites storage and the number of calculations required to update the inverse of the Jacobian as reasons for using this approach. Ross and Seider use the liquid-phase-splitting algorithm of Gautam and Seider (1979). This approach differs from the Boston and Shah (1979) stability test, as the Rand method is utilized to minimize the Gibbs free energy. The introduction of Murphree tray efficiencies forced Ross and Seider (1981) to solve the material balance and equilibrium relationships simultaneously and also to develop special equations to allow the calculations of temperatures. Ross and Seider also changed the method of testing for convergence. Boston and Sullivan (1974) and Boston and Shah (1979) tested convergence of the outer loop by examining the changes in the coefficients of the physical property models. Ross (1979) notes that this makes it difficult to distinguish between convergence and effects due to damping. Ross and Seider (1981),therefore, compare the temperatures and average liquid mole fractions used to calculate the parameters for the approximate physical property models with the values calculated at the end of the inner loop. The test for convergence is thus unaffected by damping, as the conditions at the end of the inner loop are purely a function of the compositions and temperatures at which the parameters for the physical property models are calculated. Ross and Seider (1981) solved example 2 proposed by Block and Hegner (1976) and also looked at a second example by replacing the propanol of this example with ethanol. Schuil and Bool (1985) also used the inside-out method to look at these two examples. These authors, however, used the mixed K-value approach in conjunction with the General Mass Balance (GMB) package, which is described by Russell (1980). This package uses the inside-out algorithm of Boston and Sullivan (1974). They report comparable results with those obtained by Block and Hegner (1976) and Ross (1979) and point out the implementation of mixed K values required much less work than writing a new program. Schuil and Bool used the LLE calculation method described by Prauznitz et al. (1980) to account for phase splitting. Homotopy-ContinuationMethods. The last decade or so has seen increasing interest in the application of homotopy-continuation methods in solving chemical engineering problems. The role of continuation in engineering analysis has been recently reviewed by Seydel and Hlavacek (1987) and highlights the advantage of continuation-based methods when dealing with problems that have multiple solutions. The application of homotopy-continuation methods to the solution of distillation problems has received some interest in recent years (see, for example, Wayburn and
1368 Ind. Eng. Chem. Res., Vol. 29, No. 7 , 1990
Seader (1984) and Vickery and Taylor (1986)). The basic concept behind the method is the creation of homotopy, which is simply the deformation of a set of equations that are difficult to solve, f ( x ) = 0, into a set whose solution is known or easily found, g(x) = 0. The homotopy equations, H ( x , t ) , are generated by some combination of f ( x ) , g(x), and a homotopy parameter, t , such that
H ( x , t ) = H [ f ( x ) g ( x ) , t I= 0 (18) The key to this deformation is that at t = 0 the homotopy reduces to the simple system of equations, g(x), and equals the original equation f ( x ) , set, when t = 1. That is, H(x0,O) = g(x0) = 0 (19)
H(xf,l)= f ( x 9 = 0
(20) Here x represents the vector of independent variables, xo the initial x and a solution to g ( x ) = 0, and xfthe final x and a solution to f ( x ) = 0. The form of the combination of f ( x ) , g ( x ) , and t that makes up the homotopy varies from application to application, and each has its own advantages and disadvantages as discussed by Seydel and Hlavacek (1987). Recently, Kovach and Seider (1987a) extended the homotopy-continuation method of Allgower and Georg (1980) to simulate heterogeneous azeotropic distillation columns. This method uses convex linear homotopy H ( x , t ) = t f ( x ) + (1- t ) g ( x ) = 0 (21) which is so named because H(x,t) is linear in the parameter t. Kovach and Seider used a g ( x ) function that describes a linear decay of the errors in the components of f ( x ) from their initial value to zero. This occurs with the so-called Newton's homotopy, which is created when g(x) is set as g(x) = f ( x )
- f(xO)= 0
(22)
To simplify matters, Kovach and Seider worked with the homotopy parameter X = 1 - t , which reduces the homotopy equations to H(x,X)= f ( x ) - Xf(x0) = 0
(23)
where H(xo,l) = f ( x o ) - f ( x o ) = 0
at t = 0, X = 1
(24)
H(rf,O)= f ( x 9 = 0 at t = 1, X = 0 (25) The locus of roots to eq 23 from X = 1, with x = xo, to X = 0 and x = xf defines the homotopy path or curve. Seydel and Hlavacek (1987) note that continuation means the tracing of this solution curve from the initial starting point to the final solution. Given the extent of the problem, only the final result may be of interest, whereas, if desired, continuation also provides information on intermediate results encountered on the way to this solution. This is particularly useful in the regions of multiple solutions. To track the homotopy of eq 23, continuation methods differentiate the equations and convert them into initial value problems. Kovach and Seider (1987a,b) noted that difficulties can arise with the tracking at turning or limit points in the homotopy path. Kovach and Seider further noted that, when the solutions of f ( r ) = 0 are close to the limit points, Newton-based methods become unreliable, as the Jacobian matrix of f ( x ) becomes singular at the turning point. Kovach and Seider (1987a,b) point out that many applications of homotopy-continuation used in the literature only seek to locate the final solution, xf,and do not necessarily follow the homotopy path too closely. The simulation of azeotropic distillation, however, presents a
problem in which limit points are frequently encountered, and the path must be closely followed so as to avoid cycling and to locate multiple solutions. To achieve this, the authors extended the algorithm of Allgower and Georg (1980),which used a Euler predictor step, in the direction of the tangent vector, and a Newton corrector step, in the direction orthogonal to the tangent vector, to speed the return of x to the homotopy path. The original Allgower and Georg algorithm utilized only two Newton correction steps. Kovach and Seider noted that this is not enough to track the homotopy path around limit points and that in azeotropic distillation the curve can be complex and the trace may end up on another branch. To avoid this, Kovach and Seider repeat the Newton corrections until quadratic convergence is obtained and convergence tolerances are satisfied. This, however, means that limit points are approached more closely, and hence, uninvertible near-singular Jacobian matrices are encountered before the path-following parameters can be changed to avoid the singularity. Kovach and Seider, therefore, proposed a series of tests designed to indicate the approach of a limit point, enabling the path-folowing parameter to be changed to the variable having the greatest rate of change during the preceding step. To solve the distillation problem, Kovach and Seider ordered the equations in the same manner as Waybum and Seader (1984) with material balances first, followed by the energy balances and then the equilibrium equations. After the iteration variables are updated by either the Euler predictor or Newton corrector steps, the stream enthalpies are calculated and the liquid phases are checked for phase splitting by using a LLE homotopy algorithm. This represents an important element of the simulation method, as the authors found that the original method of Gautam and Seider (1979) did not always produce reliable results. If a stable phase is detected, the second-phase flow rate is added to the first and dropped from the iteration vector. Kovach and Seider point out that this addition and deletion of the second liquid phase to the variable vector results in no discontinuities in the mass balance residuals. Furthermore, the correction to the component flow rates in the second liquid phase are determined by using the full set of linearlized equations and when the overall liquid composition approaches the binodal curve, the impact of the VLE equations is minimal. These properties avoid instabilities associated with the addition and deletion of the second liquid phase. Kovach and Seider used their homotopy-continuation algorithm to solve the ethanol-water-benzene azeotropic distillation problem proposed by Prokopakis and Seider (1983). They detected five steady-state solutions over a very narrow range of operating conditions and showed that three solutions existed over a slightly larger range. A three-phase application was demonstrated by simulating a tower used to dehydrate sec-butanol. Kovach and Seider (1987b) performed experiments to obtain plant data which they compared with solutions obtained by the homotopy-continuation method. It is perhaps interesting to note that, during the process of trying to simulate the plant data, the authors obtained two steady-state solutions. The first was easily found by using a Newton-based iteration method and had a single liquid phase on all trays. However, the product compositions predicted by this solution were not in agreement with the observed plant values. The second solution was reported to be only found with parameterization by using the homotopy-continuation method. This solution predicted two liquid phases on 70% of the trays, and product compositions were in close
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1369 agreement with the experimental values. It should be noted, however, that Kovach and Seider (1987b) did not measure plate compositions but were constrainted to analyzing the product streams and measuring several key temperatures. Therefore, there is no way of knowing if the actual column was operating with 70% of the trays in the three-phase region.
Finite Element Collocation Methods Swartz and Stewart (1987) recently extended the orthogonal collocation approach of Stewart et al. (1985) to cover multiphase distillation. The original work by Stewart and colleagues (1985) involved representing the fractionation system as a series of interconnected models or finite elements, each of which correspond to a physical section of the column. The motive behind this approach is to reduce the order of the modeling equations and avoid the large stiff systems of nonlinear algebraic and differential equations associated with the so-called full-order simulations. The basic idea of this technique is to approximate the states in each module (finite element) by Lagrangian polynomials whose nodal values are determined by orthogonal collocation. This approximation reduces the problem to one of a lower order. The polynomials are fitted to the process model at points designed to give a good approximation of the full set of stages. Stewart et ai. (1985) presented a new least-squares criterion for the selection of the grid points, with the resulting collocation method shown to be accurate, converging to the full stagewise solution as the number of grid points is increased. Swartz and Stewart (1986) extended this approach to investigate the optimal design of distillation systems. This was achieved by lifting the restriction on module length and treating it as a continuous variable, thus extending the discrete solution set. The authors point out that distillation design problems are consequently reduced to standard nonlinear programming problems that greatly facilitate optimization calculations. With the extension of continuous model length variables, Swartz and Stewart (1987) expanded the method to include multiphase systems. This was accomplished by using a separate finite element to represent each multiphase region. This enables a distinction to be made between the physical column sections, as determined by feed and draw-off streams, and the collocation modules, where sections containing phase discontinuities are represented by the corresponding number of finite elements. The adjustable module lengths are treated as continuous variables and their s u m constrained to be consistent with the physical dimensions of the column. Swartz and Stewart (1987) show how the module boundary location can be considered as boundaries of the multiphase region and included as a variable. The location of these boundaries is such that the additional liquid phase is either just beginning to form or has just disappeared. Swartz and Stewart noted that this newly forming or vanishing phase is required to be in equilibrium with the other liquid phases and hence is analogous to the bubble point condition. With this equilibrium condition, the authors developed a series of equations describing either the formation or disappearance of a new phase. The solution procedure adopted by Swartz and Stewart (1987) involved obtaining the initial distribution of breakpoints by the calculation of a two-phase solution to the problem. The stability test of Fournier and Boston (1981), which is based on the work of Boston and Shah (1979), is then applied to the liquid phase at each collocation point and at the liquid entrance and exit positions
for each module. Column sections containing phase discontinuities are then subdivided into the indicated modules. Swartz and Stewart used interpolation to give the states at new collocation points, and the multiphase system of model equations was then solved by using a damped Newton method. Each adjustable module length is bounded from above and below the correspondingnumber of collocation points. The authors demonstrated their technique by solving example 2 of Block and Hegner (1976) and the acrylonitrile-acetonitrile-water problem proposed by Buzzi Ferraris and Morbidelli (1982). Summary of Attempted Three-phase Distillation Examples. Table I presents a summary of the systems and distillation parameters that make up the three-phase distillation applications found in the literature. In all but one case (number 12), the separations involve water. By far, the most often attempted example has been the butanol-water-propanol system (number 3) originally proposed by Block and Hegner (1976). It would appear that most of the simulation algorithms discussed above have been used to find a solution for this example, largely based on the predictions of NRTL thermodynamics. The earlier work of Pratt (1947) (number 10) represents the group of workers who, around the period after World War 11, presented various azeotropic column designs based on hand calculations, which accounted for liquid-phase splitting. Pratt (1947) looked at two examples that showed that the top two stages of the proposed columns would lie within the three-phase region. These columns do not appear to have been simulated by using any of the modern methods. Many of the examples shown in Table I look at the effect of withdrawing a side stream (numbers 3,6,7,8,12, and 13). This usually takes the form of removing the aqueous phase from the column once it is detected. Several three-phase azeotropic columns have been stimulated by using most of the simulation methods (numbers 1, 2, 11, 14, and 15). These examples vary depending on the systems involved and whether or not a mixed reflux or only the entrainer phase from the decanter is sent back to the column. Other variations in the azeotropic examples involve changing the feed plate location, product rates, and reflux ratios. Most of the earlier examples investigated appear to have been limited to ternary systems, possibly because of the difficulty in obtaining appropriate equilibrium data for mixtures containing more components. Recent studies have investigated multicomponent systems with up to four and five components, with the equilibrium data either being specially measured (Kovach and Seider (1988)-for use in number 15) or generated by using the group contribution method of UNIFAC (Baden (1984), number 12).
Liquid-Phase Splitting A fundamental thermodynamic criteria for equilibrium between phase is the equality of chemical potentials for each component in all phases. That is, for two phases, I and 11, 1; = /.ti1
i = 1, ..., n,
(26)
It can be shown (see, for example, Arlt et al. (1982)) that this condition reduces to the well-known isoactivity criterion for liquid-liquid equilibrium, which states uf = up
i = 1,
..a,
n,
(27)
The activity of component i (ai)can be expressed in terms of an activity coefficient (ri),which is defined by the relationship
1370 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 Table I. Summary of AttemDted Three-Phass Distillation Examsles feed plate and flow q value no. system 1 ethanol (1) F, = 100 mol/h q ' = 1.0 water (2) F2 = 100 mol/h ethyl acetate (3) q = 1.0 F6 = 100 mol/h q = 1.0 2
2-propanol (1) water (2) benzene (3)
F6 = 100 mol/h q = 1.0
F6 = 100 mol/h q = 1.0 3
butanol (1) water (2) propanol (3)
F4 = 50 mol/h q = 1.0
F4 = 50 mol/ h q = 1.0
7
butanol (1) water (2) butyl acetate (3) butanol (1) water (2) ethanol (3) propylene (1) benzene (2) n-hexane (3) water (4) acetone (1) chloroform (2) water (3)
FI = 50 mol/ h q = 1.0
F7 = 50 mol/h q = 1.0 F3 = 200 mol/h q = 1.0
F3 = 50 mol/h q = 0.0
F4 = 100 mol/h q = 1.0 butanol (1) water (2) butyl acetate (3) acrylonitrile (1) acetonitrile (2) water (3) acetonitrile (1) trichloroethylene (2) water (3)
F3 = 160 mol/h q = 1.0
F12= 100 mol/h q = 1.0 Fl = 104.9 mol/h q = 1.0
F4 = 100 mol/h q = 1.0 F , = 32.24 mol/h q = 1.0 F, = 100 mol/h q = 1.0 11
benzene (1) water (2) ethanol (3)
F, = 42.101 mol/h q = 1.0 F6 = 100 mol/h q = 1.0
12
13
14
propane (1) butane (2) pentane (3) methanol (4) hydrogen sulfide (5) water (1) acetone (2) ethanol (3) butanol (4) ethanol (1) water (2) cyclohexane (3)
sec-butyl alcohol di-sec-butyl ether (2) water (3) methyl ethyl ketone (4) butylenes (5)
23
= 0.500 = 0.415 = 0.085
21
= 0.500
2, 2;
0.415 2 3 = 0.085 zI = 0.75 22 = 0.25 23 = 0.00 t l = 0.75 22 = 0.25 2 3 = 0.00 z I = 0.13 z2 = 0.65 23 = 0.22 tl = 0.13 22 0.65 23 = 0.22 t, = 0.24 z2 = 0.30 23 = 0.46 z1 = 0.03 z 2 = 0.75 23 = 0.22 t i = 0.45 22 = 0.30 t g = 0.20 z4 = 0.05 tl = 0.167 z 2 = 0.333 23 = 0.500 2 , = 0.6 22 = 0.2 z3 = 0.3 z1 = 0.500 22 = 0.3125 23 = 0.1875 z1 = 0.40 22 = 0.35 23 = 0.25 z1 = 0.450 z2 = 0.521 23 = 0.029 z1 = 0.6848 22 = 0.0002 23 = 0.3150 z1 = 0.450 z 2 = 0.521 23 = 0.029 21 = 0.89985 z2 = 0.000 17 23 = 0.09998 z1 = 0.7376 z2 = 0.0403 23 = 0.2221
bottoms flow 91.5 mol/h
distillate flow 8.5 mol/h
91.5 mol/h
8.5 mol/h
75.0 mol/h
25.0 mol/h
75.0 mol/h
25.0 mol/h
21.0 mol/h
29.0 mol/h
6.53 mol/h
29.0 mol/h
34.84 mol/h
15.6 mol/h (aq decanter phase)
36.0 mol/h
14.0 mol/h
130 mol/h
70.0 mol/h
110 mol/h
40.0 mol/h
120 mol/h
40.0 mol/h
39.2 mol/h
60.8 mol/h
66.1 mol/h
138.8 mol/h
89.59 mol/h
42.65 mol/h
84.501 mol/h
57.6 mol/h
8235.258 mol/h
115.0 mol/h
27.1 kg/h
17.6 kg/h
1.57167 kmol/h
2.36293 kmol/h (to decanter)
69.97 kmol/h
138.94 kmol/h (leaving decanter)
2 %=
22
= 0.000 = 0.11
23
= 0.89
21
F12= 8480.258 mol/h
z1 =
q = 1.0
22
23
0.00048
= 0.00426 = 0.00267
0.977 19 0.01540 z1 = 0.0504 w t 22 = 0.0136 wt 2 3 = 0.0043 wt L, = 0.0317 wt z1 = 0.8555 zz = 0.1445 z, = 25 =
F1, = 100 kg/h q = 1.0 Fl = 2.7287 kmol/h q = 1.0
23
F1 = 1.0259 kmol/h q = 1.0
15
feed composition
F7 = 208.91 kmol/h q = 1.0
= O.oo00
z1 =
0.0871
= 0.0089 23 = 0.9040 21 = 0.3490 22
0.0072 0.6379 0.0013 = 0.0046
zz = z3 = z4 = 25
(1) Block and Hegner (1976). (2) Buzzi Ferraris and Morbidelli (1981). (3) Buzzi Ferraris and Morbidelli (1982). (4) Ross and Seider (1981). (5) Kinoahita et al. (1983a,b). (6)Schuil and Boo1 (1985). (7) Walraven and van Rompay (1987). (8) Swartz and Stewart (1987). (9) Ross (1979). (10) Boston and Shah (1979). (11) Bril et al. (1975). (12) Pratt (1947). (13) Baden (1984). (14) Pucci et al. (1986). (15) Baumgartner et al. (1985). (16) Kovach and Seider (1987a,b). a
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1371 no. of plates (excl condenser; incl reboiler) 35 (+ decanter)
none
35 (+ decanter)
none
12.0 (0% phase)
10 (+ decanter)
(a) 12.0 (b) 25.0 (c) 13.8 (org phase) 3.0
liquid-phase thermo NRTL
preeeure, atm 1.00
rep 11
total
NRTL
1.00
11
none
total
NRTL
1.00
11
15 (+ decanter)
none
total
NRTL
1.00
11
11
none
total
NRTL
1.00
1-8
3.0
11
S, = 14.47mol/h
total
NRTL
1 .00
1,4
0.91
7 (+ decanter)
none
total
NRTL
1 .00
1
3.0
11
none
total
modified van Law
1.00
4,6,9
2.5
5
(a) none (b) aq phase on any tray
total
Margules
1.00
10
2.5
10
(a) none (b) aq phase on any tray
total
UNIQUAC
1.00
10
2.5
5
(a) none (b) aq phase on any tray
total
NRTL
1.00
10
27.1
18
none
total
NRTL
1.00
3, 8
0.614
12
none
total
exptl
1.00
12
2.3
13
none
total
exptl
1.00
12
5.47
25
none
partial
UNIFAC
1.00
13
7.0
20
Ss= 50 mol/h; S, = 60 mol/h
partial
UNIFAC
4.75
13
0.43
29
Sm
total
NRTL
1.00
14
4.5
32
none
total
UNIQUAC
1.00
15
1.7 (reflux from decanter)
32
none
total
UNIQUAC
1.27 (top); 1.39 (bottom)
16
reflux ratio 25.0 (org phase)
side streams
55.3 kg/h
condenser type total
1372 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990
i = 1, ..., n,
ai = yixi
(28)
Equation 27, therefore, becomes $xi =
i = 1, ..., n,
ryxp
(29)
where x : represents the mole fraction and 7;the activity coefficient of component i in phase I. The activity coefficient is, in turn, a function of composition and temperature Ti = Yi(xl,xZ,.*.rxnc,T)
(30)
and can be calculated by using an appropriate model. The liquid-liquid flash problem that is required to be solved for three-phase distillation is to calculate the equilibrium phase compositions and phase fractions for a given overall composition (zi) and system temperature (T).If two equilibrium phases exist for these conditions, then the problem involves 2n, + 1 equilibrium and mass balance equations:
xi$ $Xi
- xpyy =
+ (1 - s)xf
0
i = 1, ...) n,
= zi
i = 1, ..., n,
(31) (32)
n.
2x/ =1 i=l
(33)
The unknowns are the phase compositions ( x f and x;’) and the fraction of the overall mixture occupied by phase I (s). The summation of phase I mole fraction (eq 33) has been included arbitrarily. This then excludes the corresponding sum for phase I1 as it is implicit in the equation set since the feed mole fractions (q)are constrained to sum to one. In principle, it is possible for systems that exist as two liquid phases a t equilibrium to find values of xf,xi’, and s that satisfy eq 31-33. However, this is complicated considerably by the fact that there are multiple solutions to this equation set. Classical thermodynamic analysis states that a stable equilibrium will be attained when the total Gibbs free energy of the system attains a global minimum. A necessary but insufficient condition for this is that the Gibbs surface is stationary a t the equilibrium point. The criterion of (26) is derived from this stationary condition only and, therefore, represent a necessary and, unfortunately, insufficient test for stable equilibrium. The stationary condition will admit local minima, as well as maxima in the Gibbs surface, as equilibrium states, and therefore, (26) could detect an equilibrium that is far from the global minimum in Gibbs free energy. This phenomenon is linked to the stability of the proposed equilibrium phases, and it is important to recognize that, for threephase distillation, metastable phases represent unstable phases. Adding to the problem of computing LLE is the trivial solution. Heidemann (1983) notes that when two phases are described by the same thermodynamic model, as is the case for LLE, it is possible that the equilibrium computation may fail and converge to the so-called “trivial solution”, with both phases having the same composition and density. When this happens, the sizes of the phases are indefinite, and the phase fraction parameter(s) can take any value. Further, since the phases are identical, the free energy is unchanged when mattter is transferred from one phase to the other. The trivial solution poses a severe problem in LLE calculations. In practice, it is found that the domain of convergence for the trivial solution is large, and whether this solution or the correct solution is reached depends
upon the computational procedure used and the way the phase split is initiated. There are two broad approaches that can be taken when calculating LLE. The first involves searching for a solution to the equation set (31)-(33), or some other equivalent formulation, from an initial estimate of the phase split. The initialization must be such that the calculation technique is provided with a good opportunity to converge on the stable equilibrium phases, if they exist, and avoid the trivial solution. If the trivial solution is found from this or a series of starting points, the mixture is considered to be homogeneous. The second approach involves performing a stability test to decide whether or not the mixture is capable of existing as two liquid phases before attempting to solve the phase split equations. The first approach has been used by many authors who have employed a variety of different convergence and initialization schemes to arrive at a solution. The computational merits of many of these schemes have been critically evaluated by Ohanomah and Thompson (1984) and more recently by Swank and Mullins (1986). One of the first techniques used was by Henley and Rosen (1969), who arranged the equations so that they could be solved in a direct iteration scheme. While this approach works well for VLE calculations, the solution for LLE is approached slowly, and some authors have investigated various acceleration methods to promote convergence. Null (1970) used successive substitution with Newton’s method, whereas Bril’ et al. (1973d) opted for an analytical criterion on the initial choice of liquid-liquid K values. Prausnitz et al. (1980) used a scheme involving initial convergence with Newton’s method, which was accelerated by using the technique described by Wegstein (1955). Soares et a1 (1971) have looked at a Newton-Raphson method for solving the equation set, and Sumberova et al. (1977),and more recently Kovach and Seider (1987a),have implemented homotopy-type methods that imbed a parameter in the equation set. Other workers have formulated the problem to allow a search for a minimum value of Gibbs free energy, given an initial estimate of the phase split. Guffey and Wehe (1972),Dluzniewski and Adler (19721, Marina and Tassios (19731, Heidemann (1974), Heidemann and Maudhane (1975), Gautam and Seider (19791, Michelsen (1982a,b), and Soares et al. (1982) have all looked at various ways of using minimization. Some of the numerical methods suggested in the literature include Murray’s method (Murray, 19721, the Broyden-Fletcher-Goldfarb-Shanno method (Dennis and More, 19771, and the Rand algorithm (Gautam and Seider, 1979). The key to the success or failure of all the above methods is the initialization of the phase split. Many techniques require an advance knowledge of the so-called “key components” that initiate the split. Nelson (1987) argues that in many cases this does not cause a serious problem because “a process engineer almost always knows which components in his slate will dominate a second liquid phase”. The usual procedure is to set the initial estimate of the two liquid phases to be pure (or nearly pure) in the key components. This hopefully represents a starting point far enough away from the feed composition to avoid converging on the trivial solution. There have been some attempts to automate the procedure of selecting the key components. Prausnitz et al. (1980) select the two components, present in the system in significant quantities, that have the lowest binary sol-
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1373 ubilities. The method then sets the initial phases to be almost pure with only 2% solubility of the opposing components in each phase. Boston and Shah (1979) devised an initial guess algorithm that begins by identifying the two key components based on what the authors call the “maximum effective infinite dilution activity”. This is defined as the product of the component’s mole fraction and infinite dilution activity coefficient. The first phase indicator is selected as the component with the largest effective infiiite dilution activity. The second is then chosen by identifying the component with the largest infinite dilution activity in a binary mixture with the first phase indicator. The components are then distributed between the phases according to the mass balance and equilibrium requirements of the system. Gautam and Seider (1979) also developed an initial split algorithm which is similar to that of Boston and Shah (1979). Instead of basing the decision on the infinite dilution activity, Gautam and Seider (1979) use the component’s actual activity. The first key component is chosen as the component with the largest activity in the singleliquid-phase mixture. The second key component is then selected as the component with the largest activity in the binary mixture with the first. Correct proportionality is maintained by using the relative concentrations of the components in the original mixture. The two key components are then distributed between the trial phases by neglecting the other components and solving the binary isoactivity and material balance equations. The remaining components are then distributed by using the material balance and isoactivity equations. The underlying problem with methods that only search for a solution to the equation set (31)-(33) is that there is no a priori information available on the true number of phases present at equilibrium. It is, therefore, conceivable that a method may fail to find the correct two-liquid-phase solution by always converging on the trivial solution and hence incorrectly assume the mixture is homogeneous. This problem can become even more complicated when the trivial solution is recognized as the true solution in the metastable region. Gautam and Seider (1979) removed this problem by checking the Gibbs free energy of all the feasible trial phases created by using their initialization procedure described above. The two trial phases that have the lowest Gibbs free energy are taken as candidates to replace the original mixture, and the Gibbs free energy of the split phase system is compared to the initial mixture free energy. If no decrease in Gibbs free energy in obtained after a few iterations with these trial phases, a new trial phase is selected with the next lowest Gibbs free energy and the process repeated. If no decrease is achieved after examining all combinations of source and trial phases, the phase splits are considered unsuccessful and the mixture is considered stable. This is similar to the second approach for the LLE problem where a stability test is used to decide whether or not phase splitting will occur before converging the phase split equations. An important part of any stability test, however, is the ability to identify correctly metastable regions as unstable. Some approximate and, therefore, restricted techniques have been proposed. Maurer and Prausnitz (1979) used the convexity condition to exclude intrinsically unstable regions by examining the condition of the determinant of the matrix of second composition derivatives of G. An approximate procedure is then used to determine the bi-
nodal curve by restricting the calculations to systems where only one binary pair has a miscibility gap. Metastability is accounted for by stepping toward the midpoint of the binary immiscible pair and testing for intrinsic instability. If, after the step, the mixture is unstable, the original mixture is considered to be metastable and is split into two phases. If the mixture remains intrinsically stable, the original mixture is considered to lie outside the metastable region and is thus identified as stable. Maurer and Prasnitz (1979) developed a procedure for choosing the appropriate step length; however, the method would appear restricted to mixtures with only a few components. Pham and Doherty (1990) are soon to publish a very simple technique for determining stability based on the system’s temperature. The authors realized that, for some of the systems encountered in three-phase distillation, there is a unique maximum temperature at which a given two-liquid-phasemixture ceases to exist as a heterogeneous system. Above this temperature, the binodal curve has contracted, leaving the mixture outside the unstable region. Pham and Doherty (1988) discuss the computational aspects of finding this maximum temperature, which then allows for a comparison with the system temperature to determine stability. The authors note, however, that the technique is not general and can only be used for certain types of mixtures. The more fundamental approach to stability testing is to examine the Gibbs free energy of the system to see if it can be reduced. This is essentially the approach used by Gautam and Seider (1979) and a similar method was adopted by Boston and Shah (1979) and later by Fournier and Boston (1981). These authors used the initial phase split algorithm of Boston and Shah (1979) to calculate the Gibbs energy of a split phase system, which was compared to the original homogeneous mixture free energy. If after 10 iterations of their algorithm the free energy was not reduced, the original mixture was considered stable. Arguably the most useful stability technique to appear in the literature is the tangent plane stability test. This method offers scope to test for stability based on a system’s Gibbs free energy and to simultaneously generate initial conditions for the phase split calculations should an unstable mixture be detected. This technique has been used in this work and is discussed below.
Tangent Plane Stability Analysis Heidemann (1983) points out that a phase whose stability is in question must be assured to be a t the lowest possible Gibbs free energy before instability can be ruled out. A concise mathematical way exists to express this criterion. If the questionable phase has mole fractions zi and chemical potentials poi, then this phase must be compared with every other possible trial phase having the whole range of mole fractions yi and chemical potentials piCy). If, for all other phases at the given temperature and pressure of the system, nc
F(Y)=
CYi[pi(Y) - 1LpI 2
i=l
0
(34)
then in the original mixture zi is stable. The basic significance of (34) is that F(y) represents the vertical distance from the tangent hyperplane to the molar Gibbs energy surface at composition z to the energy surface at the test phase composition y . This is illustrated in Figure 2 for a binary mixture. The stability of the original mixture (2) requires that, for all possible values of y, the tangent at z at no point lies above the Gibbs curve. That is, if a y can be found where this condition is violated, the
1374 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990
Figure 2. Tangent plane stability analysis.
mixture is unstable as shown in Figure 2. An important aspect of this criterion is that metastable regions are recognized as unstable. Baker and Pierce (1981) proved this theorm, which was originally derived by Gibbs, but gave no numerical implementation for the method. Recently, Eubank and Barrufet (1988) have investigated some simple implications for a binary mixture; however, the most extensive treatment of this stability test is due to Michelsen (1982a). Michelsen presented a variety of techniques for locating trial compositions of y that would indicate instability, and the method of stationary point comparisons has been implemented in this work. The basic idea behind this technique is that stability can be checked by evaluating the left-hand side of (34) at its stationary points. If the condition is satisfied for all stationary points, then F(y) will be satisfied because all minima of F(y) are located in the interior of the permissible region (Le., yi 2 0 and Cyzlyi = 1). The key advantage of this technique is that the method also produces excellent initial estimates for the phase split calculations as reported by Swank and Mullins (1986). To formulate the method for liquid-liquid calculations, the chemical potential is expressed as pi = p y
+ RT In (xiyi)
(35)
where p y is the component chemical potential in the standard state. This allows (34) to be rewritten in terms of a new function, g(y), since temperature is constant: F(Y) g(Y) = nm -
t Figure 3. Tangent plane stationary point method.
This is illustrated for an unstable binary system in Figure 3. The original mixture is thus stable provided K is nonnegative at all stationary points. Further, K will equal zero at the trivial solution ( 2 ) as shown in Figure 3. Michelsen (1982a) notes that by introducing the new variables
Yi = y i c K
i = 1, ..., n,
eq 39 can be written as In Yi + In yi - hi = 0
i = 1, ..., n,
(40) (41)
The new independent variables, Yi,can then formally be interpreted as mole numbers, giving Yc =
By choosing the same standard states (usually the pure liquid), (36) collapses to n,
g(Y) = Cy,Lln (YlYJ 1=1
-
h,I
(37)
where h, = In iz,rP) (38) Michelsen (1982a) shows that finding the stationary points of (37) reduces to solving i = 1, ..., n, In (y,y,) - hi = K (39) where K is independent of i. Michelsen notes that at the stationary points of g(y) the tangent hyperplane to the energy surface is parallel to the hyperplane at z, with K representing the vertical distance between the two planes.
Yi n, i = 1, '.., n,
c Yj
(42)
j=1
In summary, the method reduces to finding the solutions at (411, and the stability is verified provided a t all stationary points K 1 0, which corresponds to xy;lYi I1. That is, ne
(1) if at all solutions of (41) CYi I 1, i=l
the mixture is STABLE (43) or nc
(2) if at some solution(s) of (41) C Y i > 1, i=l
the mixture is UNSTABLE (44)
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1375 points for the phase split solutions I and 11. Figure 4c also shows the stationary saddle point at the trivial solution
A) GIBBS FREE ENERGY OF MIXING
(2).
Michelsen (1982a) examines a variety of techniques to implement the tangent plane method including a direct iteration scheme that solves (41) by using In Ylt+1) = hi - In it) t = iteration counter (45)
~p)
where is calculated by using y f ) ,which was evaluated by (42) a t the previous iteration (Le., Y?). After each iteration step, Michelsen calculates the value of a modified distance function, g*(Y), and a convergence variable, r(Y):
PROPANOL
nc
B) CONTOURS OF G(Y) FOR A STABLE MIXTURE
----
g*(Y) = 1
+ i = l Yi[ln Yi + In yi - hi - 11
POSITIVE
r(Y) =
NEGATIVE
2g* -
P
(46) (47)
where
WATER
HEPTANE
PROPANOL C ) CONTOURS OF G(Y) FOR UNSTABLE MIXTURE
- POSTIVE NEGATIVE
---aTIE-LINE .
WATER
.
HEPTANE
Figure 4. Ternary example of tangent plane stability analysis.
If negative stationary points for g(y) are detected, then these positions represent excellent initial estimates of the phase split. This can be seen in the simple binary case illustrated in Figure 3 where the stationary points Cysp) are close to the final split compositions I and 11. The initial estimates of the phase compositions are computed with (42) and the values for Yi calculated at the solutions to (41). Swank and Mullins (1986) found that in all cases tested these estimates provided reliable starting points for phase split calculations. The significance of the tangent plane stability test for multicomponent mixtures is illustrated in Figure 4 for a mixture of n-heptane, 1-propanol, and water. Figure 4a shows the Gibbs free energy of mixing surface for a hypothetical single phase, and it is clear that for low concentrations of propanol a single liquid phase would be unstable. Indeed, for very low values of propanol, the Gibbs free energy of mixing is greater than zero. Figure 4b shows the contours of gCy) for a stable mixture composition z. The values of g(y) are positive everywhere since the tangent plane always lies below the surface. The only minimum occurs at the trivial solution, z. Figure 4c, on the other hand, shows the contours of g(y) for an unstable mixture. The values of gCy) are negative in some regions, as the tangent plane lies above the Gibbs surface. Two negative minima are apparent Cysp), and this figure illustrates the good estimates provided by the stationary
Michelsen points out that, as Y approaches the trivial solution, r approaches one. The technique used by Michelsen was to converge simultaneously on the solutions to (41) by starting with as many initial estimates as there are components in the system. Each initial guess for the trial phase is taken to be a pure component. The value of r is then calculated at each iteration, and if its value exceeds 0.8, the solution is considered to be converging on the trivial solution, and the search is abandoned. The approach used in this work differs slightly from that above and uses the Newton-Raphson method to obtain solutions to (41) from two internally generated starting points. The values of Y(t+l)are calculated by using Y(t+l)= yt + AY (49) where AY represents the correction vector for Yt calculated from (50) JAY = -D Here J is the Jacobian matrix housing the partial derivatives of (41) and D is the vector of the residuals. That is, Di = In Yi + In y i - h, (51) Equation 50 is solved by standard Gaussian elimination to give the correction vector, which is used in (49) to generate new values for the iteration variables (YJ. The residuals are then calculated by (51), and if a convergence condition is not met, a new iteration is initiated. Rather than performing the above procedure from a pure trial phase for each component, two starting points were developed to provide a good opportunity of converging on the desired stationary point roots. In most cases, each of the liquid phases is dominated by one of the components. This is the idea of the key components discussed earlier. Since the stationary points are usually in close proximity to the phase-split compositions, identification of the key components would reduce the search from all the components to only two. Thus, searching for the roots of (41) can be considered similar to the LLE calculation itself, and therefore, any of the previously discussed methods of phase initialization would be suitable. However, further inspection of (41) and (42)
1376 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990
reveals that, if an approximate guess for the stationary point trial phase is made as a pure component in the key component A, then (41) reduces to In YA = h A (52) This occurs because, for a pure component trial phase, equals one and the corresponding log of the activity coefficient equals zero. By use of (38), (52) reduces to
YA = zAyX
(53)
which implies that the number of moles in pure component trial A is equal to the activity of component A in the original mixture. Heidemann (1983) points out that an often overlooked result of inequality (34) is that a phase is unstable if a component's activity is greater than one. This corresponds to the observation of (53) since if component A has an activity greater than one, then the sum of the mole numbers, YA (since the trial phase is pure in A), is greater than one, and the mixture is detected as unstable according to (44). It then follows that one of the key components to initialize a search for the stationary point is the component with the largest activity. This is similar to Gautam and Seider (1979). Numerical experience has shown that a reliable starting point for the Newton-Raphson scheme is obtained by setting the mole number of the identified key component ( YA)equal to that component's activity in the original mixture, while all the remaining mole numbers are set close to zero. Furthermore, a good selection of a second starting point to attempt to find the other root of (41), if it exists, is to use the component with the second largest activity in a similar way. This choice, however, can cause some problems in certain instances, for example, a two-liquid-phase mixture involving small quantities of several hydrocarbons, just enough to form a second liquid phase, in a dilute ethanol-water solution. The activity of the hydrocarbon components will be large, because of the large, almost infinitely dilute, activity coefficients. Therefore, the search for the second root may find the same root as the first root search (which indicates instability), because the second largest activity identifies a hydrocarbon and not water. In this case it would be necessary to continue the root search from each component in decreasing the order of activity. The procedure used in this work to identify instability can therefore be summarized as follows: 1. Calculate the activities for all the components in the original mixture. 2. Select the two largest activities and identify the components. 3. Set the initial estimates of the stationary point trial phase mole numbers equal to zero except for the components identified in 2, which are set equal to their respective activities in the original mixture. 4. Start Newton-Raphson iterations from each initial estimate. If the value of r given by (47) exceeds 0.8, assume the root is converging on the trivial solution. 5. If both roots converge on the trivial solution or correspond to CblYi C 1, consider the original mixture stable. 6. If one or both of the roots are found where xF;lYi > 1,consider the mixture unstable. Normalize converged roots using (42) to provide starting points for the liquidliquid phase-split calculation. Numerical experience has shown that the above procedure is reliable for most mixtures and allows the Newton-Raphson iterations to converge in normally less than 10 iterations. Furthermore, a good estimate of the phase
fraction, s, can be obtained by using the key component identified with the largest activity, say A, in the following way: (54)
where xf\ and xf: are the mole fractions of A in the initial phase estimates at the stationary points calculated above.
Phase-Split Calculation Given the initial estimates of the phase split from the stability test, (31)-(33) can then be solved to provide the true equilibrium values. The method used in the work was to solve all the equations simultaneously by using the Newton-Raphson method. As with the stability test, the equations are written as deviation functions: (55) Di = xfyf - xflyfl i = 1, ..., n,
Di+", = zi - sxf - (1- s)xfl
i = 1, ...,n,
(56) (57)
The variable vector is given as XT
= ( x ~..., , X ~ ~ , X..., ~xI ~ , , s )
(58)
which is updated according to
+
= x ( ~ ) Ax
t = iteration counter
The correction vector Ax is calculated by using JAx = -D
(59) (60)
which is solved by standard Gaussian elimination. The iterations are continued until the mean square error condition for the residuals (D) is satisfied for a tolerance of 10-6.
Experience has shown that the step length corrections
(Ax)calculated by (60) can be too large, thus forcing the iterations to diverge or converge on the trivial solution. To avoid this, the maximum corrections are limited to xmax for the mole fractions and smax for the phase fraction. Values used in this work were typically xmax = 0.10
(61)
smax = 0.15
(62)
Using this scheme and the initial conditions of the stability test usually allow the liquid-phase split to be calculated in less than six iterations. It is possible that one of the two stationary points located for an identified unstable phase could be detected as converging on the trivial solution. This usually occurs for mixtures that are close to the binodal curve and, hence, close to one of the equilibrium phases (and stationary points). One option would be to increase the limit for the trivial solution indicator, r, from 0.8 to near 1 in the hope that the stationary point root is detected. Numerical experience, however, has shown that a good estimate of the composition for this phase is given by the following: 1. Set the mole fraction for the key component whose root search had been detected as approaching the trivial solution equal to one. 2. Set the mole fraction of the key component whose root search converged to a negative stationary point (and therefore detected an unstable phase) equal to zero. 3. Set the remaining component mole fractions equal to their values in the original mixture.
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1377 This estimate, along with the other stationary point, has been found to provide the Newton-Raphson phase split calculation with a good opportunity of avoiding the trivial solution, and converging on the true solution that usually has a very small (or near unity) value of s. The initial estimate of s can either be conservative at 0.5 or closer to 0 or 1 depending on which of the key component trial phases approached the trivial stationary point. This scheme has proved reliable even for very small values of s. However, if the trivial solution is approached by the phase split calculations, other options that can be tried are to reset the offending phase equal to the feed composition (setting s close to 0) or to begin from a pure component phase in the identified key component (setting s equal to 0.5).
New Simulation Algorithm The simulation method used in this work is based upon the Naphtali and Sandholm (1971) algorithm as described by Furzer (1986) and Cairns and Furzer (1988). Basing the iterations on the average liquid mole fraction, eq 14, then for constant molal overflow, we can write the deviation function for stage n as
derivatives of the deviation functions with respect to variables of composition and temperature. The elements of the B matrices correspond to the partial derivatives of the stage deviation functions with respect to the variables for that stage. The C matrix elements are derivatives with respect to the variables of the stage below, and the A matrix are elements with respect to the variables of the stage above. For example, the B matrix for a general stage n is given by
B, =
n = 1, ..., N dimension = (n, + l ) ( n ,
For the top stage (1) and the reboiler ( N ) ,we have
+ 1)
The individual elements of block A are described by Furzer (1986) and simply represent the ratio of the liquid flow rate leaving the stage above to the flow rate of liquid leaving the stage of interest. The B and C matrix derivatives are obtained by differentiating (63)-(68) with respect to the average liquid mole fractions and stage temperatures. For example, the B matrix derivatives for the general stage n, other than the top and reboiler stages, are = 1,...,n,; j = 1, ..., n,
i
aFn,i aYn,i = - -Vn ax+ L,
= 1, ..., n, (67)
1
i = 1, ..., n,
(73) (74)
while the C matrix derivatives are simply The Naphtali-Sandholm method groups these equations by stage and uses the Newton-Raphson iteration scheme xk+l = xk Ax k = iteration counter (69)
+
to update the variables housed in the matrix x. The correction matrix Ax is obtained from JAx = -F (70) where J is the Jacobian matrix, holding the deviation function partial derivatives, and F is the matrix of deviation function residuals. The Naphtali-Sandholm method allows the Jacobian to be set up in block tridiagonal form, and Furzer (1986) describes a technique that exploits this structure to calculate Ax calculated from (70) by simple Gaussian elimination. The Jacobian can be represented by rB,
c,
1
where A, B, and C are submatrices that house the partial
aFn,i
-=--
afn+l,j
V n + l dYn+l,i Ln
.
z = 1 , ...,n,; j = 1 ,..., n,
afn+l,j
(76)
At this stage, the advantage of using the average mole fractions as variables becomes clear. It allows the problem to be formulated in the same manner as the single-liquid-phase case with the important block tridiagonal structure of (71) being preserved. The execution of the algorithm by a computer program can be summarized by the flow sheets of Figures 5 and 6. The main program begins by reading the input data, which includes the necessary thermodynamic model parameters. A thermodynamic routine is then called to evaluate the appropriate vapor compositions and derivatives for each stage. For each stage the stability routine is called to assess if phase splitting will occur. If the liquid is considered stable, the standard two-phase y values and y derivatives
1378 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990
where to1 = 1.0 x 10-it"'
The integer value of itol is supplied by the input data file and typically has a value of six. If the convergence criterion has not been met, the B and C matrices are calculated by using (73)-(77) and the block elimination technique of Furzer (1986) is used to reduce (70) to a form where the correction matrix can be obtained by back-substitution. Both the block elimination and back-substitution steps require the use of a Gaussian elimination subroutine. The iteration variables are then updated according to (69); however, maximum projection limits are applied to the corrections. The limits are supplied by the input data file as ATMAx and AXMAX,and if (80) lax1 > AXMAX
7' caicv1ate deuiat1on function$
I
Calculate residue
I
write intermadiate print out i t rewired
I NO
write convergence message
lATl > ATMAX
Gaussian
Figure 5. Main program flow sheet.
0 START
i
f
(81)
for any interaction, the projections are reset to the maximum values. Furthermore, if the projections of the Newton-Raphson scheme produce negative mole fractions, these are reset close to zero. Cairns and Furzer (1988) have shown that this technique provides a robust and flexible way of constraining the Newton-Raphson iterations. Calculation of Y Values and Partial Derivatives. The liquid-liquid routine returns with a full thermodynamic description of both liquid phases, which includes the compositions and phase fractions, as well as the phase activity coefficients and their partial derivatives with respect to composition and equilibrium temperature. This information is necessary to enable the calculation of the partial derivatives of the vapor composition with respect to the overall liquid composition and equilibrium temperature. The calculation of the vapor composition in equilibrium with two liquid phases, themselves at equilibrium, can be based on either of the two liquid-phase K values. Choosing, arbitarily, phase 1, gives yn,i = Kl,n,i~l,n,i i = 1, ..., n,; n = 1, ..., N (82)
Calculate B and C
PerCom block e l m n a t i o n w i n g G a u ~ s i a nroutine
(79)
where
Calculare K-values LLE 31ven X P i t i a l
stable
C a l c u l a t e y-values
e s t i m a t e s from atability matine
Calculate y-derivatives
Equation 83 can be substituted into (82) and the result differentiated with respect to the overall liquid composition to give dYn,i
- Pi,i
d(Yl,n,iXl,n,i)
axnBi P, Calculate K-values
= 1, ..., 12,;
n = 1, ..., N, j = 1, ..., n, (84)
It is important to note that both the phase activity coefficients and mole fractions are functions of the overall liquid composition. Equation 84 can be expanded by using the product rule to give
routine
Calculate y - v a l u e s Ca1cula:e
afn3j
.
2
y-derivatives
Figure 6. Thermodynamic routine flow sheet.
are calculated. Otherwise, the liquid-liquid equilibrium routine is called to solve the isoactivity criterion given the initial estimates provided by the stability routine. The deviation functions are then calculated, and the convergence test given by Furzer (1986) is used to check that the residue of the deviation functions is less than the acceptable error tolerance. Convergence is satisfied if N
(78)
Xl,n,i
d f n ,j
+
dx 1,n,i P Y 1 , n . i afn, j
]
(85)
The phase 1 activity coefficient is in turn dependent on the phase composition and temperature: Yl,n,i
=
Yl,n,i(Xl,n,i,Xl,n,2,...,Xl,n,n,,Tn)
i = 1, ..., n,
(86)
Therefore, the activity coefficient derivative can be expanded to give
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1379 by forming the matrix equation
Ax = B
The derivatives of the phase activity coefficient with respect to phase composition are the standard derivatives obtained from the thermodynamic model and are available from the liquid-liquid routine. The problem is thus reduced to obtaining the derivative of the phase compositions with respect to the overall liquid composition, viz. ax,,n,i/aan,j i = 1, ..., n,; j = I, ..., n, (88) To calculate these values, it is necessary to make use of the equations that have been solved by the liquid-liquid routine. Recalling (31)-(33) and differentiating (31) with respect to fn, gives a(Xl,n,iYl,n,i) afn, j
-
a(x2,nj~2,n,i) afn, j
=0
(95) in which the elements of A are the coefficients of the left-hand side of the equation set, the elements of B are the values of the right-hand side of the equation set, and the elements of x are the unknowns. Since all the elements of A and B can be calculated, x can be found by standard Gaussian elimination. This is illustrated in the matrix equation shown in Figure 7 for a ternary system. Once the derivatives of the phase 1compositions with respect to the overall composition have been calculated, they can be substituted into (87) to give the activity coefficient derivatives. Both of these sets of derivatives are in turn used in (85) to evaluate the y value derivatives. Temperature Derivatives. Differentiating (82) with respect to stage temperature gives
i = 1, ..., n,; j = 1, ..., n, (89)
Expanding this equation by using the product rule and substituting (87) gives, upon rearrangement,
a5ln.i
a(Kl,n,ixl,n,i)
aT n
aTn
-=
(96)
and substituting (83) gives
is a function of temperature, The vapor pressure which is calculated by the Antoine equation In
PB,,i
= A, -
+T A2
A3
Expanding (97) gives j = 1, ..., n, (90)
This set represents (n, X n,) equations in 2(n, X n,) unknowns. The unknowns are the partial derivatives of the phase compositions with respect to the overall compositions, viz.
axl,n,i azn,j
a~2,n,i
-
i = 1, ..., n,; j = 1, ..., n, (91)
The remaining equations are obtained by differentiating (32) and (33). It is important to note, however, that the phase fraction (s,) is also a function of the overall composition and temperature. Therefore, extra n, unknowns have been added, namely, dsn/aan,j j = 1, ..., n, (92) This makes the total number of unknowns equal to n,(2n, + 1). This is matched by the (n,X n,) equations of (90) and the n,(n, + 1) equations generated when (32) and (33) are differentiated with respect to the overall liquid-phase mole fractions. The differentiation of (32) with respect to f,,,,yields
-aYn.i dT,
Xl,n.iYl,n,i a P k , i --P, aT"
+
The partial derivative of P i , i is easily obtained from (98) as
and the phase 1activity coefficient derivatives are calculated from the thermodynamic model and are returned from the liquid-liquid routine. Therefore, the only unknowns in (99) are the partial derivatives of phase 1composition with respect to temperature. These can be found in a similar manner to the compositional derivatives by differentiating (31)-(33) with respect to temperature to produce 2n, + 1 equations in 2n, + 1 unknowns. The unknowns are
An example of the matrix equation for a ternary system is given in Figure 8. The summation equation (34) when differentiated simply gives
Therefore, in summary, (901, (93), and (94) represent n,(2n, + 1)equations that can be solved for the nc(2nc+ 1)unknowns given by (91) and (92). This can be accomplished
Conclusions A new algorithm has been developed for multicomponent three-phase azeotropic distillation. The method uses an extended Naphtali-Sandholm method that accounts for liquid-phase stability testing and liquid-phase splitting. The equations that describe the system have been set with f , the liquid-phase mole fractions, as the primary variables. This change in the identification of the variables leads to important and significant advances in the structure of
1380 Ind. Eng. Chem. Res., Vol. 29, No. 7 , 1990
c
0
0
1 0
0
0
0
-
-
0
c
0
o
0
I
,
II
I
1
-e:
*"
"c.
0
0
0
x
x
0
0
I
0
0
0
0
P
-
O
0
c? 3
x
C
t
".
-
C
v;
N Y
0 .
I
I
v;
i
Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 1381 the equations, the formation of partial derivatives, the use of analytical partial derivatives from modern liquid-phase thermodynamic models, and improvements in the numerical procedures. The new algorithm has an excellent robustness characteristic and is more powerful in handling phase-stability analysis than the extensive list of published methods that have been reviewed. The predicted temperature and composition profiles along a multicomponent three-phase azeotropic distillation column have been compared with the first extensive set o{,experimental data given by Cairns and Furzer (1990) in part 1 of this series. Part 3 of this series will use the new algorithm to solve some of the examples discussed in this paper and look a t the issue of multiple solutions in azeotropic distillation.
Nomenclature a = activity f = fugacity, Pa D = distillate flow rate, kmol s-l F = feed flow rate or deviation variable, kmol s-l g = homotopy equation, or function h = liquid enthalpy, J kmol-’ H = vapor enthalpy or homotopy equation, J kmol-I J = Jacobian matrix K = K value L = liquid flow rate, kmol s-l N = number of theoretical stages n, = number of components P 8 = saturated vapor pressure, Pa P = total pressure, Pa
R = reflux ratio or gas constant s = phase fraction t = homotopy parameter T = temperature, K V = vapor flow rate, kmol x = liquid composition y = vapor composition Y = new modified vapour composition z = feed mixture composition y = liquid-phase activity coefficient Literature Cited Allgower, E.; Georg, K. Simplicial and Continuation Methods for Approximating Fixed Points and Solutions to Systems of Equations. SIAM Rev. 1980,22,28. Amundson, N. R.; Pontinen, A. J. Multicomponent Distillation Calculations on a Large Digital Computer. Ind. Eng. Chem. 1958, 50, 730. Arlt, W.; Grenzheuser, P.; Sorensen, J. M. Calculation of LiquidLiquid Equilibrium. German Chem. Eng. 1982,5,87. Baden, N. Comparison of Methods for Distillation Column Calculations. Ph.D. Thesis, Danmarks Tekniske Hojskole, 1984. Baden, N.; Michelsen, M. L. Computer Methods for Steady-State Simulation of Distillation Columns. Inst. Chem. Eng. Symp Ser. 1987,104,A425. Baker, L. E.; Pierce, A. C. Gibbs Energy Analysis of Phase Equilibria. SPE/DOE Second Joint Symposium on Enhanced Oil Recovery, Tulsa, OK, 1981. Baumgartner, A. M.; Swindells, R. J.; Fell, C. J. D.; Henry, R. D.; Wiley, D. E. Comparison of Design Packages for Ethanol Dehydration Column Design. 13th Australian Chemical Engineering Conference, 1985;Paper CAC. Block, U.; Hegner, B. Development and Application of a Simulation Model for Three-phase Distillation. AIChE J . 1976,22(31,532. Boston, J. F.; Britt, H. I. A Radically Different Formulation and Solution of the Sinele-Staee . Enp. - Flash Problem. C o m ~Chem. 1979,2, 109. Boston. J. F.: Shah. V. B. An Aleorithm for Rieoroue Distillation Calculations with Two Liquid-Phases. Pap; presented at the AIChE Meeting, Houston, TX, April 1979. Boston, J. F.; Sullivan, S. L. An Improved Algorithm for Solving Mass Balance Equations in Multistage Separation Processes. Can. J . Chem. Eng. 1972,50,663. I
Boston, J. F.; Sullivan, S. L. A New Class of Solution Methods for Multicomponent, Multistage, Separation Processes. Can. J. Chem. Eng. 1974,52,1. Bril’, Zh. A.; Mozzhukhin, A. S.; Petrova, T. V.; Petlyuk, F. B.; Serafimov, L. A. Computer Modelling of Liquid-Liquid-Vapour Equilibria in Multicomponent Mixtures. Part I: Liquid-Vapour Equilibrium in Binary Mixtures. Russ. J. Phys. Chem. 1973a,47 (lo),1466. Bril’, Zh. A.; Mozzhukhin, A. S.; Petrova, T. V.; Petlyuk, F. B.; Serafimov, L. A. Computer Modelling of Liquid-Liquid-Vapour Equilibria in Multicomponent Mixtures. Part 11: Liquid-Liquid-Vapour Equilibrium in Binary Mixtures. Russ. J. Phys. Chem. 1973b,47(lo), 1469. Bril’, Zh. A.; Mozzhukhin, A. S.; Petlyuk, F. B.; Serafimov, L. A. Computer Modelling of Liquid-Liquid-Vapour Equilibria in Multicomponent Mixtures. Part 111: Heterogeneous Liquid Equilibrium Diagrams for Ternary Mixtures. Russ. J . Phys. Chem. 1973c,47 (ll),1556. Bril’, Zh. A.; Mozzhukhin, A. S.; Petlyuk, F. B.; Serafimov, L. A. Computer Modelling of Liquid-Liquid-Vapour Equilibria in Multicomponent Mixtures. Part IV: Calculations on the Liquid-Liquid Equilibrium. Russ. J . Phys. Chem. 1973d,47 (ll), 1558. Bril’, Zh. A.; Mozzhukhin, A. S.; Petlyuk, F. B.; Serafimov, L. A. Computer Modelling of Liquid-Liquid-Vapour Equilibria in Multicomponent Mixtures. Part V Calculations of the Equilibrium in Multiphase Systems. Russ. J . Phys. Chem. 1974a,48 (5), 660. Bril’, Zh. A.; Mozzhukhin, A. S.; Petlyuk, F. B.; Serafimov, L. A. Mathematical Modelling of Multicomponent Mixture Rectification with Separation of the Liquid Phase on the Column Plates. Theo. Found. Chem. Eng. 1974b,8 (3),327. Bril’, Zh. A.; Mozzhukhin, A. S.; Petlyuk, F. B.; Serafimov, L. A. Mathematical Simulation and Investigation of the Heteroazeotropic Fractionation Process. Theo. Found. Chem. Eng. 1975,9 (6),761. Bril’, Zh. A.; Mozzhukhin, A. S.; Timofeev, V. S.; Serafimov, L. A. Computer Modelling of Liquid-Liquid-Vapour Equilibria in Multicomponent Mixtures. Part VI: Boundaries and Configuration of Bundles of Distillation and Rectification Lines for Heteroazeotropic Mixtures. Russ. J . Phys. Chem. 1976,50(3),349. Bril’, Zh. A.; Mozzhukhin, A. S.; Petlyuk, F. B.; Serafimov, L. A. Investigation of Optimal Conditions of Heteroazeotropic Rectification. Theo. Found. Chem. Eng. 1977, I 1 (6),675. Buzzi Ferraris, G.; Casapollo, A. Risoluzione di Sistemi di Equazioni Algebriche non Lineari e lor0 Applicazione a Problemi del’Ingegneria Chimica. Ind. Chim. It. 1972,8, 261. Buzzi Ferraris, G.; Morbidelli, M. Distillation Models for Two Partially Immiscible Liquids. AIChE J. 1981,27 (6),881. Buzzi Ferraris, G.; Morbidelli, M. An Approximate Mathematical Model for Three-phase Multistaged Separators. AIChE J . 1982, 28 (l),49. Cairns, B. P.; Furzer, I. A. Rapid-Robust Calculational Procedures for Multicomponent Distillation Problems. Chem. Eng. Aust. 1988,13 (2),13. Cairns, B. P.; Furzer, I. A. Multicomponent Three-phase Azeotropic Distillation. 1. Extensive Experimental Data and Simulation Results. Ind. Eng. Chem. Res. 1990,preceding paper in this issue. Dennis, J. E.; More, J. J. Quasi-Newton Methods, Motivation and Theory. SIAM Rev. 1977,19,46. Dluzniewski, J. H.; Adler, S. B. Calculation of Complex Reaction and/or Phase Equilibria Problems. Inst. Chem. Eng. Symp. Ser. 1972,5,4:21. Eubank, P. T.; Barrufet, M. A. A Simple Algorithm for Calculation of Phase Separation. Chem. Eng. Educ. 1988,22 (l),36. Fournier, R. L.; Boston, J. F. A Quasi-Newton Algorithm for Solving Multiphase Equilibrium Flash Problems. Chem. Eng. Commun. 1981,8 (4-6),305. Furzer, I. A. Distillation for Uniuersity Students; Furzer, I. A,, Publisher; Department of Chemical Engineering, University of Sydney: Sydney, Australia, 1986. Gautam, R.; Seider, W. D. Computation of Phase and Chemical Equilibria. Part I: Local and Constrain Minima in Gibbs Free Energy. Part 11: Phase Splitting. AIChE J. 1979,25 (6),991. Guffey, C. G.; Wehe, A. H. Calculation of Multicomponent LiquidLiquid Equilibrium with Renon’s and Black’s Activity Equations. AIChE J. 1972, 18 (5),913. Guo, T.-M. The Simulation of Three-phase Flash and Distillation Processes. Proc. World Congr. Chem. Eng. 3rd (Tokyo) 1986, Paper 121-301,769.
1382 Ind. Eng. Chem. Res., Vol. 29, No. 7, 1990 Hanson, D. N.; Duffin, J. H.; Somerville, G. F. Computation of Multistage Separation Processes; Reinhold: New York, 1962. Heidemann, R. A. Three-phase Equilibria Using Equations of State. AIChE J . 1974, 20 (51, 847. Heidemann, R. A. Computation of High Pressure Phase Equilibria. Fluid Phase Equilib. 1983, 14, 55. Heidemann, R. A.; Maudhane, J. M. Ternary Liquid-Liquid Equilibri: The van Laar Equation. Chem. Eng. Sci. 1975, 30, 425. Henley, E. J.; Rosen, E. M. Material and Energy Balance Computations; Wiley: New York, 1969. Kinoshita, M.; Hashimoto, I.; Takamutsu, T. A New Simulation Procedure for Multicomponent Distillation Column Processing Nonideal Solutions or Reactive Solutions. J . Chem. Eng. J p n . 1983a, 16 (51, 370. Kinoshita, M.; Hashimoto, I.; Takamutsu, T. A Simulation Procedure for Multicomponent Distillation Column within which Three Phases of Vapour and Two Partially Immiscible Liquids are Present. J . Chem. Eng. J p n . 1983b, 16 (61, 515. Kovach, J. W., 111; Seider, W. D. Heterogeneous Azeotropic Distillation-Homotopy-ContinuationMethods. Comput. Chem. Eng. 1987a, 11 (61, 593. Kovach, J. W., 111; Seider, W. D. Heterogeneous Azeotropic Distillation Experimental and Simulation Results. AIChE J. 1987b, 33 (81, 1300. Kovach, J. W., 111; Seider, W. D. Vapour-Liquid and Liquid-Liquid Equilibria for the System Sec-Butanol, Di-Secbutyl Ether, Water. J . Chem. Eng. Data 1988,33 (11, 16. Marina, J. M.; Tassios, D. P. Prediction of Ternary Liquid-Liquid Equilibrium from Binary Data. Znd. Eng. Chem. Proces Des. Deu. 1973, 12, 271. Maurer, G.; Prausnitz, J. M. Thermodynamics of Multicomponent Liquid-Liquid-Vapour Equilibria for Distillation Column Design. Inst. Chem. Eng. Symp Ser. 1979,56, 1.3/41. Michelsen, M. L. The Isothermal Flash Problem. Part I: Stability. Fluid Phase Equilib. 1982a, 9 (11, 1. Michelsen, M. L. The Isothermal Flash Problem. Part 11: PhaseSplit Calculation. Fluid Phase Equilib. 1982b, 9 (l), 21. Murray, W. Methods for Unconstrained Optimisation; Academic Press: New York, 1972. Myers, A. L.; Seider, W. D. Introduction to Chemical Engineering and Computer Calculations; Prentice-Hall: Englewood Cliffs, NJ, 1976. Naphtali, L. M.; Sandholm, D. P. Multicomponent Separation Calculations by Linearization. AIChE J . 1971 17 (11, 1148. Nelder, J. A.; Mead, R. A Simplex Method for Function Minimization. Comput. J . 1965, 7, 308. Nelson, P. A. Rapid Phase Determination in Multiple-Phase Flash Calculations. Comput. Chem. Eng. 1987, I 1 (6), 581. Niedzwiecki, J. L.; Springer, R. D.; Wolfe, R. G. Multicomponent Distillation in the Presence of Free Water. Chem. Eng. Prog. 1980, 76, 57. Null, H. R. Phase Equilibrium in Process Design; Wiley: New York, 1970. Ohanomah, M. 0.;Thompson, D. W. Computation of Multicomponent Phase Equilibria. Part I 1 Liquid-Liquid and Solid-Liquid Equilibria. Comput. Chem. Eng. 1984, 8 (3/4), 157-162. Pham, H. N.; Doherty, M. F. Design and Synthesis of Heterogeneous Azeotropic Distillation Sequences. 1. Heterogeneous Phase Diagrams. Chem. Eng. sci. 1990, in press. Pratt, H. R. C. Continuous Purification and Azeotropic Dehydration of Acetonitrile Produced by the Catalytic Acetic Acid-Ammonia Reaction. Trans. Inst. Chem. Eng. 1947, 25, 43.
Prausnitz, J. M.; Anderson, T.; Grens, E.; Eckert, C.; Hsieh, R.; OConnell, T. Computer Calculations for Multicomponent Vapour-liquid and liquid-Liquid Equilibria; Prentice-Hall: Englewood Cliffs, NJ, 1980. Prokopakis, G. J.; Seider, W. D.; Ross, B. A. Azeotropic Distillation Towers with Two Liquid Phases. Foundations of ComputerAided Chemical Process Design; Mah, R. s. H., Seider, W. D., Eds.; AIChE: Washington, DC, 1981; Vol. 2, p 239. Prokopakis, G. J.; Seider, W. D. Feasible Specifications in Azeotropic Distillation Towers. AIChE J. 1983, 29 (l),49. Pucci, A.; Mikitenko, P.; Asselineau, L. Three-phase Distillation. Simulation and Appfication to the Separation of Fermentation Products. Chem. Eng. Sci. 1986, 41 (3), 485. Ross, B. A. Simulation of Three-phase Distillation Towers. Ph.D. Dissertation, University of Pennsylvania, University Park, 1979. Ross, B. A.; Seider, W. D. Simulation of Three-phase Distillation Towers. Comput. Chem. Eng. 1981, 5 (11, 7. Russell, P. A. Chemical Process Simulation-The GMB System. Comput. Chem. Eng. 1980,4, 167. Schuil, J. A.; Bool, K. K. Three-phase Flash and Distillation. Comput. Chem. Eng. 1985, 9 (3), 295. Seydel, R.; Hlavacek, V. Role of Continuation in Engineering Analysis. Review Article No 24. Chem. Eng. Sci. 1987, 42 (6), 1281. Soares, L. de J. S.; Parker, L. G.; Ellis, S. R. M. Proc. Znt. Solvent Extract. Conf. 1971, 7. Soares, M. E.; Meoina, A. G.; McDermott, C.; Ashton, N. Three Phase Flash Calculations Using Free Energy Minimisation. Chem. Eng. Sci. 1982, 37 (41, 521. Stewart, W. E.; Levien, K. L.; Morari, M. Simulation of Fractionation by Orthogonal Collocation. Chem. Eng. Sci. 1985,40 (3), 409. Sumberova, V.; Hlavacek, V.; Holodnick, M. Calculation of LiquidLiquid Equilibria Problems by Continuous Analogy of NewtonKantorovich Method. Chem. Eng. Sci. 1977, 32, 783-785. Swank, D. J.; Mullins, J. C. Evaluation of Methods for Calculating Liquid-Liquid Phase-Splitting. Fluid Phase Equilibr. 1986,30, 101-110. Swartz, C. L. E.; Stewart, W. E. A Collocation Approach to Distillation Column Design. AIChE J . 1986, 32 (111, 1832. Swartz, C. L. E.; Stewart, W. E. Finite-element Steady State Simulation of Multiphase Distillation. AIChE J. 1987, 33 (121, 1977. Vickery, D. J.;Taylor, R. Path Following Approaches to the Solution of Multicomponent, Multistage Separation Problems. AZChE J . 1986, 32 (41, 547. Vickery, D. J.; Taylor, R. Multicomponent Separation Process Calculations-A Review. Technische Universiteit: Delft, 1987. Walraven, F. F. Y.; van Rompay, P. A Short-Cut Algorithm for Three-phase Distillation Towers. CHISA '87, Prague, 1987. Wang, J. C.; Henke, G. E. Tridiagonal Matrix for Distillation. Hydrocarbon Process. 1966, 45 (8), 155. Wayburn, T . L.; Seader, J. D. Solution of Systems of Interlinked Distillation Columns by Differential Homotopy-Continuation Methods. Foundations of Computer-Aided Chemical Process Design; Westerberg, A. W., Chen, H. H., Eds.; Wiley: New York, 1984; p 765. Wegstein, J. H. Commun. ACM 1955, I, 9. Whitehouse, P. A. Ph.D. Thesis in Chemical Engineering, University of Manchester Institute of Science and Technology, 1964. Receiued for reuiew April 25, 1989 Revised manuscript received October 4, 1989 Accepted October 27, 1989
