European Symposiumon ComputerAidedProcess Engineering- 11 R. Gani and S.B. Jorgensen(Editors) 9 2001 ElsevierScienceB.V. All rights reserved.
EQUISTAR:
153
R e l i a b l e S o f t w a r e for D e s i g n o f N o n i d e a l a n d R e a c t i v e S y s t e m s
S. T. Harding and C. A. Floudas I Department of Chemical Engineering, Princeton University Princeton, NJ 08544, USA
1
Introduction
Many commercial packages for process design and thermodynamic calculations are available today. Typically, these packages allow a broad variety of design problems to be solved and a large number of thermodynamic models from which to choose. Many design alternatives can be screened relatively quickly using such applications, and they greatly reduce the effort needed to design complicated processes. However, the solution methods employed by these packages frequently fail for systems that exhibit complex behavior. The reason for this failure is that local solution techniques are used to solve the nonlinear equilibrium equations that arise in the problem formulations. In this paper, a new tool for robustly and efficiently solving process design and thermodynamics problems is presented. This new tool is called EQUISTAR, which stands for E Q U i l i b r i u m Solution Toolkit for Azeotropic and Reactive Distillation Design. EQUISTAR addresses the need for a computational tool that can be reliably applied to highly nonideal and reactive systems and can solve the most commonly occurring thermodynamics problems. In section 2, the capabilities and design of EQUISTAR are presented. The global optimization algorithms that EQUISTAR employs are outlined in section 3.
2
O v e r v i e w of E Q U I S T A R
Local solution approaches have two features that make them attractive, especially for commercial implementations: they are fast, and they are relatively easy to implement. In order to achieve the kind of reliability that global optimization methods provide, one usually pays a premium in speed. In addition, global approaches are substantially more difficult to implement both in terms of the analysis of the mathematical equations, and in the actual programming effort. However, recent developments in global optimization approaches take a large step toward the practical implementation of global methods for process design. By analyzing the structure of the mathematical equations, it is possible to identify properties that can be exploited to greatly reduce the computational effort for guaranteed reliability. 2.1
EQUISTAR
Capabilities
EQUISTAR is a versatile package that incorporates global optimization methods in order to provide reliable solutions of many thermodynamic equilibrium problems that arise in the design and simulation of chemical processes. In addition, EQUISTAR can use this "toolkit" of robust equilibrium solution techniques to determine the design specifications for reactive and nonreactive distillation columns. The user is not forced to use the rigorous global optimization approach for any of the problems that EQUISTAR solves. In some cases, one may wish to do 1Author to whom all correspondence should be addressed.
154 a quick local search for a solution. EQUISTAR gives the user the option of specifying whether the problem is to be solved to global optimiality, or local optimality, and the user can specify the number of local searches that are performed. EQUISTAR can solve the following thermodynamic equilibrium problems: 1) phase and chemical equilibrium, 2) Gibbs free energy minimization, 3) phase stability (through Gibbs tangent plane distance minimization), 4) finding all homogeneous reactive and non-reactive azeotropes, 5) finding all heterogeneous reactive and non-reactive azeotropes, 6) isothermal or reactive flash calculation, 7) reactive or nonreactive bubble point calculation, and 8) reactive or nonreactive dew point calculation. The equilibrium conditions in each of these problems are only necessary conditions for the global equilibrium solution. Therefore, the solution that is obtained may correspond to a thermodynamically unstable system. EQUISTAR allows the user to specify whether or not to check the stability of each equilibrium solution that is obtained. The stability check is performed by solving the tangent plane distance minimization problem to global optimality, or until a negative tangent plane distance is located. In addition to solving stand-alone thermodynamics problems, EQUISTAR incorporates its solution algorithms into the capability of solving reactive and non-reactive distillation design problems. EQUISTAR provides the user a choice of algorithms for reactive or non-reactive distillation design algorithms, 1) a modification of the Inside-Out Algorithm, and 2) a modification of the Bubble-Point algorithm. EQUISTAR allows the user to choose from a wide range of thermodynamic models for representing the system's physical behavior. A number of equations of state axe available: 1) the Redlich-Kwong equation, 2) the Soave-modified Redlich-Kwong equation, 3) the PengRobinson equation, and 4) the van der Waals equation. In addition, several activity coefficient equations are available: 1) the Wilson equation, 2) the NRTL equation, 3) the UNIQUAC equation, and 4) the UNIFAC group-contribution method. 2.2
EQUISTAR
Software Design
The EQUISTAR program primarily consists of optimization problem formulations and highlevel optimization algorithms. The program is written in C. These formulations and algorithms are based on novel analysis by several authors and are described in section 3. EQUISTAR automatically generates problem formulations in the c~BB problem format. c~BB is a global optimization approach for solving general twice-continuously differentiable nonlinear programming problems developed by [2, 1]. c~BB is based on a branch-and-bound framework coupled with novel convex underestimators and the program manages the variable branching and the formulation of the upper and lower bounding problems. In order to solve the actual optimization problems, MINOPT is called. MINOPT is a Mixed-Integer Nonlinear OPTimization modeling language and solver developed by [10]. MINOPT converts the formulation of the optimization problem into a format that can be sent to an equation solver. MINOPT has interfaces with a number of linear, nonlinear, mixed-integer linear, mixed-integer nonlinear, and differential and algebraic equation solvers. Depending upon the type of problem that is passed, MINOPT converts the formulation into the correct format, sends it to the appropriate solver, and passes the solution back to the program that called it. Through this structure, the implementation of EQUISTAR and the formulation of its problems are independent of the local optimization method. The algorithms within EQUISTAR possess their own interdependencies, as shown in figure 1. The middle layer of the figure are the basic global optimization algorithms for solving
155
equilibrium problems. Note that each of these algorithms may call the phase stability problem to verify the stability of the solution. These algorithms call c~BB as the general global optimization solver. At the top of the figure are the highest level algorithms: distillation design and phase and chemical equilibrium. Each of these top level algorithms require the repeated solution of the equilibrium algorithms.
Figure 1: Relationship between EQUISTAR components
3
R e v i e w of S o l u t i o n M e t h o d s
Each of the problem types addressed by EQUISTAR has its own solution algorithm. section provides a summary of each of the algorithms. 3.1
This
Gibbs Free Energy Minimization
The minimization of the Gibbs free energy of the system is a fundamental approach for determining the equilibrium state of a system. A necessary and sufficient condition for equilibrium is that the Gibbs free energy of a system at constant temperature and pressure be at its global minimum. The general problem referred to as the Gibbs Free Energy Minimization Problem (GMIN) is defined as follows" Given N components with initial moles {sT1, sT,... ,n T } participating in up to P potential phases and R chemical reactions at constant temperature and pressure, find the mol vector n that minimizes the value of the Gibbs free energy function and satisfies the material balance constraints. The algorithm that EQUISTAR uses to determine the global minimum Gibbs free energy is based on the approach developed by [8] and is outlined below. 1. The user provides the system temperature and pressure and overall composition and the number and type of phases. 2. The GMIN problem is solved locally to generate an upper bound. 3. A branching variable is chosen and the current domain is partitioned by bisecting the bounds of the branching variable.
156
4. In each new domain a convex lower bounding problem is solved and the domain is discarded if the solution is greater than the current best upper bound. 5. Return to Step 2 and repeat until the best upper and best lower bounds converge. 6. The solution of the problem provides the composition of each phase.
3.2
Tangent Plane Distance Minimization
[3] and [11] have proved that a necessary and sufficient condition for a candidate equilibrium solution to be the true equilibrium solution is that the tangent plane distance function be nonnegative for all phases in the candidate solution. The tangent plane distance function is defined as the distance between the Gibbs free energy surface for the new phase, and the tangent plane to the Gibbs energy surface constructed at the candidate equilibrium solution. Based on the work of [9] and [4], the tangent plane distance minimization problem (TED) is solved in EQUISTAR using the following algorithm: 1. The user candidate 2. The T P D 3. Check the
provides the system temperature and pressure, and the composition of the phase. problem is solved locally to generate an upper bound. current best upper bound. I f it is less than zero, then stop the algorithm
because the candidate phase is unstable. 4. A branching variable is chosen and the current domain is partitioned by bisecting the bounds of the branching variable. 5. In each new domain a convex lower bounding problem is solved and the domain is discarded if the solution is greater than the current best upper bound. 6. Return to Step 2 and repeat until the best upper and best lower bounds converge or the best upper bound becomes negative. 7. The solution of the problem determines the stability or instability of the candidate phase. Enclosing All Azeotropes
3.3
The phenomenon of azeotropy occurs in many industrial applications. Azeotropes restrict the amount of separation of a multicomponent mixture that can be achieved by distillation. The ability to predict whether a given mixture will form one or more azeotropes and to calculate the conditions and compositions of each azeotrope is essential if one wants to model separation processes. An azeotrope is defined as a liquid mixture that boils at a constant temperature where the composition of the vapor phase is identical to the composition of the boiling liquid. When the boiling liquid contains a single phase, this phenomenon is called a homogeneous azeotrope. If the liquid consists of two or more phases it is classified as a heterogeneous azeotrope. Azeotropes may also occur in systems where one or more chemical reactions are occurring. These are called reactive azeotropes, and may be classified as homogeneous reactive azeotropes or heterogeneous reactive azeotropes, depending upon the number of liquid phases.
3.3.1
Enclosing All Homogeneous Azeotropes
The algorithm presented below for the location of all homogeneous non-reactive and reactive azeotropes is based on the work of [6]. 1. The user provides the system pressure.
157
2. A branching variable is chosen and the current domain is partitioned by bisecting the bounds of the branching variable. 3. In each new domain a convex lower bounding problem is solved and t h e d o m a i n is d i s c a r d e d if the solution is greater t h a n zero. 4. Return to Step 2 and repeat until all domains have been eliminated, or the size of the remaining domains are within a given tolerance. 5. The solution of the problem determines the temperature and composition of all homogeneous azeotropes in the system.
3.3.2
Enclosing All Heterogeneous Azeotropes
The algorithm presented below for the location of all heterogeneous non-reactive and reactive azeotropes is based on the work of [5]. 1. The user provides the system pressure. 2. A branching variable is chosen and the current domain is partitioned by bisecting the bounds of the branching variable. 3. I n e a c h n e w d o m a i n , t h e p o s s i b i l i t y of a t r i v i a l s o l u t i o n is c h e c k e d . If a t r i v i a l solution is possible, the d o m a i n is k e p t , b u t S t e p 4 is skipped for the d o m a i n . 4. In each new domain a convex lower bounding problem is solved and the domain is discarded if the solution is greater than zero. 5. Return to Step 2 and repeat until all domains have been eliminated, or the size of the remaining domains are within a given tolerance. 6. The solution of the problem determines the temperature and composition of all heterogeneous azeotropes in the system. 3.4
Flash Calculation
Unlike Gibbs free energy minimization, the flash calculation takes an equation-solving approach to the determination of phase and chemical equilibria. The solution of the isothermal flash and reactive flash problems in EQUISTAR provides all compositions that satisfy the flash equations. The algorithm is based on the approach for finding all solutions to nonlinear systems of equations developed by [7]. 1. The user provides the system temperature and pressure and feed rate and composition. 2. A branching variable is chosen and the current domain is partitioned by bisecting the bounds of the branching variable. 3. In each new domain a convex lower bounding problem is solved and the domain is discarded if the solution is greater than zero. 4. Return to Step 2 and repeat until all domains have been eliminated, or the size of the remaining domains are within a given tolerance. 5. The problem solution specifies the composition and flowrate of the vapor and liquid phases. 3.5
Bubble
Point and Dew Point Calculation
The calculation of bubble point temperatures and dew point temperatures is a natural extension of the flash calculation. These calculations are commonly encountered in the design and simulation of distillation columns. The bubble and dew point calculations are phase a n d / o r chemical equilibrium calculations.
158 In these formulations, the phase equilibrium condition is represented as the equality of chemical potentials for all components. The difference between the bubble and dew point problems and the flash problem is that the composition of only one phase is being determined, and the temperature of the equilibrium state is being determined. 1. The user provides the system pressure and liquid (vapor) composition. 2. A branching variable is chosen and the current domain is partitioned by bisecting the bounds of the branching variable. 3. In each new domain a convex lower bounding problem is solved and the domain is discarded if the solution is greater than zero. 4. Return to Step 2 and repeat until all domains have been eliminated, or the size of the remaining domains are within a given tolerance. 5. The solution of the problem determines the bubble (dew) temperature and the composition of the vapor (liquid) phase.
4
Conclusion
Based on significant developments in global optimization approaches developed over the past several years, EQUISTAR provides a suite of algorithms for the reliable and efficient solution of process design and thermodynamic equilibrium problems.
References [1] Adjiman C.S., Androulakis I.P., and Floudas C.A., 1998b, A global optimization method, aBB, for general twice-differentiable N L P s - II. Implementation and computational results. Comput. Chem. Eng. 22, 1159-1179. [2] Adjiman C.S., Dallwig S., Floudas C.A., and Neumaier A., 1998a, A global optimization method, aBB, for general twice--differentiable N L P s - I. Theoretical advances. Comput. Chem. Eng. 22, 1137-1158. [3] Baker L., Pierce A., and Luks K., 1982, Gibbs energy analysis of phase equlibria. Soc. Petro. Eng. J. (p. 731). [4] Harding S. and Floudas C., 2000a, Phase stability with cubic equations of state: A global optimization approach. AIChE J. 46, 1422-1440. [5] Harding S. and Floudas C., 2000b, Locating all heterogeneous and reactive azeotropes in multicomponent systems. I~EC Res. 39, 1576-1595. [6] Harding S.T., Maranas C.D., McDonald C.M., and Floudas C.A., 1997, Locating all azeotropes in homogeneous azeotropic systems. ISJEC Res. 36, 160-178. [7] Maranas C.D. and Floudas C.A., 1995, Finding all solutions of nonlinearly constrained systems of equations. Journal of Global Optimization 7, 153-182. [8] McDonald C. and Floudas C., 1994, Decomposition based and branch and bound global optimization approaches for the phase equilibrium problem. Journal of Global Optimization 5, 205-251. [9] McDonald C. and Floudas C., 1995a, Global optimization for the phase stability problem. AIChE J. 41, 1798. [10] Schweiger C.A. and Floudas C.A., 1998c, MINOPT A Modeling Language and Algorithmic Framework for Linear, Mixed-Integer, Nonlinear, Dynamic, and Mixed-Integer Nonlinear Optimization. Kluwer Academic Publishers, in preparation. [11] Smith J., Missen R., and Smith W., 1993, General optimality criteria for multiphase multireaction chemical equilibrium. AIChE J. 39, 707.