Устранение расхождений при расчете точек росы и кипения в многокомпонентных смесях

A Complementarity-based Vapor-Liquid Equilibrium Formulation for Equation-Oriented Simulation and Optimization Vibhav Dabadghao1, Jaffer Ghouse 2,3, John Eslick2,3, Andrew Lee 2,3, Anthony Burgard2, David Miller2, and Lorenz Biegler 1 1Carnegie Mellon University, Pittsburgh, PA 2National Energy Technology Laboratory, Pittsburgh, PA 3NETL Support Contractor, Pittsburgh, PA Thursday 1st December, 2022 Abstract Vapor-Liquid Equilibrium (VLE) is a cornerstone of computer aided process en- gineering (CAPE). Embedded within process system models, VLE calculations are inherently procedural with nonsmooth behavior that frequently require discrete deci- sions. Traditionally, these features resist the incorporation of VLE within efficient, large-scale equation-oriented (EO) process simulation and optimization strategies. On the other hand, recent reformulation of VLE models through the incorporation of complementarity constraints has broadened its scope to deal seamlessly with phase transitions and even supercritical excursions in process simulation and optimization. In this study, we extend these VLE complementarity models to EO frameworks where procedural thermodynamic property libraries are still required. Here, we develop an efficient, non-intrusive, and intuitive “square-flash” equation system that has been im- plemented within the IDAES Integrated Platform (IDAES-IP). The effectiveness of this modular approach is demonstrated on case studies for non-ideal flash calculations and distillation optimization with disappearing phases and supercritical transitions. 1 Introduction Vapor-liquid equilibrium (VLE) calculations are an important part of most process applica- tions. They especially appear in large-scale processes with heat exchanger networks, distilla- tion columns, emulsification and natural gas network pipelines, and in dynamic simulations and equation-oriented (EO) flowsheet optimization. Given their wide-ranging applications, the importance of performing VLE (or flash) calculations reliably and efficiently cannot be overstated. Flash calculations based on non-ideal thermodynamic relationships are well- known to be challenging to incorporate and solve in EO flowsheet optimization. These 1 algorithms experience difficulties with VLE systems, especially in operating conditions cor- responding to the existence of a single phase, since the flash equations are only valid in the two-phase region between the bubble point and the dew point. The literature on formulation and calculation of single-stage flash can be broadly classified into two categories: formulations relying on sequential modular (SM) approaches involving iterative solution procedures and those based on EO approaches involving simultaneous so- lution of the equation system. Rachford and Rice [1] proposed an iterative algorithm in which the vapor fraction is determined using a recursive binary search over its entire domain of values from zero to one. Subsequently, several formulations based on the Rachford-Rice reformulation were proposed that relied on nested-loop type algorithms. These formulations are computationally expensive, especially when phase equilibria and enthalpy departures are defined by complex models. Boston and Britt [2] introduced a radically different class of algorithms known as inside-out algorithms. These algorithms are based on an itera- tive nested-loop strategy in which the calculation of the flash equations is separated from the computation of thermodynamic properties. A key idea used in these algorithms is the introduction of nonphysical variables that have low sensitivity to primitive variables such as temperature, pressure, and phase fractions. Because of this low dependence, their per- formance is largely insensitive to the initial guess of the iterative procedure. In practice, however, it is understood that a substantial amount of initialization is required to begin the flash calculations. Compared to earlier works, inside-out algorithms are computationally more tractable, broadly applicable, and agnostic to the thermodynamic models involved as well as whether the mixture is wide-, moderate- or narrow-boiling. Inside-out algorithms and their modifications have been widely used in commercial software as they are easy to construct and computationally inexpensive. However, these algorithms are only applicable to situations in which the solution of the flash calculation lies in the two-phase vapor-liquid region. Single-phase solutions often occur during process simulation or optimization, where the converged solution indicates an all-liquid or an all-vapor result. Instead, heuristic meth- ods are typically used to find the correct single-phase result in such situations, and these add to the computational complexity and can lead to numerical failure. Watson et al. [3] proposed a modification of Boston and Britt’s inside-out algorithm in which the existence of single-phase results is handled systematically using a nonsmooth reformulation of the flash problem. The modified formulation involves the addition of a nonsmooth equation in the in- ner loop to relax equilibrium relations when the operating conditions indicate a single-phase result. The outer loop is solved using fixed-point iterations, while the solution of the inner loop requires nonsmooth equation solvers such as the semismooth Newton method [4] or the LP-Newton method [5]. In SM approaches, modules are solved sequentially using black-box models, and each module can be solved using custom numerical strategies, as is done in [3]. Modules in a flowsheet 2 are relatively simple to build and debug because of customized procedures. However, SM approaches suffer from the drawback of having to solve modules repeatedly in a single so- lution procedure because the intermediate failure of a module can be fatal for a controller during online optimization. Furthermore, derivatives of the implicit models are not directly available and have to be estimated using finite difference approximations leading to inac- curacies and round-off errors. In EO approaches, the model equations corresponding to either process units or entire flowsheets are solved simultaneously. The accuracy and overall performance of optimization algorithms are not subject to rounding errors in gradients due to the availability of exact derivatives. This approach requires precise problem definition, and the significantly larger problem size compared to the SM approach necessitates careful initialization and a robust nonlinear equation solver. As such, compared to the inside-out algorithms, EO algorithms for the flash problem do not require either heuristics or special- ized solution procedures to handle single-phase results. They can be systematically used to leverage properties of cubic equations of state but can be readily extended to incorporate any thermodynamic model. The EO approach is especially preferred for complex flowsheets with many nested recycle loops and implicit design specifications [6]. For the single-stage flash problem, fully EO formulations have been developed that auto- matically handle direct transitions to the correct number of phases. These approaches are based on the necessary conditions for Gibbs free energy minimization of the flash prob- lem, with complementarity conditions that reflect the number of phases. Initially, the flash formulation was developed through relaxation of the pressure specification and LP-Newton methods in [7] and then extended to smoothed optimality conditions in [8]. Both approaches allow the simulation of distillation columns with dry and vaporless trays. Also, these for- mulations were extended to the optimal design of distillation systems in [9]. Phase changes in VLE systems can also be handled using mathematical programs with complementarity constraints (MPCCs1). Complementarities are relationships between variables where either one (or both) must be at its bound [10]. These constraints are often embedded in opti- mization problems to handle discrete decisions with an EO approach, as SM approaches have significant limitations in process optimization [11, 12]. Formal MPCC formulations were adopted for the flash problem in [13, 14] where several MPCC algorithms were demon- strated and compared for optimization of distillation columns. Dynamic optimization of distillation columns with phase transitions using a simultaneous discretization and solution strategy for the MPCC problem was described in [15]. Extending these results, Kamath et al. [6] proposed an optimization-based EO formulation with an emphasis on thermodynamic models based on cubic equations of state. In this work, the problem of handling single-phase results is reformulated using complementarity constraints. Moreover, the properties of the roots of the cubic equation of state and its first and second derivatives are leveraged to 1MPCC and the related mathematical programs with equilibrium constraints (MPEC) are often used interchangeably in the literature. 3 handle disappearing phases by embedding these properties in the EO formulation to identify and isolate the appropriate roots. A correction factor is added to the equilibrium constraint to relax the equation when single-phase results are encountered. The corrector is tied to slack variables which, in turn, are tied to vapor and liquid flowrates via complementarity constraints. These complementarities are reformulated using penalty terms in an objective function, and the formulation is recast as an optimization problem. Similar optimization- based EO formulations have been extended to support optimization of complex flowsheets [16, 17]. Alternatively, Burgard et al. [18] proposed a smooth, square formulation for single-stage flash that does not require solving an optimization problem or adding penalty terms to an objective function. In this formulation, artificial variables representing mole fractions in “shadow” streams are used to calculate the bubble and dew points for the mixture. An equilibrium temperature is introduced and restricted to the two-phase region. Equilibrium relations, which are invalid in the single-phase region, are always calculated at this temper- ature to avoid failures due to non-convergence in the case of single-phase results. Moreover, smoothing methods for complementarity constraints are used to ensure the problem remains differentiable throughout its domain. The formulation is also extended to incorporate super- critical conditions via complementarity constraints to ensure that vapor flows vanish above the critical pressure. In this work, we present an improved EO formulation for the single-stage flash as well as extensions to distillation models. This formulation is a smooth, square model that elimi- nates the calculation of separate bubble and dew points to define the two-phase region and significantly reduces the computational cost. Instead, a “pseudo” equilibrium temperature is defined using complementarity constraints with the vapor and liquid flowrates. The com- plementarities ensure that the appropriate slack variables and/or phase flowrates become zero. Furthermore, in the case of single-phase results, nonzero slack variables ensure that the VLE relations are always evaluated in the two-phase region. Moreover, the formulation identifies the existence of the appropriate phase as dictated by the operating conditions us- ing derivatives of the cubic equation of state (CEOS). The properties of the derivatives at the roots are embedded in the formulation directly. Moreover, the extension to supercritical conditions, based on complementarity constraints, is much simpler than in previous works. The derivation of the proposed square-flash formulation is presented in detail in Section 2, along with new complementarity formulations for cubic equations of state and supercritical extensions. Section 3 discusses the implementation of this formulation within the Institute for the Design of Advanced Energy Systems Integrated Platform (IDAES-IP) [19]. Here, attention is focused on the use of external functions for evaluating the roots of the CEOS as well as unit initialization strategies within IDAES-IP. Section 4 presents results of case studies for two VLE systems under a broad set of specifications, as well as distillation optimization 4 with disappearing phases. Section 5 summarizes the paper and suggests areas for future work. 2 Flash equations Consider the general flash problem where we assume the standard form for calculating the liquid and vapor phase flowrates and compositions for a mixture of components in set C as given in equation (1). F =V +L (1.1) Fzi =Lxi +Vyi, i ∈C (1.2) 0 = X i∈C (yi−xi) (1.3) yi =Ki (T,P, x, y)xi, i ∈C (1.4) FhF =VhV +LhL−Q (1.5) where F , L, and V are the flowrates of the inlet, outlet liquid, and outlet vapor streams, respectively, while z, x, and y are their mole fractions. The equilibrium ratios, Ki, are defined by an equation of state and calculated by external functions with arguments of temperature, pressure, and the liquid and vapor phase compositions. Similarly, hF , hL, and hV are the enthalpies of the inlet, outlet liquid, and outlet vapor streams, respectively, while Q is the heat duty. As such, formulation (1) implicitly assumes that the solution of the flash problem lies in the two-phase vapor-liquid region. This is because equation (1.4) is only valid in the two-phase region, and the formulation will lead to numerical failure if the operating conditions indicate a single-phase result. This often requires careful treatment using heuristics or other specialized solution procedures to converge to the correct result. As a result, the challenges of developing and implementing a VLE formulation that is complete in its applicability as well as tractable for its numerical solution can be summarized as follows.

• The formulation must be valid for both single-phase and two-phase results.
• The formulation should lead to a square set of equations that accommodates a modular framework.
• It should account for the entire range of operating conditions, including supercritical conditions.
• The appearance or disappearance of phases must be handled in a way that is amenable for numerical solution algorithms. We subsequently present various aspects of the proposed formulation that help address these challenges. 5 2.1 Handling single- and two-phase conditions In this study, we specifically focus on VLE with CEOS. We leverage properties of the CEOS roots, which take the form: f(Z) =Z 3 +bZ 2 +cZ +d = 0 (2) Kamath et al. [6] describe the conditions that must hold for both single-phase and two-phase conditions in the form of equality and inequality constraints. These can be summarized in constraint set (3). f(ZV ) = 0, f′(ZV )≥ 0, f′′(ZV )≥ 0 (3.1) f(ZL) = 0, f′(ZL)≥ 0, f′′(ZL)≤ 0 (3.2) where ZL and ZV are the liquid and vapor roots, and f′ and f′′ are the first and second derivatives of the cubic equation of state, respectively. Imposing the constraints in (3.1) and (3.2) along with the VLE model based on the cubic equation of state allows us to distinguish between the vapor and liquid root. 2.2 Improved formulation for flash calculations Traditionally, flash calculations have been performed using SM approaches. Here, condi- tional checks are applied to determine whether the system is in single-phase or two-phase by calculating the bubble and dew points and comparing them with the system state. The equilibrium expression is applied only if the system is in two-phase. On the other hand, it has been reported that the use of cubic equations of state with such iterative approaches can lead to false convergence, and several improvements have been suggested to overcome these chal- lenges [20–23]. However, these post-processing conditional checks introduce non-smoothness to the problem and are not suitable for large-scale flowsheet optimization. We propose an EO-based formulation that overcomes these drawbacks in a systematic way and can automatically handle different phase conditions without requiring specialized solu- tion procedures. In this work, we replace the specified flash temperature,T in equation (1.4), with an equilibrium pseudo temperature, T . yi =Ki

T,P, x, y
xi, i ∈C (4) The relation between T and T is given by the following complementarities: T−T =sV−sL (5.1) 0≤sV⊥V≥ 0 (5.2) 0≤sL⊥L≥ 0 (5.3) 6 Here, sV and sL are slack variables associated with the vapor and liquid stream flowrates, respectively. When the feed conditions imply an all-liquid result, V = 0, sV ≥ 0, L > 0 and sL = 0. Thus, from equation (5.1), T = T +sV will ensure that T is at the bubble temperature where equation (4) is valid. Similarly, for an all-vapor result, L = 0, sL≥ 0, V > 0 andsV = 0;T =T−sL will forceT to be at the dewpoint where equation (4) is valid. Finally, in the case of a two-phase result, both sV and sL vanish, and T = T as desired. Note that equation (5) also implies that 0 ≤sL⊥sV≥ 0. For cubic equations of state, the second derivative properties, equations (3.1) and (3.2), can be used to distinguish between the liquid and vapor roots. These properties are enforced using the following complementarities: f′′(ZV ) =g+ V−g− V (6.1) f′′(ZL) =g+ L−g− L (6.2) 0≤g− V⊥V≥ 0 (6.3) 0≤g+ L⊥L≥ 0 (6.4) Again, g+ V , g− V , g+ L , and g− L are non-negative slack variables associated with the vapor and liquid stream flowrates. For an all-liquid result, g+ L = 0 and f′′(ZL) =−g− L≤ 0. Similarly, for an all-vapor result, g− V = 0 and f′′(ZV ) = g+ V ≥ 0. Finally, for a two-phase result, g− V = g+ L = 0 and the derivative properties of both roots in equations (6.1) and (6.2) are satisfied as per equations (3.1) and (3.2), respectively. In the case of a single real root, a pseudo-root corresponding to the missing phase is determined from a translated cubic equation described in Section 3.1. 2.3 Extension to supercritical conditions The flash formulation can be further improved to automatically handle conditions in which a stream is above its critical point. This is important as we cannot know a priori whether the operating point belongs to the subcritical or the supercritical region, and systematic handling of these cases can avoid failures in flowsheet optimization. First, we estimate the critical temperature and pressure of the mixture. This can be done from the following equations, as described in further detail by Dowling and Biegler [16]. ΩA =am(z) Pc T 2 cR2 (7.1) ΩB =bm(z) Pc TcR (7.2) where ΩA and ΩB are constants specific to the cubic equation of state; am and bm are interaction parameters of the mixture and are dependent on the feed composition, z;Pc and Tc are the critical pressure and temperature of the mixture; and R is the gas constant. As 7 per Burgard et al. [18], we introduce artificial pressure variables P + and P− to determine whether the operating pressure is below or above the critical point. P + = 1 2  P−Pc + q (P−Pc)2 +ε2  ≈ max(0,P −Pc) (8.1) P− = 1 2  Pc−P + q (Pc−P )2 +ε2  ≈ max(0,Pc−P ) (8.2) Here,ε is a smoothing parameter. To ensure phase equilibrium is never evaluated above the critical point, a subcritical pseudo-pressure, P , is introduced and assigned as P =P−P +− ε 4 (9) and the equilibrium equation (4) is modified and replaced by equation (10). yi =Ki

T,P, x, y
xi, i ∈C (10) The complementarity (11) ensures that the vapor flowrate vanishes in the supercritical region, as the mixture is assumed to exhibit liquid-like behavior. 0≤P +⊥V≥ 0 (11) Finally, since equations (5.3) and (6.4) both complement L, they can be combined and replaced with equation (12). 0≤g+ L +sL⊥L≥ 0 (12) 2.4 Smooth approximations of complementarities The above set of equations describes vapor liquid equilibrium for any set of feed specifications, accounting for both subcritical and supercritical conditions. The non-smoothness of the complementarity constraints is handled by replacing them with smooth approximations, i.e., for scalars ξ,η : 0≤ξ⊥η≥ 0⇐⇒ξ− max(0,ξ−η) = 0≈ 1 2  ξ +η− p (ξ−η)2 +ε2  = 0 This property allows us to handle single-phase solutions in equation-oriented simulation and optimization frameworks. The proposed formulation for flash calculations is presented in equation system (13). F =L +V (13.1) Fzi =Lxi +Vyi, i ∈C (13.2) 8 FhF =LhL (x,T,P ) +VhV (y,T,P )−Q (13.3) yi =Ki

T,P, x, y
xi, i ∈C (13.4) 0 = X i∈C (yi−xi) (13.5) T−T =sV−sL (13.6) 0 = 1 2  sV +V− q (sV−V )2 +ε2  (13.7) f′′(ZL) =g+ L−g− L (13.8) f′′(ZV ) =g+ V−g− V (13.9) 0 = 1 2  g− V +V− q g− V−V
2 +ε2  (13.10) 0 = 1 2  g+ L +sL +L− q g+ L +sL−L
2 +ε2  (13.11) ΩA =am(z) Pc T 2 cR2 (13.12) ΩB =bm(z) Pc TcR (13.13) P + = 1 2  P−Pc + p (P−Pc)2 +ε2  (13.14) P− = 1 2  Pc−P + p (Pc−P )2 +ε2  (13.15) P =P−P +− ε 4 (13.16) 0 = 1 2  P + +V− q (P +−V )2 +ε2  (13.17) Formulation (13) is equivalent to the KKT conditions for free energy minimization of the VLE system. Further details of this derivation are shown in the Appendix A.4. This sys- tem of equations represents the square-flash formulation that handles phase transitions and extensions to the supercritical region. This set of equations can be incorporated in modular form within most CAPE platforms, as shown in Section 3 for IDAES-IP. Detailed discussion on the properties of the roots of the cubic equation of state is presented in appendices A.2 and A.3. 3 Integration into IDAES-IP This section describes the implementation of the proposed square-flash formulation within IDAES-IP and focuses on its integration with CEOS thermodynamic packages. In particular, Section 3.1 describes the calculation of modified cubic roots to deal with disappearing phases in the square-flash and how external functions are leveraged for this. Section 3.2 then discusses the initialization strategies that have been implemented for unit models with VLE. 9 3.1 External Functions to Compute Modified Cubic Roots In a purely equation-oriented framework, it may be difficult to solve for the three roots of a cubic and to correctly select the liquid and vapor (and the unused middle) root. To improve robustness, the analytical cubic formula can be used instead to directly calculate the roots of a cubic equation. However, this also requires the use of conditional statements to correctly handle the single and triple real root cases. In this study, external functions were written to calculate the liquid and vapor roots of a cubic directly using the analytic cubic formula. This function provides first and second derivatives and a pseudo-root when a real root does not exist. The procedure was implemented as an AMPL user-defined function [24]. Two external functions provide the vapor and liquid roots for equation (14) if they exist. The functions have the form Z(b,c,d ). The first and second derivatives of Z with respect to b, c, and d are calculated implicitly and provided to the solver. f(Z) = 0 =Z 3 +bZ 2 +cZ +d (14) When the real liquid or vapor root does not exist, a pseudo-root is provided to improve solution robustness. The turning points, i.e., the points at which f′(Z) = 0, are used to generate an extension of the cubic function and are given by equation (15). t1,t 2 =−2b∓ √ 4b2− 12c 6 (15) When the vapor root does not exist, the original cubic is translated so that the lower turning point of a translated cubic equation is assigned to the higher turning point of the original cubic, and the lower root of the translated cubic is defined as the pseudo vapor root. This is illustrated in Figure 1a, in which the vapor root of the original cubic equation does not exist. The new curve corresponding to the missing vapor phase is the translated cubic, and its corresponding root is the pseudo vapor root. This pseudo vapor root is the lower root of equation (16). fph(Z) = (Z +t1−t2)3 +b(Z +t1−t2)2 +c(Z +t1−t2) +d +f(t2)−f(t1) (16) Similarly, when the liquid root does not exist, the original cubic equation is translated so that the higher turning point of the translated cubic is assigned to the lower turning point of the original cubic, and the higher root of the translated cubic is used as the pseudo liquid root. This is shown in Figure 1b. In this case, the new curve corresponding to the missing liquid phase is the translated cubic, and the corresponding root is the pseudo liquid root. This root is the lower root of equation (17). fpl(Z) = (Z−t1 +t2)3 +b(Z−t1 +t2)2 +c(Z−t1 +t2) +d +f(t1)−f(t2) (17) 10 −2 −1 0 1 2 3 4 Z −4 −2 0 2 4 6 f (Z) Liquid Vapor Cubic Liquid root Vapor root (a) Vapor root does not exist −2 −1 0 1 2 3 4 Z −4 −2 0 2 4 6 f (Z) Liquid Vapor Cubic Liquid root Vapor root (b) Liquid root does not exist Figure 1. Pseudo roots are calculated for the missing phase using the translated cubic When the cubic has no turning points and only one root, the same root is returned for the low and high root. This is illustrated in Figures 2a and 2b. In these cases, the translated functions −3 −2 −1 0 1 2 3 Z −4 −2 0 2 4 6 f (Z) Liquid Vapor Cubic Liquid root (a) Vapor root does not exist −3 −2 −1 0 1 2 3 Z −4 −2 0 2 4 6 f (Z) Liquid Vapor Cubic Vapor root (b) Liquid root does not exist Figure 2. In the absence of turning points, translated curves coincide with the original cubic coincide with the original cubic function, and our pseudo-root scheme provides continuous functions that always return a real root. We note that these functions are not smooth at the point corresponding to f(Z) = f′(Z) = 0. Nevertheless, this pseudo-root approach still greatly improves solver performance, because the nonsmoothness affects neither the complementarity equations nor input-output relations in the VLE block. 3.2 Procedural Initialization in IDAES-IP Lee et al. [19] describe in detail the building blocks within the IDAES-IP and how the block hierarchy is leveraged to initialize an EO flowsheet. The block hierarchical nature enables 11 solving sub-models in isolation, as long as the sub-model equations that need to be initialized have zero degrees of freedom. This capability of activating and deactivating constraints or entire model blocks from the problem to solve in isolation is extensively leveraged for the set of highly nonlinear equations involved in unit models (or flowsheets) with flash calculations. Section 3.2.1 presents the custom initialization routine built within the state blocks in IDAES-IP, specifically, for thermophysical property calculations using the CEOS. Section 3.2.2 describes the initialization routines for the unit model blocks , specifically, the flash and distillation tray unit models that have VLE calculations. 3.2.1 Initialization of Thermophysical Properties The state blocks are implemented as self-contained sub-models representing a single material state, i.e., point in space and time, along with any equations for calculating the thermo- physical properties as a function of the state. The state is typically represented by state variables such as flow rate, temperature, pressure, and composition. Given an initial guess, each material state is fully defined and can be solved in isolation. Additionally, the phase equilibrium problem can be separated from the calculation of the other thermophysical prop- erties. This is analogous to the approach proposed by Boston and Britt [2] for inside-out type approaches. The initialization routine for CEOS is illustrated in Figure 3. Calculate 𝑇𝑇𝑐𝑐, 𝑃𝑃𝑐𝑐 Calculate 𝑇𝑇, 𝑠𝑠𝑉𝑉, 𝑠𝑠𝐿𝐿 Estimate 𝑇𝑇Bubble, 𝑇𝑇Dew Calculate 𝑃𝑃, 𝑃𝑃+, 𝑃𝑃− Calculate interaction parameters and 𝑓𝑓𝑓𝑓 at 𝑇𝑇, 𝑃𝑃 Calculate 𝑔𝑔𝑉𝑉 +, 𝑔𝑔𝐿𝐿 +, 𝑔𝑔𝑉𝑉 −, 𝑔𝑔𝐿𝐿 − Estimate 𝐿𝐿, 𝑉𝑉 Solve phase equilibrium subsystem Solve flash formulation Figure 3. Initialization routine for the thermophysical properties The critical point, which depends only on the mixture composition, is calculated by solving 12 equation set (7). Subsequently, the artificial pressure variables are calculated in equations (8) and (9). As discussed earlier, the calculation of bubble and dew points is obviated from the proposed formulation. However, during initialization, these are estimated based on ideal fluid thermodynamics so as to provide an initial guess for the equilibrium temperature, T . The equilibrium temperature is then initialized as T = (Tbubble +Tdew)/2, where Tbubble and Tdew are the estimated bubble and dew temperatures, respectively. The temperature slacks sV andsL are initialized using equation set (5). The estimated equilibrium temperature and pressure are used to calculate the interaction parameters of the cubic equation of state and, subsequently, its second derivativef′′. The cubic complementarity slacks are then initialized using equation set (6). The temperature slacks are used to determine whether the mixture is single-phase or two-phase. For a two-phase mixture, the phase flow rates are estimated using a linear interpolation of the phase fraction between the bubble and dew temperatures. As mentioned earlier, the phase equilibrium problem comprising equations (13.1)–(13.5) is separated from the other calculations and solved to refine initial guesses and converge the primitive equations of the flash problem. Using this refined initialization, equation set (13) is solved for the inlet state, and the process is repeated for the outlet state. 3.2.2 Initialization of Unit Operations Unit Model Blocks in the IDAES-IP are essentially composed of inlet and outlet state blocks that are linked through material, energy, and momentum balances. In addition, there are performance equations for the corresponding unit model. One of the first steps in initializing a unit model is initializing the state blocks and then moving on to initialize the conservation and performance equations. The initialization of unit models can be summarized as follows:

1. Initialize the inlet states at the initial guesses provided or at the inlet specification provided.
2. Determine initial guesses for intermediate and outlet states based on inlet states.
3. Initialize intermediate and outlet states at initial guesses.
4. Solve the material, energy, and momentum balances. This may require a number of incremental steps depending on the complexity of the unit operation. Two examples for initialization of unit operations that involve vapor-liquid equilibrium are discussed below. These are implemented for the case studies discussed in Section 4. 3.2.2.1 Example of Initializing a Flash Operation One of the most common process unit operations is the Flash drum, which consists of a single vessel with one inlet and two outlets. Within IDAES-IP, flash operations are represented as a unit model with a single inlet state block and a single outlet mixed equilibrium state block, which is then partitioned into vapor and liquid flows, and users specify the operating 13 state of the vessel by fixing two of the state variables (e.g., P,T ) for the mixed equilibrium state. When initializing the unit model, the inlet condition provided is used as an initial guess. Next, a set of initial guesses for the outlet state is calculated based on the inlet state provided. For the flash formulation, this means guesses for T,P of the equilibrium state are required. We can leverage some process knowledge to guess these conditions by recognizing that a well-designed flash operation will operate somewhere between the bubble and dew points in order to achieve a phase separation. However, bubble and dew points are difficult to calculate and are often not known a priori. As a fallback from this, we can recognize that, in order to minimize exergy losses, the outlet state should be as close as possible to the inlet state. Thus, it is a reasonable guess that the unknown state variable will be close to that of the inlet. With this information, we now have reasonable initial guesses for the outlet state and can initialize at these conditions. 3.2.2.2 Example of Initializing a Distillation Column Distillation columns are another common unit operation involving vapor-liquid equilibrium, and due to the highly coupled nature of the system, with recirculating vapor and liquid flows, have proven to be one of the most challenging unit operations to model and initialize robustly. Within IDAES-IP, the tray column unit model is represented as a system of interconnected sub-models as shown in Figure 4 and includes the following: C R F CondenserVin Distillate Reflux Vapor ReboilerLin Bottoms Vapor Reboil V D B Feed Tray Lin Lout Vout Vin F Tray 1 Lin Lout Vout Vin Tray n Lin Lout Vout Vin Inlet ports with individual state blocks Outlet ports with a single shared state block Figure 4. Block structure of the distillation tray unit model in IDAES

• A feed tray with three inlets (feed, liquid, and vapor) and two outlets (vapor and 14 liquid).
• A series of trays above the feed tray (rectification section) with two inlets (liquid and vapor) and two outlets (liquid and vapor).
• A series of trays below the feed tray (stripping section) with two inlets (liquid and vapor) and two outlets (liquid and vapor).
• An optional condenser, which can be either a partial or a total condenser, above the top tray that recycles a portion of the vapor from the top tray back as liquid.
• An optional reboiler below the bottom tray that recycles a portion of the liquid from the bottom tray back as vapor. Similar to the flash model, for the sub-models listed above, the inlet ports have an associated state block, but there is only a single state block for the outlet mixed-phase. The outlet liquid and vapor ports are then populated by creating references to the appropriate phase from this mixed-phase state. This structure allows solving for each tray, condenser, or reboiler independently (as long as the inlet conditions are fixed or an initial guess is provided) but at the cost of carrying repetitive state blocks when the state conditions are the same. For a column with a single feed, 10 trays, a condenser, and a reboiler, the total number of state blocks with VLE equations is equal to 35. Therefore, to successfully initialize the tray column, it is imperative to have an initialization routine for the state blocks that is fast and robust, as outlined in 3.2.1. The first step in initializing a distillation column model is to initialize the feed tray. In a well-designed column, the conditions on the feed tray should be very close to those of the feed stream and, hence, we use the conditions of the feed stream for the initial guesses for the material on the feed tray. Initializing the feed tray then gives us initial guesses for the state of the vapor and liquid streams leaving the feed tray. The flows of liquid and vapor streams entering the feed tray are set to a very small fraction of the feed flow, and all other state variables are set to the same values of the feed. Next, the rectification section is initialized starting with the condenser (if present). As an initial guess, the vapor leaving the feed tray is used as the initial guess for the inlet to the condenser, which is then initialized to give an initial guess for the reflux stream being recycled to the top tray. All the trays in the rectification section are then initialized using the conditions of the vapor stream leaving the feed tray for the vapor inlet and the conditions of the reflux stream for the liquid inlet (i.e., all rectification trays are initialized at the same conditions). Once all the rectification trays have been initialized, the entire rectification section (not including the feed tray and condenser) is solved using the initial guesses for the vapor stream from the feed tray and reflux stream from the condenser. A similar approach is then used to initialize the stripping section, starting with the reboiler 15 Figure 5. Hierarchical flow of information during initialization and using the guesses for the liquid stream leaving the feed tray. The boil-up stream from the reboiler is then used as the initial guess for the vapor inlet on each of the stripping trays and the liquid stream from the feed tray for the liquid inlets. All the stripping trays are initialized at these conditions, after which the entire stripping section is solved (again without the reboiler or feed tray). Note that for the feed tray and trays in the rectification section and stripping section, the initialization is applied to the state blocks first, followed by the mass, energy, and momentum balances for the tray. With the feed tray and rectification and stripping sections now initialized, the entire tray column is initialized next, keeping the initial guesses for the reflux and boilup streams. Fi- nally, the model is solved with the condenser included, followed by the full system, including the reboiler. 3.2.3 Link to Flowsheet Initialization Successful initialization of an EO-flowsheet depends on the successful initialization of the unit models and associated state blocks. The initialization routines described earlier for the state blocks and the unit model blocks are critical pieces in the flow of state information when initializing a flowsheet. The hierarchical flow of information in IDAES-IP is illustrated in Figure 5, where steps 1, 2, and 3 show how information on state cascades from flowsheet-level to state blocks, and steps 4, 5, and 6 show how successful initialization of the sub-blocks cascades state information to the next unit in the flowsheet. The initialization routines for the sub-blocks need to be robust and fast, especially for the state blocks with the thermophysical properties because of the number of state blocks that will exist in a large flowsheet. Section 4 shows the advantages of the proposed VLE formulation in terms of initialization speed. 16 4 Computational Results We demonstrate our proposed single-stage flash formulation with thermodynamics based on Peng-Robinson CEOS for two different component systems. Specifically, we consider flash calculations based on the equation-oriented framework described in the equation system (13) for a benzene-toluene binary mixture as well as a nitrogen-oxygen-argon ternary mixture. The latter represents the component mixture involved in cryogenic air separation units. We compare the performance of our proposed approached for these cases with the approach of Burgard et al. [18] in terms of tractability and speed of solution. It is important to note that while methods from earlier works can be applied to these examples, they are mixed inequality/equality systems that require intrusion with VLE models, which makes them ill- suited for modular thermodynamic property routines. On the other hand, the approaches in this work and the work of Burgard et al. [18] are non-intrusive methods due to their modular nature, as the smooth, square flash model can be embedded directly in process models for any thermodynamic property model without altering the structure of the process model or requiring specialized SM or optimization-based subroutines. Besides flash calculations, we demonstrate application of the proposed VLE formulation in optimizing the operation of a tray column in IDAES-IP. The proposed formulation along with case studies will be available at https://github.com/IDAES/idaes-pse/tree/main/idaes. 4.1 P Q-flash simulations for a benzene-toluene mixture In the first case study, we consider an equimolar feed of a binary mixture comprising benzene and toluene. The temperature of the feed is 368 K, and we consider a range of pressures from 0.1Pc in the subcritical region to 1.2Pc in the supercritical region. PQ -flash simulations are performed at each pressure for a parametric sweep of fifty heat duties ranging from 0 J to 70,000 J. Figure 6 shows the variation of the vapor fraction at the flash outlet with temperature for each of these simulations. At low duties, the mixture is entirely in the liquid phase. As more heat is supplied and the bubble point of the mixture is attained, the mixture becomes two-phase until its dew point is reached. As more heat is supplied beyond the dew point, the mixture exists entirely in the vapor phase. As the PQ -flash sweeps are performed at higher pressures, the two-phase condition is realized at higher temperatures. Beyond the critical pressure of the mixture, the underlying assumption is that the behavior of the mixture is liquid-like. This is also illustrated in Figure 6, where the vapor fraction of the mixture is zero when the flash pressure is 1.2 Pc for any heat duty supplied. Similarly, Figure 7 shows the corresponding variation of the heat duty with the flash temperature. As pressure is increased, the two-phase region, demarcated by points of discontinuity, keeps shrinking until it eventually disappears at supercritical pressures. To demonstrate the reliability, tractability, and speed of convergence of the proposed flash formulation in equation system (13), we compare its performance with the formulation in 17 400 450 500 550 600 650 Temperature [K] 0.0 0.2 0.4 0.6 0.8 1.0 Vapor fraction [-] 0.1 Pc 0.2 Pc 0.3 Pc 0.4 Pc 0.5 Pc 0.6 Pc 0.7 Pc 0.8 Pc 0.9 Pc 0.95 Pc 0.99 Pc 0.999 Pc 1.2 Pc Figure 6. Variation of vapor fraction with temperature for the benzene-toluene mixture 400 450 500 550 600 650 Temperature [K] 0 10000 20000 30000 40000 50000 60000 70000Heat duty [J] 0.1 Pc 0.2 Pc 0.3 Pc 0.4 Pc 0.5 Pc 0.6 Pc 0.7 Pc 0.8 Pc 0.9 Pc 0.95 Pc 0.99 Pc 0.999 Pc 1.2 Pc Figure 7. Temperature variation with heat duty for the benzene-toluene mixture Burgard et al. [18]. We perform PQ -flash simulations with parametric sweeps of heat duty for pressures ranging from 0.1 Pc to 2.0Pc and compare the two formulations based on three metrics: successful convergence, speed of initialization, and speed of simulation. The simulations are performed on an Intel ® Core™ i7 CPU @ 2.11 GHz processor with 24 GB memory. Figure 8 compares the convergence of the simulations for the two formulations using 18 the large-scale nonlinear equation solver IPOPT [25]. Each square in the grid represents a 0 7368 14737 22105 29474 36842 44211 51579 58947 66316 Heat duty [J] 2.0 1.5 1.2 1.0 0.999 0.99 0.95 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 P/P c [-] Burgard et al. (2018) 0 7368 14737 22105 29474 36842 44211 51579 58947 66316 Heat duty [J] This work Unsuccessful Successful Figure 8. PQ -flash simulation convergence comparison for the benzene-toluene mixture particular instance ofPQ -flash. The size of a problem instance is 45 equations and variables, and the average solve time is 0.004 CPU seconds. As illustrated, a grey square indicates that the simulation was completed successfully, while a black square implies the solver did not converge to the solution. We note that the formulation in Burgard et al. [18] fails to converge for pressures either close to or above the critical pressure, whereas the proposed formulation converges successfully for all subcritical and supercritical pressures. We can attribute the superior performance of the proposed formulation, especially near the critical point, to the complementarities associated with the cubic roots in equation system (13). Moreover, the calculation of the bubble and dew points in the previous work is absent in our proposed method, leading to lower computational complexity. Additionally, in the earlier work, when we are close to the critical point, the bubble and dew temperatures and composition converge to the same values leading to degeneracy in the problem. This is naturally avoided in the proposed formulation. Figure 9 compares the speed of the two formulations in their solution and initialization stages, respectively. The proposed formulation is almost twice as fast as the earlier work in solving the PQ -flash simulations and more than twice as fast during the initialization stage. This can be directly attributed to lower complexity in the proposed formulation as the embedded calculation of the bubble and dew points is obviated. 4.2 P Q-flash simulations for an air separation mixture This case study involves a ternary system comprising nitrogen (N 2), oxygen (O2), and argon (Ar) with feed mole fractions of 0.60, 0.35, and 0.05, respectively. The temperature of the feed is 95 K. We consider the same range of pressures as in the previous case. The simulations 19 0 20 40 60Time [s] Time required to solve 50 P Q-flash simulations 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99 0.999 1.0 1.2 1.5 2.0 P/P c [-] 0 5 10Time [s] Time required to initialize P Q-flash Burgard et al. (2018) This work Figure 9. Comparison of simulation and initialization speeds for the benzene-toluene mixture are performed for each pressure for fifty heat duties ranging from -1000 J to 7000 J. As done for the benzene-toluene mixture, Figure 10 shows the variation of the vapor fraction at the flash outlet with temperature, and Figure 11 shows the variation of heat duty with temperature. The same performance metrics are used to compare the proposed formulation with the formulation proposed by Burgard et al. [18] for the ternary mixture. Figure 12 compares the convergence for the two formulations. The problem size of an instance of PQ - flash is 50 equations and variables, and the average solve time is 0.006 CPU seconds. Again, the proposed formulation successfully converges in all cases using IPOPT, while the earlier work fails in the neighborhood of the critical point. Figure 13 illustrates the solution and initialization speeds of the two formulations. 4.3 Optimization of a tray column It is not uncommon to encounter disappearing phases when performing equation-oriented optimization either as part of the optimization iterations or at the optimal solution. In this case study, we demonstrate the numerical robustness of the proposed VLE formulation when used in equation-oriented optimization of a tray column where the disappearing phases will not result in discontinuities and infeasibilities. Here, we use the tray column model from the IDAES-IP unit model library with the Peng-Robinson CEOS property package for the 20 80 100 120 140 160 180 Temperature [K] 0.0 0.2 0.4 0.6 0.8 1.0 Vapor fraction [-] 0.1 Pc 0.2 Pc 0.3 Pc 0.4 Pc 0.5 Pc 0.6 Pc 0.7 Pc 0.8 Pc 0.9 Pc 0.95 Pc 0.99 Pc 0.999 Pc 1.2 Pc Figure 10. Vapor fraction variation with temperature for the air separation mixture 80 100 120 140 160 180 Temperature [K] −1000 0 1000 2000 3000 4000 5000 6000 7000 Heat duty [J] 0.1 Pc 0.2 Pc 0.3 Pc 0.4 Pc 0.5 Pc 0.6 Pc 0.7 Pc 0.8 Pc 0.9 Pc 0.95 Pc 0.99 Pc 0.999 Pc 1.2 Pc Figure 11. Temperature variation with heat duty for the air separation mixture nitrogen-oxygen-argon component system described above. The unit model consists of the material balance, equilibrium condition, summation condition, and heat balance (MESH) equations augmented by the smoothed complementarities in (13). Specifically, the optimiza- tion problem minimizes the reboiler duty subject to three different purity specifications:
• case 1: xDistillate N2 ≥ 0.90 and xBottoms O2 ≥ 0.75 21 -1000 -158 684 1526 2368 3211 4053 4895 5737 6579 Heat duty [J] 2.0 1.5 1.2 1.0 0.999 0.99 0.95 0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1 P/P c [-] Burgard et al. (2018) -1000 -158 684 1526 2368 3211 4053 4895 5737 6579 Heat duty [J] This work Unsuccessful Successful Figure 12. PQ -flash simulation convergence comparison for the air separation mixture 0 25 50 75Time [s] Time required to solve 50 P Q-flash simulations 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.95 0.99 0.999 1.0 1.2 1.5 2.0 P/P c [-] 0.0 2.5 5.0 7.5 Time [s] Time required to initialize P Q-flash Burgard et al. (2018) This work Figure 13. Comparison of simulation and initialization speeds for the air separation mixture
• case 2: xBottoms O2 ≥ 0.60
• case 3: xBottoms O2 ≥ 0.44 22 For these three cases, the number of trays in the column is fixed at 30, with the first tray at the top of the column and tray 30 at the bottom of the column. The 30 trays do not include the condenser and reboiler. The feed is introduced at tray 15, and the condenser specifications are set to a total condenser and at the same pressure as the feed. The decision variables for this optimization problem are the reflux ratio in the condenser and the boilup ratio in the reboiler; both variables are bounded between 0 and 5. Each case comprises approximately 2200 equations and variables and is solved on an Intel® Core™ i9 CPU @ 2.60 GHz processor with 32 GB memory. The average computation time for the optimization problems, not including initialization, is approximately 12 CPU seconds. The initialization strategy for all three cases is the same and follows the sequence outlined in Section 3.2.2.2 where the feed inlet is used to initialize the trays starting with the feed tray first. The optimal reboiler duty, reflux ratio, and boilup ratio are shown in Table 1, and Figure 14 shows the liquid and vapor flow rates along the column. For case 1, the high-purity speci- fications result in a higher reboiler duty, and non-zero liquid and vapor flows are observed across the column as reflux and vapor are recirculated through the column to improve sepa- ration. For case 2, where the purity specification for the nitrogen in the distillate is removed and reduced for the oxygen in the bottoms product, the optimal reflux ratio drops to zero and the reboiler duty decreases. As no reflux is coming into the top tray, there is no liquid flow for the trays in the rectification section, i.e., trays 1 to 14. The zero flows are more pronounced in case 3, as shown in Figure 14, where the purity specification for oxygen was set to the same condition that would be obtained if the two-phase feed were separated in a single flash drum. This low-purity specification can be attained simply by letting the vapor and liquid from the feed tray be the distillate and bottoms product streams, respectively. There is no requirement to recirculate reflux or vapor through the column, which leads to zero liquid flows for trays in the rectification section and zero vapor flows for trays in the stripping section. These case studies show the versatility of the proposed VLE formulation for handling disappearing phases that can be encountered during EO-based distillation opti- mization. Since the smooth, square formulation in the proposed approach can be embedded directly in EO-based optimization without altering the structure of the model, it can be adopted to generic EO-based flowsheets and applied with different thermodynamic property routines. In contrast to optimization-based approaches such as [6], it doesn’t require an objective function to treat the underlying nonlinear complementarity problem. 5 Conclusions Newton-based solution strategies are essential for large-scale simulation and optimization problems, and modern NLP solvers can now handle optimization problems with millions of variables. These solvers rely on smooth models as well as exact first and second derivatives. Moreover, several optimization modeling frameworks, including IDAES-IP, provide these derivatives automatically. 23 Table 1. State variables, heat duties, and optimal reflux and boilup ratios for the three cases Case 1 Case 2 Case 3 xDistillate N2 ≥ 0.90 xBottoms O2 ≥ 0.60 xBottoms O2 ≥ 0.44 xBottoms O2 ≥ 0.75 Distillate Bottoms Distillate Bottoms Distillate Bottoms Flow rate (mol/s) 5.9079 4.0921 7.5021 2.4979 5.3866 4.6134 xAr 0.02706 0.08312 0.04321 0.07039 0.04321 0.05793 xN2 0.90000 0.16688 0.69003 0.32961 0.69003 0.49488 xO2 0.07294 0.75000 0.26676 0.60000 0.26675 0.44720 Temperature (K) 101.15 111.48 103.49 108.51 103.49 106.00 Pressure (MPa) 0.79 0.79 0.79 0.79 0.79 0.79 Heat duty (kJ) -61.14 34.56 -37.41 10.59 -26.86 2.38× 10−6 Optimal reflux or boilup ratio 1.224 1.542 1.09× 10−6 0.813 2.11×10−6 1.03× 10−5 10 20 30 Tray number 6 8 10 12Flow rate [mol/s] Liquid Vapor (a) Case 1 10 20 30 Tray number 0 2 4 6Flow rate [mol/s] Liquid Vapor (b) Case 2 10 20 30 Tray number 0 2 4 6Flow rate [mol/s] Liquid Vapor (c) Case 3 Figure 14. Liquid and vapor flows in the column However, VLE calculations derived from free energy minimization present a challenge for these optimization solvers, as they are inherently procedural, in order to accommodate non- smooth behavior that requires discrete decisions. Alternatively, these nonsmooth elements can be handled as complementarity conditions, and the corresponding optimization prob- lem can be formulated as an MPCC. A popular MPCC approach that is compatible with large-scale NLP solvers is through the application of a successive smoothing strategy. This MPCC approach has been the subject of several optimization studies that focus on VLE and distillation as described in Section 1. To generalize this approach beyond optimization modeling platforms, this study develops a “square-flash” formulation based on complementarity conditions. This leads to a modular framework that generally applies to most CAPE simulation platforms. The approach allows a direct interface to general thermodynamic packages, as is demonstrated here on CEOS models. In particular, we provide a detailed description of how the square-flash formulation integrates into the architecture of the IDAES-IP platform and lends itself to conventional initialization strategies, leading to efficient and robust simulation and optimization. Our proposed approach is tested on two VLE case studies with pressure sweeps that seamlessly 24 transition from liquid to two-phase to vapor to supercritical regimes. The approach is also demonstrated on a ternary distillation that is able to converge with dry trays below the minimum reflux ratio. These cases cannot be obtained with current commercial tools as they are not equipped to relax the flash equations dynamically in response to phase changes. Future work will extend our proposed square-flash approach to explore process systems with more complex multi-phase thermodynamic properties. We will also consider distributed systems such as multi-phase pipeline and complex heat exchanger models. Our approach can also be applied to simultaneous dynamic optimization strategies for multi-phase systems, where VLE represents a large portion of the dynamic model. Finally, we note that for these optimization problems, the modular square-flash approach lends itself to parallel solution strategies, which will also be explored in the future. Acknowledgments This work was conducted as part of the Institute for the Design of Advanced Energy Systems (IDAES) with support through the Simulation-Based Engineering, Crosscutting Research Program and the Solid Oxide Fuel Cell Program’s Integrated Energy Systems thrust within the U.S. Department of Energy’s Office of Fossil Energy and Carbon Management. Disclaimer: This project was funded by the United States Department of Energy, National Energy Technology Laboratory, in part, through a site support contract. Neither the United States Government nor any agency thereof, nor any of their employees, nor the support con- tractor, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any informa- tion, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or ser- vice by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. 25 References [1] H. H. Rachford and J. D. Rice. Procedure for use of electronic digital computers in cal- culating flash vaporization hydrocarbon equilibrium. Journal of Petroleum Technology , 4(10):19–3, 1952. [2] J. F. Boston and H. I. Britt. A radically different formulation and solution of the single-stage flash problem. Computers & Chemical Engineering , 2(2-3):109–122, 1978. [3] H. A. J. Watson, M. Vikse, T. Gundersen, and P. I. Barton. Reliable flash calculations: Part 1. nonsmooth inside-out algorithms. Industrial & Engineering Chemistry Research, 56(4):960–973, 2017. [4] L. Qi and J. Sun. A nonsmooth version of Newton’s method.Mathematical programming, 58(1):353–367, 1993. [5] F. Facchinei, A. Fischer, and M. Herrich. An LP-Newton method: nonsmooth equations, KKT systems, and nonisolated solutions. Mathematical Programming, 146(1):1–36,

[6] R. S. Kamath, L. T. Biegler, and I. E. Grossmann. An equation-oriented approach for handling thermodynamics based on cubic equation of state in process optimization. Computers & Chemical Engineering , 34(12):2085–2096, 2010. [7] L.G. Bullard and L. T. Biegler. Iterative linear programming strategies for nonsmooth simulation: Phase equilibrium applications. Computers and Chemical Engineering , 17, 1:95, 1993. [8] V. Gopal and L. T. Biegler. Smoothing methods for complementarity problems in process engineering. AIChE journal , 45(7):1535–1547, 1999. [9] Y-D Lang and L. T. Biegler. A distributed stream method for tray optimization. AIChE J., 48, 3:582, 2002. [10] L. T. Biegler. Nonlinear programming: concepts, algorithms, and applications to chem- ical processes. SIAM, 2010. [11] Z-Q Luo, J-S Pang, and D. Ralph. Mathematical programs with equilibrium constraints. Cambridge University Press, 1996. [12] M. C. Ferris and J-S Pang. Complementarity and variational problems: State of the art. Proceedings of the International Conference on Complementarity Problems SIAM, Philadelphia, 1997. [13] A. Raghunathan and L. T. Biegler. MPEC Formulations and Algorithms in Process Engineering. Computers and Chemical Engineering , 27:1381–1392, 2003. 26 [14] B. T. Baumrucker, J. G. Renfro, and L. T. Biegler. MPEC Problem Formulations in Chemical Engineering Applications. Computers and Chemical Engineering , 32:2903– 2913, 2008. [15] A. U. Raghunathan, M. S. Diaz, and L. T. Biegler. An MPEC formulation for dynamic optimization of distillation operations. Computers & chemical engineering, 28(10):2037– 2052, 2004. [16] A. W. Dowling and L. T. Biegler. A framework for efficient large scale equation-oriented flowsheet optimization. Computers & Chemical Engineering , 72:3–20, 2015. [17] A. W. Dowling, C. Balwani, Q. Gao, and L. T. Biegler. Optimization of sub-ambient separation systems with embedded cubic equation of state thermodynamic models and complementarity constraints. Computers & Chemical Engineering , 81:323–343, 2015. [18] A. P. Burgard, J. P. Eason, J. C. Eslick, J. H. Ghouse, A. Lee, L. T. Biegler, and D. C. Miller. A smooth, square flash formulation for equation-oriented flowsheet optimization. In Computer Aided Chemical Engineering , volume 44, pages 871–876. Elsevier, 2018. [19] A. Lee, J. H. Ghouse, J. C. Eslick, C. D. Laird, J. D. Siirola, M. A. Zamarripa, D. Gunter, J. H. Shinn, A. W. Dowling, D. Bhattacharyya, L. T. Biegler, A. P. Burgard, and D. C. Miller. The idaes process modeling framework and model library—flexibility for process simulation and optimization. Journal of Advanced Manufacturing and Pro- cessing, 3(3):e10095, 2021. [20] B. E. Poling, J. M. Prausnitz, and J. P. O’Connell. Properties of gases and liquids . McGraw-Hill Education, 2001. [21] T. Gundersen. Numerical aspects of the implementation of cubic equations of state in flash calculation routines. Computers & Chemical Engineering , 6(3):245–255, 1982. [22] F. M. Khoury. Calculate the right density. Hydrocarbon Processing, 6:155–157, 1978. [23] E. Zhao and S. Saha. Applications of complex domain in vapor- liquid equilibrium calculations using a cubic equation of state. Industrial & engineering chemistry research, 37(5):1625–1633, 1998. [24] D. M. Gay. Hooking Your Solver to AMPL . AMPL Optimization, Inc., Albuquerque, NM, 2017. [25] A. W¨ achter and L. T. Biegler. On the implementation of an interior-point filter line- search algorithm for large-scale nonlinear programming. Mathematical programming, 106(1):25–57, 2006. 27 Appendix A Analysis of the nonlinear complementarity problem (NCP) formulation for flash calculations In this section, we analyze the NCP problem described in Section 2. Specifically, we show that the proposed formulation in equation set (13) follows from minimizing the Gibbs free energy of the system. A.1 Basic NCP formulation The complementarity formulation maintains the structure of the external CEOS functions. A temperature relaxation T is introduced that models disappearing phases: F =L +V (18.1) Fzi =Lxi +Vyi, i ∈C (18.2) FhF =LhL (x,T,P ) +VhV (y,T,P )−Q (18.3) yi =Ki

T,P, x, y
xi (18.4) 0 = X i∈C yi−xi (18.5) T−T =T +−T− (18.6) 0≤T +⊥V≥ 0 (18.7) 0≤T−⊥L≥ 0 (18.8) Replacing complementarity constraints with smooth approximations yields: F =L +V (19.1) Fzi =Lxi +Vyi, i ∈C (19.2) FhF =LhL (x,T,P ) +VhV (y,T,P )−Q (19.3) yi =Ki

T,P, x, y
xi (19.4) 0 = X i∈C yi−xi (19.5) T−T =T +−T− (19.6) 0 =T + +V− q (T +−V )2 +ε2 (19.7) 0 =T− +L− q (T−−L)2 +ε2 (19.8) 28 A.2 Cubic Equations of State Cubic equations of state are desirable models for VLE because they represent physical phe- nomena with reasonable accuracy and require only moderate computational complexity for a variety of non-polar and slightly polar compounds (with appropriate mixing rules and interaction parameters). For instance, they have been used for decades in the oil and gas industries for hydrocarbon systems. The most general cubic equation of state is given by f(Z) =Z 3− (1 +B−uB)Z 2 + (A +wB 2−uB−uB2)Z −AB−wB 2−wB 3 = 0 (20) with parameters A = am(z)P R2T 2 and B = bm(z)P RT and mixture parameters am(z) and bm(z) dependent on the flash feed mole fractions z. In addition, u = 1,w = 0 for Soave-Redlich-Kwong (SRK) and u = 2,w = −1 for Peng- Robinson (PR) models. When the cubic roots are determined, we have three cases:

• Three roots: ZL <Z s <Z V with f′(ZL),f′(ZV )> 0 and f′′(ZL)< 0,f′′(ZV )> 0.
• Two roots: ZL <Z V with andf′′(ZL)< 0,f′′(ZV )> 0 and eitherf′(ZL)> 0,f′(ZV ) = 0 or f′(ZL) = 0,f′(ZV )> 0.
• Single root: For the single real root, we have f′′(Z)< 0 or f′′(Z)> 0. (f(Z) = 0 and f′′(Z) = 0 implies f(Z) is not a cubic equation.) and we define ZL =ZV . Thus, the incorporation of CEOS into the NCP includes inequality constraints used to select the appropriate root for each phase ( ZL for liquid, ZV for vapor) with: f′(ZL) = 3Z 2 L− 2(1 +B−uB)Zl +A +wB 2−uB−uB2≥ 0 (21.1) f′′(ZL) = 6ZL− 2(1 +B−uB)< 0 (21.2) f′(ZV ) = 3Z 2 V− 2(1 +B−uB)ZV +A +wB 2−uB−uB2≥ 0 (21.3) f′′(ZV ) = 6ZV− 2(1 +B−uB)> 0 (21.4) Note that f(Z) = 0 and f′′(Z) = 0 defines an inflection point that cannot represent the extreme roots (ZL,ZV ), so strict inequalities apply for (21.2) and (21.4). This leads to the following NCP with CEOS constraints. F =L +V (22.1) 29 Fzi =Lxi +Vyi, i ∈C (22.2) FhF =LhL (x,T,P ) +VhV (y,T,P )−Q (22.3) yi =Ki

T,P, x, y
xi (22.4) 0 = X i∈C yi−xi (22.5) T−T =sV−sL (22.6) 0≤sV⊥V≥ 0 (22.7) 0≤sL⊥L≥ 0 (22.8) f′′(ZL) =g+ L−g− L (22.9) f′′(ZV ) =g+ V−g− V (22.10) 0≤g− V⊥V≥ 0 (22.11) 0≤g+ L⊥L≥ 0 (22.12) Since (22.8) and (22.12) both complement L, they can be combined to form 0 ≤g+ L +sL⊥ L≥ 0, which appears in the final square-flash formulation (13). A.3 Properties of CEOS formulation This section presents the properties that describe the relationship between roots of the cubic function and the phase flows for different cases. Property 1: ZL̸=ZV =⇒L,V ≥ 0. Proof: Since ZL̸=ZV implies f′′(ZL)< 0,f′′(ZV )> 0 we have f′′(ZL) =g+ L−g− L < 0 and g− L > 0, and also f′′(ZV ) =g+ V−g− V > 0 and g+ V > 0. a) Now consider that g+ L > 0. Then from (22) L = 0,V = F and ZL = ZV , which leads to f′′(ZL) = f′′(ZV ) > 0. This is a contradiction, which implies that g+ L = 0, f′′(ZL) < 0 but L = 0 may still hold, i.e., at the dew point. b) Similarly, consider g− V > 0. Then, again from (22), V = 0,L = F and ZV = ZV , which leads to f′′(ZV ) = f′′(ZL) < 0, a contradiction. This implies that g− V = 0, f′′(ZV ) > 0 and V = 0 may still hold, i.e., at the bubble point. Property 2: L,V > 0 =⇒ZL̸=ZV . Proof: From (22) L,V > 0 implies that g+ L,g− V = 0. Since neither f′′(ZL) nor f′′(ZV ) can be zero, we must have from (22) that f′′(ZL) =−g− L < 0 and f′′(ZV ) = g+ V > 0, yielding ZL̸=ZV . Property 3: ZL =ZV =⇒L = 0,V =F or V =F,L = 0. Proof: ZL =ZV implies that f′′(ZL) =f′′(ZV ) and we first consider f′′(ZL) =f′′(ZV )> 0. Here f′′(ZL) = g+ L−g− L > 0, g+ L > g− L≥ 0 and L = 0 from (22). On the other hand, when f′′(ZL) =f′′(ZV )< 0, then f′′(ZV ) =g+ V−g− V < 0, g− V >g + V≥ 0 and V = 0. 30 A.4 Justification for Nested Formulation Flash separations operate at vapor-liquid equilibrium (i.e., boiling mixtures) between so- called dew and bubble point conditions. Outside of these ranges, one of the phases disap- pears and equilibrium no longer holds between the phases. To model these phase changes, Gibbs minimization at a fixed temperature T and pressure P is written as a constrained optimization problem of the form: min li,vi G(T,P,l i,vi) = NCX i=1 li ¯GL i + NCX i=1 vi ¯GV i (23) s.t. li≥ 0, v i≥ 0 li +vi =mT i > 0, i= 1,...,NC where NC refers to the the number of chemical components with index i, G(T,P,l i,vi) is the total Gibbs free energy ¯GL i = ¯Gig i (T,P ) +RT ln(fL i ) ¯GV i = ¯Gig i (T,P ) +RT ln(fV i ) ¯Gig i is the ideal gas free energy per mole for component i, ϕL i and ϕV i are the mixture liquid and vapor fugacities for component i, li and vi are the moles of component i in the liquid and vapor phase, R is the gas constant, and mT i are the total moles of component i. The first order KKT conditions for this problem are given by: ¯Gig i (T,P ) +RT ln(ϕL i ) + (NCX j=1

lj ∂ ¯GL j ∂li +vj ∂ ¯GV j ∂li !) −λi−γi = 0 (24.1) ¯Gig i (T,P ) +RT ln(ϕV i ) + (NCX j=1

lj ∂ ¯GL j ∂vi +vj ∂ ¯GV j ∂vi !) −fi−γi = 0 (24.2) 0≤λi⊥li≥ 0 (24.3) 0≤fi⊥vi≥ 0 (24.4) li +vi =mT i, i= 1,...,NC. (24.5) The bracketed terms in (24.1) and (24.2) are equal to zero from the Gibbs-Duhem equation. Subtracting (24.1) from (24.2) leads to: RT ln(ϕV i /ϕL i )−fi +λi = 0. 31 Moreover, defining ϕV i =ϕV i (T,P,y i)yi, ϕL i =ϕL i (T,P,x i)xi and Ki =ϕL i/ϕV i where ϕL i and ϕV i are fugacity coefficients, leads to: yi = exp fi−λi RT  Kixi. (25) From (24.3) and (24.4) we can deduce 0 ≤ λi⊥ fi≥ 0 and by defining βi = exp(fi−λi RT ) we have fi > 0 =⇒ βi > 1, λ i > 0 =⇒ βi < 1 leading to the following complementarity system: yi =βiKixi (26.1) βi− 1 =β+ i −β− i (26.2) 0≤li⊥β− i ≥ 0 (26.3) 0≤vi⊥β+ i ≥ 0. (26.4) In this manner, phase existence can be determined within the context of an MPCC. If a slack variable (β+ i orβ− i ) is positive, either the corresponding liquid or vapor phase is absent and βi̸= 1 relaxes the phase equilibrium condition, as required in (24). Now consider the unconstrained case: min li,vi G(T,P,l i,vi) = NCX i=1 li ¯GL i + NCX i=1 vi ¯GV i (27) s.t. li +vi =mT i, i= 1,...,NC with the KKT conditions: ¯Gig i (T,P ) +RT ln(ϕL i )−γi = 0 (28.1) ¯Gig i (T,P ) +RT ln(ϕV i )−γi = 0 (28.2) li +vi =mT i, i= 1,...,NC. (28.3) Subtracting (28.1) from (28.2) leads to RT ln(ϕV i /ϕL i ) = 0, which applies for equilibrium at a temperature, say T , between bubble and dew point conditions, thus leading to: yi = Ki(T,P, x, y)xi. Now using the same compositions as in (25), and assuming that dKi dT > 0, we can compare (25) with the temperature relaxation and write: yi = exp fi−λi RT  Ki(T,P, x, y)xi =βiKi(T,P, x, y)xi (29) 32 =Ki(T,P, x, y)xi = Ki(T,P, x, y) Ki(T,P, x, y)Ki(T,P, x, y)xi (30) and assign βi =Ki(T,P, x, y)/Ki(T,P, x, y). Since dKi dT > 0, we can replace βi− 1 by T−T and derive the following complementarity system: yi =Ki(T,P, x, y)xi (31.1) ¯T−T =T +−T− (31.2) 0≤li⊥T−≥ 0 (31.3) 0≤vi⊥T +≥ 0. (31.4) In this manner, phase existence can be determined within the context of an MPCC. If a slack variable (T + orT−) is positive, either the corresponding vapor or liquid phase is absent and ¯T̸= T relaxes the phase equilibrium condition. Note that this is equivalent to the system derived in equation set (18). 33

Источник

Источник

Теги

• #AspenHYSYS
• #PhaseEnvelope
• #Troubleshooting
• #Термодинамика

Источник: https://industrialmonitordirect.com/blogs/knowledgebase/aspen-hysys-dew-point-vs-bubble-point-calculation-discrepancies

#AspenHYSYS #PhaseEnvelope #Troubleshooting #Термодинамика