A Unifying Approach for the Calculation of ... - ACS Publications

Subscriber access provided by UNIV OF NEBRASKA - LINCOLN

Article

A unifying approach for the calculation of azeotropes and pinch points in homogeneous and heterogeneous mixtures. Mirko Skiborowski, Jürgen Bausa, and Wolfgang Marquardt Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.6b01303 • Publication Date (Web): 03 Jun 2016 Downloaded from http://pubs.acs.org on June 12, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

A unifying approach for the calculation of azeotropes and pinch points in homogeneous and heterogeneous mixtures. Mirko Skiborowski,∗,† Jürgen Bausa,†,‡ and Wolfgang Marquardt†,¶ Department of Biochemical and Chemical Engineering, Laboratory of Fluid Separations, TU Dortmund University, 44227 Dortmund, Germany E-mail: [email protected]

Abstract This paper presents a unifying approach for a reliable calculation of all azeotropes and pinch points for homogeneous and heterogeneous non-reactive mixtures. The approach builds on a reformulation of the pinch-equation system and an efficient continuation algorithm including bifurcation analysis. The calculation of azeotropes is based on an analogy of univolatility curves and the pinch branches for pure component products. The resulting method can efficiently and reliably calculate all homogeneous and heterogeneous non-reactive azeotropes. The only exception would be the existence of isolated univolatility curves inside the composition space, which to the best knowledge of the authors has never been reported. Since the method can further determine the ∗

To whom correspondence should be addressed Department of Biochemical and Chemical Engineering, Laboratory of Fluid Separations, TU Dortmund University, 44227 Dortmund, Germany ‡ AVT - Process System Engineering, RWTH Aachen University, 52064 Aachen, Germany, Current address: Covestro Deutschland AG, 51365 Leverkusen, Germany ¶ AVT - Process System Engineering, RWTH Aachen University, 52064 Aachen, Germany, Current address: Forschungszentrum Juelich GmbH, 52425 Juelich, Germany †

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 58

pinch branches and pinch points based on the knowledge of the azeotropes, it presents a self-contained method for calculation of the most important information on non-ideal multi-component mixtures to support process synthesis. The application of the method is illustrated by means of several complex mixtures.

1

Introduction

Back in the 1990s, distillation was used for about 95% of all fluid separations in the chemical industry, accounting for 3% of the world-wide and even 10% of the U.S. energy consumption 1 . According to current estimates not much has changed in the last 20 years 2 ; there are approximately 40.000 distillation columns in the U.S. 3 , which account for around 40% of the total energy consumption in the chemical industry 4 . Despite the evolution of alternative separation technologies and low thermodynamic efficiency 5 , distillation is still the preferred separation technology for high capacities and purities. Conceptual design of distillation processes for azeotropic mixtures remains a complex task, which has to be based on rigorous thermodynamics and account for the inherent limitations regarding split feasibility 6 . Split feasibility primarily depends on the existence of distillation boundaries, which directly relates to the existence of azeotropes. In order to model distillation-based separation processes, it is essential to predict whether a given mixture will form one or more azeotropes and to calculate the composition, temperature and the stability information for each azeotrope 7 . The meaning of the word azeotrope, which is derived from ancient Greek and translates to ”to boil unchanged” 8 , provides the mathematical definition xi = yi∗

, i = 1, . . . , nc ,

(1)

requiring that the liquid composition x matches its equilibrium vapor composition y ∗ . This is also the case for a heterogeneous mixture, for which both liquid phases (xI and xII ) are in


Page 3 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


equilibrium with the vapor phase and add up to a mixture that equals the vapor phase, i.e., ∗ xi = (1 − ϕ)xIi + ϕxII i = yi

, i = 1, . . . , nc ,

(2)

where ϕϵ[0, 1] represents the liquid phase ratio. While these relationships look quite simple, the determination of the existence of azeotropes and the computation of existing azeotropes is quite challenging, since it requires the calculation of all solutions of a strongly nonlinear set of algebraic equations for the calculation of thermodynamic equilibrium. In case of heterogeneous mixtures, also the correct determination of phase stability has to be accounted for, in order to discriminate between vapor-liquid equilibrium (VLE) and vapor-liquid-liquid equilibrium (VLLE). Since azeotropes are the fixed points of residue curve and distillation line maps, they define the topology of these maps 9 . Consequently, the calculation of all existing azeotropes is an essential step for many process synthesis methods 10–12 , as well as for distillation shortcut methods that require the definition of feasible product specifications 13 . Significant effort has been devoted to the development of reliable algorithms for the calculation of azeotropes and R modern process simulators like ASPEN Plus⃝ provide corresponding routines.

In addition to the calculation of all azeotropes, also the fixed points of the finite difference or a differential model for column profile calculations 14 play an important role in shortcut methods for conceptual process design. These fixed points are also termed pinch points, since the driving force for mass transfer between the phases becomes infinitesimally small 15 . The pinch points are the basis for different shortcut methods for the calculation of the minimum energy demand in single-feed distillation columns for the separation of homogeneous 16–20 and heterogeneous mixtures 21,22 , extractive distillation columns 23 and reactive distillation 24 . Although impressive results have been published for these methods they are not implemented in commercial process simulators. The pinch branches, which are the loci of all pinch points, are also an essential part in different



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

synthesis and analysis methods for the determination of feasible splits. The pinch lines, which are the pinch branches that are directly connected to a product composition, are equivalent to the composition profiles in a section of a reversible distillation column 17 . While residue curves and distillation lines relate to a separation with a maximum driving force, the pinch lines relate to one with minimum driving force. Their combination presents the envelope of all potential composition profiles in a column section 15 . These ”operation leafs” 11 can be used for the determination of feasible products in single-feed 15 and two-feed distillation columns 25 . Based on a linear approximation of the operational leafs for multicomponent mixtures Thong et al. 11,26 propose an automatic procedure for the synthesis of distillation column sequences for the separation of homogeneous multicomponent azeotropic mixtures. More recently, Petlyuk et al. 27–29 revived and extended their previous work on the method of infinitely sharp splits for the identification of feasible product specifications for single and two-feed columns, also making use of pinch branches. Although impressive results have been presented with these methods, none of these approaches is available in commercial process simulators either. Although the calculation of azeotropes is possible in most commercial process simulation software, this information is not further utilized and the addressed process synthesis and design methods are generally missing. One reason for this might be the absence of an efficient and reliable method for the calculation of pinch lines and pinch points, for which only very few methods have been proposed. To overcome this limitation, we present a novel continuation-based method, which builds on the equivalence of pinch branches for pure component products and univolatility curves. The resulting method facilitates a numerically robust and efficient computation of both, azeotropes and pinch points, in homogeneous and heterogeneous mixtures. The method is neither restricted to a certain number of components, nor limited to a specific thermodynamic model. However, the accuracy of the results of course heavily depends on the appropriate selection of a thermodynamic model. At first a concise review of available methods for azeotrope and pinch point calculation in


Page 4 of 58

Page 5 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


given in Section 2. Afterwards, Section 3 demonstrates the equivalence of pinch branches and univolatility curves and how these can be utilized for the computation of azeotropes in homogeneous and heterogeneous mixtures. Section 4 presents the continuation-based methods for the calculation of pinch branches, azeotropes and pinch points. The utilization of these approaches is illustrated for several case studies in Section 5, before presenting a final conclusion and outlook in Section 6.

2

Azeotrope and pinch-point calculation

In order to highlight the added value of the presented approach, the available methods for the calculation of azeotropes and pinch points are first reviewed in this section.

2.1

Calculation of azeotropes

The detection of homogeneous binary azeotropes is rather simple. As described by Jaubert and Privat 30 , ”a positive azeotrope is simply a minimum boiling-point azeotrope that is characterized by a minimum on an isobaric Txy diagram and by a maximum on the corresponding isothermal pxy diagram”. A negative azeotrope results from the opposite extremal points. This relationship has been described by Lecat 31 almost a century ago and allows to easily identify binary azeotropes from Txy and pxy diagrams, as illustrated in Figure 1. However, a visual analysis of these diagrams is no viable option for multicomponent mixtures, since the number of binary sub-systems increases exponentially. Therefore, simple criteria have been derived to determine the existence and location of binary diagrams. The most efficient criteria for the existence of azeotropes in binary mixtures rely only on the evaluation of the K-values 32 or the activity coefficients at infinite dilution 33 . The K-values at infinite dilution can be determined from the ratio of the vapor to the liquid composition ∞ ∞ = [y2 /x2 ]) 32 . The existence of a negative azeotrope = [y1 /x1 ], K21 for the trace species (K12

(maximum-boiling azeotrope) relates to both of these K-values being smaller than 1, while



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Txy diagrams for three binary mixtures at 1atm: acetone-chloroform (minimum azeotrope), chloroform-methanol (maximum azeotrope) and benzene-hexafluorobenzene (double azeotrope). the existence of a positive azeotrope (minimum-boiling azeotrope) relates to both K-values being larger than 1. These conditions are generally satisfied for the common case of a single minimum or maximum azeotrope, as illustrated for the first two mixtures in Figure 1. However, for the rare case of double azeotropy, as illustrated for the benzene-hexafluorobenzene mixture in Figure 1, this criteria fails. This might be the reason, why the ”azeotrope search” in the simulation R software ASPEN PLUS⃝ does not detect any azeotrope for this mixture (tested with ASPEN

Plus

TM

V7.2 and V8.8), although the Txy diagram looks exactly as illustrated in Figure 1.

Jaubert and Privat 30 presents further examples of double azeotropy or saddle azeotropes, for which these criteria fail. For multicomponent mixtures and ternary or higher azeotropes such simple criteria do not even exist. Several articles were published in the last two decades dealing with the problem in different ways. Fidkowski et al. 34 developed a continuation approach for the calculation of homogeneous azeotropes. The approach is based on a homotopy from the ideal equilibrium relationship towards the nonideal equilibrium relationship and a bifurcation analysis. Eckert and Kubicek 35 show that the continuation method is in principle capable of calculating homogeneous and heterogeneous azeotropes by a simple extension of the problem formulation. Wasylkiewicz et al. 36 finally present the full extension of the continuation approach to 6 ACS Paragon Plus Environment

Page 6 of 58

Page 7 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


compute homogeneous and heterogeneous azeotropes with a more general homotopy function. Tolsma and Barton 37,38 also present an extension of the method of Fidkowski et al. 34 to heterogeneous mixtures, which is independent of the representation of the nonideality of the system and can be readily extended to handle multiple liquid phases. They include a thorough theoretical derivation including several proofs for completeness, but demonstrate application only for ternary mixtures. Another interesting approach making use of continuation is a geometrical method by Qi and Sundmacher 39 . They propose to represent the mass balance conditions for azeotropes geometrically as reaction kinetic surfaces, which are further generalized as potential singular point surfaces (PSPS) that contain all pure components and azeotropes. The PSPS can be calculated from the pure components by means of continuation. Despite the calculation of nonreactive homogeneous and heterogeneous azeotropes, also reactive and kinetic azeotropes can be determined by means of this method, since reaction stoichiometry and thermodynamic models are considered by using the transformed compositions previously introduced by Barbosa and Doherty 40 . However, this geometric method is limited to ternary mixtures due to the necessary visualization. Maier et al. 41,42 present a different approach for the reliable computation of all homogeneous azeotropes based on interval analysis. They propose the use of an interval-Netwton/generalizedbisection-algorithm in order to provide a mathematical and computational guarantee, that all azeotropes are located. Salomone and Espinosa 43 adapt the method of Maier et al. 41 and combine the sequential formulation in a directed order with the topological consistency criterion of Zharov and Serafimov 44 , in form of the general topological equations presented by Pöllmann et al. 45 . The adaption avoids the numerical verification of the non-existence of azeotropes within subsystems, corresponding to a subset of the components in the mixture, and therefore increases the computational efficiency. Developing convex underestimations for different thermodynamic models, Harding et al. 46 propose the use of a branch and bound framework in combination with successive refinement and partitioning in order to determine



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 58

all homogeneous azeotropes by means of global optimization. The approach is based on the reformulation of enclosing all solutions of nonlinear systems of constrained equations into a global optimization problem as presented by Maranas and Floudas 47 . Harding and Floudas 48 extend their approach to heterogeneous and reactive systems in a further publication, including an additional check for the possibility of trivial solutions (xIi = xII i , i = 1, ..nC ). Due to the deterministic global optimization approach, the method theoretically guarantees enclosing all azeotropes. However, the method is computational demanding and requires specific convex underestimation functions for each thermodynamic model. While the authors demonstrate the supremacy of a sequential search for each subspace, requiring the existence of each component (N-ary azeotrope), over a search in the full composition space (k-ary azeotropes), the computational times reported for the single subspaces indicate a significant exponential increase in CPU time for an increasing number of components. The authors implemented the methods in the EQUISTAR software 7 package, which is however not available anymore. In contrast to a deterministic optimization approach, Bonilla-Petriciolet et al. 49 proposed a stochastic optimization based on simulated annealing, for the determination of homogeneous azeotropes. For strongly nonideal systems they advise to use a hybrid approach starting with simulated annealing, followed by a local deterministic method for final refinement. While the approach also aims at a global optimization, there is no theoretical proof of convergence.

2.2

Calculation of pinch points

While azeotropes are the fixed points of residue curve and distillation line maps, pinch points represent the fixed points of a column section operated at finite reflux. Azeotropes and pinch points are directly related, which can easily be seen from the tray-to-tray profile in a rectifying section of a simple distillation column assuming constant molar overflow (CMO) 6

yn+1 =

r 1 xn + xD , r+1 r+1


(3)

Page 9 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


with yn+1 being the composition of the entering vapor stream, xn the composition of the emanating liquid stream, xD the top product composition and r the reflux ratio. For total reflux (r → ∞) the equation simplifies to the distillation line equation

yn+1 = xn ,

(4)

indicating that azeotropes are nothing less than the pinch points for a column section operated at total reflux. This relation is also valid without the CMO assumption and is utilized in continuation-based approaches for the calculation of pinch points, by utilizing pure components and azeotropes as initial solutions. The resulting solution branches, i.e. the locus of all pinch points for varying reflux ratios, are consequently termed pinch branches. Fidkowski et al. 50 and Julka and Doherty 51 propose the use of a pseudo-arc-length, predictor-corrector algorithm of Keller 52 in order to compute the pinch branches and pinch points for the finite difference (tray-to-tray) column model. In contrast to that, Pöllmann and Blass 53 propose the integration of a differential equation for the profile of reversible distillation, which they derive by means of implicit differentiation. Most recently, Felbab 54 presented another approach based on a differential equation, but with a different parametrization, performing the integration for the reflux ratio instead of the temperature. The author claims that the method is simpler and less computationally intensive than that of Pöllmann and Blass 53 and presents several computational examples. However, the approach is limited by the CMO assumption and the application is only demonstrated for quaternary mixtures. Other approaches for the calculation of pinch points that are not based on continuation or integration include a deflation method and a hybrid sequential niche algorithm 55 , as well as an optimization based approach with a systematic refinement of the composition space 56 . However, these approaches are computationally much less efficient. Most importantly, all listed approaches have only been presented for homogeneous mixtures.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.3

Discussion

As illustrated in the review of the available methods, multiple approaches for the computation of all azeotropes have been presented, while only few approaches for pinch point computation have been presented. Although some of the methods for azeotrope computation are applicable to heterogeneous mixtures, none of the approaches for pinch point computation has been extended to heterogeneous mixtures. In addition, the methods for azeotrope computation are generally based on sets of equations, which are different from those for pinch point computation. The article of Felbab 54 is the only one that proposes to utilize the pinch branch calculation also to determine all azeotropes. However, the article does not establish any derivation for this approach and proposes to use a multitude of product specifications, logging the results into a database of found azeotropes to increase the probability of determining all azeotropes. This paper presents a self-contained method for the calculation of all existing azeotropes, as well as pinch branches for non-reactive homogeneous and heterogeneous systems. Although the approach also utilizes the pinch branches for azeotrope calculation, the idea of the approach is rather to use the univolatility curves as a tool to detect all azeotropes. By acknowledging the equivalence between univolatility curves and pinch branches for pure component products, which is demonstrated in Section 3, the pinch branch calculation provides us with exactly this tool. Apart from this theoretical derivation, the approach differs from the method of Felbab 54 . Most importantly, it uses a continuation method for the calculation of the pinch branches and a phase-stability test, based on an embedded homotopy-continuation for the calculation of the thermodynamic equilibrium 57 , which makes it applicable to homogeneous and heterogeneous systems. A reformulation of the pinch equation systems and a continuous monitoring of the systems Jacobian determinant improve numerical robustness, avoid branch switching and allow to locate bifurcations along the pinch branches. For the calculation of azeotropes, the equivalence of pinch branches for pure component products to univolatility curves is utilized in order to find all azeotropes based on a continuation 10 ACS Paragon Plus Environment

Page 10 of 58

Page 11 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


using only the pure components as initial solutions. After all azeotropes are determined the very same continuation method is used to determine the pinch branches for an arbitrary product composition and the pinch points for a specific reflux. The presented approach can be applied in combination with any thermodynamic model for VLE or VLLE calculation.

3

Theoretical derivation

In their tutorial paper on the synthesis of ideal and non-ideal distillation-based separation systems, Westerberg and Wahnschafft 58 first described the idea for determining azeotropes in multicomponent mixtures by means of a continuation along the univolatility curves. Univolitility is based on the relative volatility

αij =

yi /yj Ki = , xi /xj Kj

(5)

which is a measure of the difference in volatility between two components i and j. The relative volatility indicates how easy or difficult a separation of these components by means of distillation is expected to be. In case of univolatility, the relative volatility is equal to one and the K-values as well as the volatility of both components are equal (Ki = Kj ). Azeotropes represent a special case of univolatility, for which the K-values of the present components are not only equal, but also equal unity. Since it is well known that the existence of azeotropes directly entails the existence of univolatility curves/surfaces connecting the azeotrope with the boundaries of the subspaces 9,59–61 , utilizing the univolatility curves for determining the azeotropes is a straight forward idea. Westerberg and Wahnschafft 58 proposed to start a continuation along the univolatility curves from binary azeotropes. However, instead of employing a specific homotopy for the continuation, they rather proposed an iterative procedure by a stepwise extension of the amount of an additional component while keeping the K-values of the components present in the initial point (current azeotrope) equal. While our approach also builds on the idea of using 11 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 58

the univolatility curves for the determination of all azeotropes, it considers a very different approach for the calculation of the univolatility curves, building on the pinch equation system.

3.1

Homogeneous mixtures

In case of a homogeneous system, the pinch equation system consists of the following mass and energy balances, summation constraints and equilibrium conditions: 0 = V · yi∗ − L · xi − D · xD,i

, i = 1, . . . , nc ,

0 = V · hV (y ∗ , T, p) − L · hL (x, T, p) − D · hD + QC , 0=

nC ∑

(6) (7)

xi − 1,

(8)

yi∗ − 1,

(9)

i=1

0=

nC ∑ i=1

0 = yi∗ − Ki (x, y ∗ , T, p) · xi

, i = 1, . . . , nc ,

(10)

Here D, xD and hD are the flow rate, molar composition and specific enthalpy of the top product. L, x and hL are the flow rate, molar composition and specific enthalpy of the liquid stream leaving the column section, while V , y ∗ and hV are the flow rate, molar composition and specific enthalpy of the vapor stream entering the column section. The vapor composition y ∗ is considered to be in equilibrium with the liquid composition x. QC represents the energy duty of the condenser. Solving this set of equations for all heat duties and a specific product composition xD results in the calculation of the pinch branches for the rectifying section with a specific product composition xD . In fact, since the energy duty is only present in the energy balance (Eq. (7)), the system of equations (6) and (8)-(10) can also be solved varying the ratio L/V to calculate all pinch branches. The corresponding heat duty can be determined afterwards by evaluating the energy balance (Eq. (7)) for a specific pinch point on the pinch branches. 12 ACS Paragon Plus Environment

Page 13 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In the special case of a pure component product, the pinch branches correspond to the binary edges connected to the pure component vertex and to all univolatility curves for the remaining components. This is a direct consequence of the mass balance for the pinch equation system (Eq. (6)). Evaluating this equation for an arbitrary pure component of a nC component mixture (w.l.o.g. we choose component 1) results in   1 x    1         0  x2   y∗       2 0 = V ·  .  − L ·  .  − D · . . .  .   .  .  .   .        0 xnC yn∗ C 

y1∗







(11)

Consequently, for all components i ̸= 1 the flow rate of the component in the liquid phase (Lxi ) equals the flow rate in the vapor phase (V yi∗ ) and the following relations hold: xi = yi∗ = 0

, ∀i ̸= 1 with xi = 0,

(12)

L yi∗ = Ki (x, y∗, T, p) = xi V

, ∀i ̸= 1 with xi > 0.

(13)

While the value of L/V is linked to the energy duty via the energy balance equation (Eq. (7)), the important result is, that in case of the pure component product each solution to the pinch equation system requires that any other component either vanishes (xi = 0) or that it has the same K-value as any other present component apart from the product. Consequently, all pinch branches for the pure product correspond to the binary edges connected to the pure product and univolatility curves connected to them. Figure 2 and 3 present an exemplary illustration of the pinch branches for a ternary system of acetone, ethanol and water at two different pressure levels. The resulting pinch branches for the three different pure components specified as product composition are illustrated separately. The thermodynamic property data for this and all other investigated mixtures is listed in the supporting material.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 2: Pinch branches for the acetone/ethanol/water (AEW) mixture at 1 atm for the three different pure component products: acetone (left), ethanol (center), water (right). Figure 2 presents the pinch branches at atmospheric pressure, at which only one binary azeotrope between water and ethanol is present. The calculated pinch branches correspond to the binary edges only, for ethanol and water as products. Only for acetone as product an additional branch is computed, that bifurcates from the binary acetone-water edge inside the ternary region and passes through the binary ethanol-water azeotrope. This branch corresponds to the univolatility curve for ethanol and water (KE = KW ).

Figure 3: Pinch branches for the acetone/ethanol/water (AEW) mixture at 10 atm for the three different pure component products: acetone (left), ethanol (center), water (right). At an elevated pressure of 10 atm, as shown in Figure 3, the system exhibits three binary and one ternary azeotrope. Therefore, for each pure product, one branch bifurcates from a binary edge inside the composition space passing through the ternary and the binary 14 ACS Paragon Plus Environment

Page 14 of 58

Page 15 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


azeotrope opposite to the product. Each of the branches corresponds to the univolatility curve of the two components that are not present within the product. The ternary azeotrope is located at the intersection of the three binary univolatility curves. Higher-order azeotropes (quaternary, quinterary, ...) are located at the intersection of the higher-order univolatility curves, which similar to the binary univolatility curves can be traced by the bifurcations on the lower order branches. Section 4 presents the continuation approach used to trace these branches for multicomponent mixtures. Further illustration is presented in the examples in Section 5.

3.2

Heterogeneous mixtures

In case of a heterogeneous mixture the pinch equation system is equal to that for a homogeneous system outside of the miscibility gap(s) and can be described by Eqs. (6)-(10). However, inside the miscibility gap the liquid phase splitting needs to be accounted for within mass and energy balances, as well as the equilibrium relationships: 0 = V · yi∗ − (1 − ϕ) · LI · xIi − ϕ · LII · xII i − D · xD,i

, i = 1, . . . , nc ,

(14)

0 = V · hV (y ∗ , T, p) − (1 − ϕ) · LI · hL (xI , T, p) − ϕ · LII · hL (xII , T, p) − D · hD + QC , 0= 0= 0=

nC ∑ i=1 nC ∑ i=1 nC ∑

(15)

xIi − 1,

(16)

xII i − 1,

(17)

yi∗ − 1,

(18)

i=1

0 = yi∗ − Ki (xI , y ∗ , T, p) · xIi

, i = 1, . . . , nc ,

(19)

0 = yi∗ − Ki (xII , y ∗ , T, p) · xII i

, i = 1, . . . , nc .

(20)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 58

Outside of a miscibility gap there is only one liquid phase (xI = xII ) and the pinch equation system (Eqs. (14)-(20)) is reduced to Eqs. (6)-(10), such that this set of equations is valid in the complete composition space. By introducing a total liquid composition (x), the mass balance can be formulated similar to the homogeneous case by Eq. (14) in conjunction with

0 = xi − (1 − ϕ) · xIi − ϕ · xII i

, i = 1, . . . , nc .

(21)

In addition, by introducing pseudo-homogeneous enthalpy values (hp,L ) 0 = hp,L − (1 − ϕ) · hL (xI , T, p) − ϕ · hL (xII , T, p),

(22)

and pseudo-homogeneous K-values (Kp ), and inserting Eqs. (19) and (20) in Eq. (21) 0 = Kp,i −

Ki (xI , y ∗ , T, p) · Ki (xII , y ∗ , T, p) (1 − ϕ) · Ki (xII , y ∗ , T, p) + ϕ · Ki (xI , y ∗ , T, p)

, i = 1, . . . , nc ,

(23)

the energy balance and equilibrium relationship can also be formulated, similar to the homogeneous case, by Eqs. (7) and (10). Consequently, the calculation of the pinch branches for a heterogeneous region can be performed similar to the homogeneous regions, using the pseudo-homogeneous K-values (Kp,i ) and varying the L/V ratio for the system of equations (6) and (8)-(10) in conjunction with equations (16)-(17), (19)-(21) and (23). By starting the analysis again with the mass balance (Eq. (6)) and evaluating this equation for an arbitrary pure component of a nC -component mixture (w.l.o.g. we take again component 1) the same results as for the homogeneous case are obtained. Consequently, for all components i ̸= 1 the flow rate of the component in the (overall) liquid phase (Lxi ) equals the flow rate in the vapor phase (V yi∗ ) and the following relations hold: xi = yi∗ = 0 yi∗ xi

= Kp,i =

L V

, ∀i ̸= 1 with xi = 0,

(24)

, ∀i ̸= 1 with xi > 0.

(25)


Page 17 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Thus, similar to the homogeneous region, all pinch branches correspond to either a binary edge or a univolatility curve. Outside of the miscibility gap(s), the pseudo-homogeneous K-values are equal to the K-values for the homogeneous system, while inside of a miscibility gap the pseudo-homogeneous K-values differ and need to be calculated from the K-values for the two liquid phases by Eq. (23). Therefore, a reliable phase stability test is required throughout the calculation of the pinch branches, in order to detect whether there is a phase split or not. However, since any point on the surface of the miscibility gap corresponds to a single liquid phase (ϕϵ{0, 1}), there is a smooth transition for the pseudo-homogeneous K-values and the total liquid composition 62,63 .

Figure 4: Pinch branches for the isopropanol/water/cyclohexane (PWC) mixture at 1 atm for the three different pure component products: isopropanol (left), water (center), cyclohexane (right). Similar to Figure 2 and 3 for the case of a homogeneous mixture, Figure 4 presents an exemplary illustration of the pinch branches for a heterogeneous ternary system of isopropanol, water and cyclohexane. The system exhibits three binary and one ternary azeotrope, whereas the ternary and the binary water-cyclohexane azeotropes are located inside of the type I miscibility gap 64 and therefore labeled as VLLE. As for the homogeneous mixture the pinch branches are presented by the binary edges, connected to the pure component product, as well as a bifurcating branch, which runs through the ternary azeotrope and the binary azeotrope on the binary edge opposite to the pure component product. The bifurcated branches represent the univolatility curves of the ternary mixture. While the K-values of the components 17 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 58

differ for the single liquid phases inside the miscibility gap, the pseudo-homogeneous Kvalues on the univolatility curve are equal for the according components. This is illustrated in Table 6 in Appendix A for the univolatility curve of isopropanol and water (Kp,P = Kp,W ), presenting exemplary compositions and K-values located on the univolatility curve. The presented derivation and the examples show that the pinch branches correspond to either binary edges or univolatility curves in case of pure component products. This is equally true for homogeneous and heterogeneous mixtures. The only difference is that inside a miscibility gap the phase splitting needs to be accounted for and pseudo-homogeneous K-values need to be considered for the equivalence in case of heterogeneous mixtures. The following section presents the implemented continuation method and the algorithms for azeotrope and pinch calculation.

4

Continuation-based calculation of pinch branches and azeotropes

In order to determine all azeotropes and pinch points a continuation approach for the pinch equation system is utilized. Based on this continuation approach the calculation of all azeotropes is performed and subsequently pinch points can be determined by the same approach utilizing the azeotropes as initial solutions.

4.1

Calculation of pinch branches

The continuation approach, which is described in the following, is restricted to non-reactive systems. However, the extension to reactive mixtures is possible by introducing additional source/sink terms, as will be discussed in an upcoming publication. The description is divided in three parts, starting with the derivation of a general and numerically robust model for the pinch equation system, followed by the presentation of a continuation approach for the calculation of the pinch branches and a final postprocessing step. 18 ACS Paragon Plus Environment

Page 19 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.1.1

Model derivation

While the theoretical derivation in Section 3 was illustrated for the balance equations of the rectifying section, the equations for the continuation approach are presented for the balance of a general column section as illustrated in Figure 5.

Figure 5: General column section and relation to rectifying and stripping sections for pinch calculation. The pinch equation system for the general column section consists of mass and energy balances, summation constraints and equilibrium conditions, taking into account normalized liquid and vapor flow rates l and v, as illustrated in Figure 5:

0 = v · yi + l · xi − zi

, i = 1, . . . , nC ,

0 = v · hV (y, T, p) + l · hL (x, T, p) − hN , 0= 0=

nC ∑ i=1 nC ∑

(26) (27)

xi − 1,

(28)

yi − 1,

(29)

i=1

0 = yi − Ki (x, y, T, p) · xi

, i = 1, . . . , nC .


(30)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 58

As already stated in Section 3, the energy balance is not required to calculate the pinch branches, but only necessary to calculate certain pinch points, depending on the net enthalpy hN . Therefore, only the mass balances, summation constraints and equilibrium conditions in Eqs. (26) and (28)-(30) need to be solved in order to calculate the pinch branches. However, in the current form, this set of equations is not well posed. While singular points (i.e. pure components and azeotropes) are always solutions to this system of equations, the normalized vapor flow rate v tends to infinity at these points. Felbab 54 consequently avoids the evaluation at singular points, starting the integration, which he performs as continuation, from pinch points at high reflux in the vicinity of the singular points. These initial solutions are computed as direct solutions of the pinch equation system at an infinitesimal offset of the singular point using a Newton method. While this procedure is rather cumbersome, the numerical difficulties at the singular points can be overcome by means of a reformulation of the pinch equation system, similar to the method described by Guddat et al. 65 for parametric optimization problems. For this reformulation, an additional parameter α is introduced and the mass balance Eq. (26) is replaced by

0 = v¯ · yi + ¯l · xi − α · zi

, i = 1, . . . , nC ,

0 = α2 + v¯2 − 1.

(31) (32)

The reformulated pinch equation system, Eqs. (28)-(32), with the reformulated flow rates v¯ = α · v and ¯l = α · l can be solved everywhere, including the singular points. While the transition from physical coordinates v and l to the new variables α, v¯ and ¯l is not unique, the mapping of (α, v¯) → v is unique. However, the characteristics of the system cause that, after transformation, every physical solution corresponds to two different solutions for the new variables. Since finally only the physical solution is of interest, the algorithm tries to circumvent multiple evaluations of a physical solution branch by keeping track of the investigated solution branches and by continuously testing if the current solution branch was


Page 21 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


already investigated for opposite signs of the new variables. In order to be applicable for homogeneous and heterogeneous mixtures, the continuation approach needs to test the liquid composition for phase stability and exchange the K-value (Ki ) in Eq. (30) with the pseudo-homogeneous K-value (Kp,i ) determined from Eq. (23). This can be accomplished by using another nested homotopy-continuation approach for the determination of the correct equilibrium solution (VLE or VLLE). Here the approach of Bausa and Marquardt 57 is used, which has already been successfully applied for the dynamic simulation of three-phase reactive batch distillation 66 and the optimization-based design of heteroazeotropic distillation processes 63 . However, any method for the determination of phase stability could be integrated, including global optimization methods for the determination of the minimum of Gibbs free energy. Refer to the articles of Kangas et al. 67 and Zhang et al. 68 for a review of global optimization approaches for phase equilibrium calculation. Similar to the application in optimization-based process design 63 , the main reason for the implementation of the homotopy continuation approach is its computational efficiency, which is an operational necessity for the multitude of evaluations required for the pinch branch calculations. The reformulated set of equations (28)-(32) corresponds to an under-determined equation system of the general form 0 = g(u), with u = (xT , y ∗ T , T, α, v¯, ¯l)T . Since the dimension of u is by one bigger than the dimension of g (dim(u) = 2nC + 4 and dim(g) = 2nC + 3), the solution set is one-dimensional and can be parameterized, e.g., using the transformed normalized vapor flow rate v¯. Note that the thermodynamic property models are inherently included and that the presented formulation is therefore valid independent of the choice of thermodynamic property models used to determine the K-values. Although the computation by means of this set of equations is applicable for heterogeneous mixtures, since Eqs. (16)-(20) and (23) are included in the pseudo-homogeneous K-values in case of liquid phase-splitting, it has proven to be numerically more efficient to perform the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 58

continuation inside a heterogeneous region by solving the system of equations

0 = v¯ · yi + ¯l · ((1 − ϕ)xIi + ϕxII ) − αzi 0=

nC ∑

, i = 1, . . . , nC ,

(33)

xIi − 1,

(34)

xII i − 1,

(35)

yi − 1,

(36)

i=1

0= 0=

nC ∑ i=1 nC ∑ i=1

0 = yi − Ki (xI , y, T, p) · xIi

, i = 1, . . . , nC .

(37)

0 = yi − Ki (xII , y, T, p) · xII i

, i = 1, . . . , nC .

(38)

0 = α2 + v¯2 − 1,

(39)

taking the single liquid phase compositions directly into account. The set of equations (33)(39) also corresponds to an under-determined equation system of the general form 0 = g(u), T T with u = (xIi , xII , y ∗ T , T, α, v¯, ¯l)T , having one degree of freedom (dim(u) = 3nC + 5 and i

dim(g) = 3nC + 4). Consequently it can be solved by the same parametrization as the set of equations (28)-(32). The following continuation approach therefore switches between the different equation sets whenever a change in phase stability is detected.

4.1.2

Continuation approach

The nonlinearity of the equation system often results in multiple roots, such that sophisticated solution strategies need to be applied. Besides interval mathematics and global optimization approaches, continuation with embedded bifurcation analysis is one possible method to efficiently calculate all roots of the system. The latter methods are computationally less demanding and capable of calculating complete solution branches and all bifurcations, if at least one point on all connected branches is known a-priori 69 .


Page 23 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6: UML activity diagram for the pinch branch continuation algorithm. A UML activity diagram 70 of the pinch branch continuation algorithm employed in this work is illustrated in Figure 6. The algorithm implements the continuation methodology described by Deuflhard and Hohmann 71 and Choi et al. 72 . In a first step, all known solutions are placed on a stack. This stack contains only the pure components, if azeotropes are calculated, and all the singular points, if the calculation of the pinch points is performed. The stack is processed iteratively. For each of the initial solutions, the pinch branch is calculated using a predictor-corrector method with an Euler-predictor and a Newton-corrector. Starting from a known solution un a standardized tangent vn to the solution branch is computed. Using the tangent and a step size s, a first linear estimate

uˆn+1 = un + s · vn


(40)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

for the next solution step is calculated (predictor step). Afterwards a Newton solver is used to determine the correct solution of the system 0 = g(un+1 ) based on the prediction uˆn+1 as initial solution (corrector step). Since the equation system is under-determined, an extended equation system needs to be used 72 . A pseudo-arc-length algorithm is used in the implementation, in which the corrector step is always orthogonal to the predictor step. There are different reasons for which the step size s is controlled. First of all, the step size is reduced in case the Newton solver fails to determine a feasible solution, which is specifically relevant in strongly nonlinear regions of the solution branch 71 . To assure that the algorithm does not switch solution branches the step size s is also reduced if a sign change of the deter∂g minant of the Jacobian of the extended system ( ∂u ) is detected. The sign change indicates

either the transition to an alternative solution branch or the occurrence of a bifurcation on the current branch. If a bifurcation is detected, the current part of the branch is stored as processed branch and the location of the bifurcation and the different directions are placed on the stack as starting points for additional branches. The direction of these branches result directly from an analysis of the kernel or null space (all eigenvectors with eigenvalue zero) of the Jacobian. The executed calculations show that the kernel is always two-dimensional. Therefore, there are always two tangents, one of which belongs to the actual solution branch and one that belongs to the bifurcating branch. Thus, there are three directions, for which additional pinch branches can be computed from the detected bifurcation. Although not highlighted in Figure 6 the step size s is also reduced in an internal loop if a change in phase-splitting is detected between two steps, in order to determine the exact intersection with the miscibility gap (ϕ = {0, 1}). Each pinch branch is pursued until the compositions reach a predefined limit. For physical meaningful solutions the compositions are obviously limited to [0, 1]. However, depending on the thermodynamic models, it is possible to evaluate the pinch branches also outside the physically meaningful region. The computation of composition profiles and pinch points outside the physically meaningful range has first been introduced by Tapp et al. 73 and Holland


Page 24 of 58

Page 25 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


et al. 74 , who show that these solutions can be of value for process synthesis and design. Section 5.4 presents an interesting example of the insight that can be gained by extending the calculation of the pinch branches to compositions outside the physically meaningful range, while Appendix B presents some further discussion on the consideration of these solutions.

4.1.3

Postprocessing

After calculation of the pinch branches the net enthalpy hn can be determined for any pinch solution by evaluating Eq. (27) for the stored solution vector u. This allows the location of pinch points associated to a specific heat duty of a distillation column section, making use of the transformations shown in Figure 5. Specific pinch point locations can easily be determined by searching the pinch branches for coinciding intervals. The exact location can be approximated by means of linear interpolation between the two known solutions and potential refinement by a Newton-corrector step. The second preprocessing step has the objective to determine the stability of the actual solution. This is of special importance for the application in pinch-based shortcut methods 16,18,23,51 and the determination of thermodynamical consistency in case of singular points (Section 4.2). To this end, the tray-to-tray recursion is interpreted as a function of the overall liquid phase composition x and pressure p. The local topology of the tray-to-tray recursion in the vicinity of the pinch points can be analyzed based on the eigenvalues and eigenvectors of the Jacobian ∂xn+1 /∂xn of the pinch equations 16,75 . The Jacobian has nC − 1 positive eigenvalues in case of a homogeneous and nC − 2 positive eigenvalues in case of a heterogeneous mixture 16,57 . The stability of the pinch points can be classified according to the values of the eigenvalues. While eigenvalues greater than unity correspond to departing profiles and an unstable eigenvector, eigenvalues smaller than unity correspond to stable ones and approaching profiles 16,53,76 . In accordance with Doherty 16,75 a pinch solution is classified as a stable node in case all eigenvectors are stable. If all eigenvectors are unstable the pinch solution corresponds to an unstable node, while the pinch solution corresponds to a sad-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

dle in case of stable and unstable eigenvectors. Since azeotropes are nothing else than the pinch solutions for a column section operated at total reflux, the same classification holds for azeotropes as well.

4.2

Calculation of azeotropes

After deriving the analogy between the pinch branches for pure component products and the univolatility curves, the presented continuation approach is now embedded into a procedure for the efficient and reliable determination of all azeotropes. This procedure is illustrated as UML activity diagram in Figure 7.

Figure 7: UML activity diagram for the azeotrope calculation algorithm. As a preliminary step a stack with pure component products and another stack with initial solutions (also pure components) for the pinch branch continuation are formed. In order to determine all azeotropes inside the composition space (0 ≤ xi ≤ 1) in principle one pure component product should suffice, assuming that higher order azeotropes are always located on the intersection of univolatility curves originating from each lower dimensional subspace. However, for the calculation of azeotropes outside the composition space (xi < 0 & xi > 1) a calculation of all univolatility curves, using all pure components for the calculation is ad26 ACS Paragon Plus Environment

Page 26 of 58

Page 27 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


vised. This generally increases the reliability, while it of course increases the computational effort. After taking one product composition from the stack the pinch branch continuation is performed, in order to determine the univolatility curves. These curves are afterwards searched for azeotropic compositions by comparing the liquid and vapor compositions, subject to some numerical inaccuracies, using the Euclidean norm:

∥x − y∥2 ≤ ϵ,

(41)

In order to account for the numerical inaccuracies and refine the azeotropic compositions and temperature, a Newton solver is used to solve the according equation systems consisting of the summation constraints, the equilibrium conditions and additional conditions for azeotropy, which utilize the knowledge from the previous step in two ways. While the known compositions and temperature values are used for initialization, the knowledge of the involved components is also used in the formulation of the azeotropy conditions, which are

0 = xi 0 = Ki (x, y, T, p) − 1

, component i is not present,

(42)

, component i is present.

(43)

This formulation prevents the solver from calculating a solution in a neighboring subspace 58 . It is similar to the sequential formulation used by Maier et al. 41 . In case of a heterogeneous mixture the pseudo-homogeneous K-value (Eq. (23)) is used in Eq. (43). After the refinement, possible duplicates are deleted. This is especially important in case more than one pure component product specification is used for the calculation of the univolatility curves. If the pure component product stack is empty, the algorithm terminates. Otherwise the pinch branch continuation is invoked with the next product composition, unless an optional topological consistency check is activated. Here, the equation proposed by Zharov and Serafimov 44,77 is used to determine the topological consistency based on the 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 58

eigenvalue information (see Section 4.1.3). A topologically consistent set of singular points satisfies the condition nC ∑

2k · (Nk+ + Sk+ − Nk− − Sk− ) = (−1)nC −1 + 1,

(44)

k=1

where Nk+ and Sk+ correspond to the number of k-component nodes and saddles with positive index, while Nk− and Sk− correspond to the number of k-component nodes and saddles with negative index. The index of the singular point is positive in case of an even number of stable eigenvectors and negative in case of an odd number. The topological consistency check was already used by Fidkowski et al. 34 for verification of the results, while Salomone and Espinosa 43 employ the consistency check to determine if the search for azeotropes has to be extended from lower dimensional subspaces to higher order regions of the composition space (e.g. into the quaternary region, after searching all ternary subspaces). Here, the topology consistency check can be used as an indicator, stopping further calculation, if a consistent set of singular points is determined. This criteria is however not sufficient for correctness. Inconsistent sets of singular points are frequently determined if the search for azeotropes is performed inside and outside of the physical composition space (cf. Appendix B). The presented approach has some distinct advantage in comparison to previous approaches. The computation of heterogeneous azeotropes is performed automatically in the pinch continuation. Also, a sequential search throughout every subsystem with an increasing number of components, as proposed by Westerberg and Wahnschafft 58 or used by Salomone and Espinosa 43 in combination with interval analysis, is not necessary. The search for azeotropes automatically continues to higher-dimensional regions of the composition space with the bifurcations along the pinch branches (univolatility curves). This also has the advantage, that no search in higher-dimensional regions is performed, if no univolatility curve originates in this region due to a bifurcation. Exceptional cases, like the double azeotrope illustrated in Figure 1, are determined without any problems.


Page 29 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


4.3

Calculation of pinch points

The calculation of pinch points builds on the presented pinch branch continuation and the azeotrope calculation. Since all singular points are solutions to the pinch equation system, independently of the specified net product and heat duties, all pure components and azeotropes are valid initial solutions for the computation of the pinch branches for an arbitrary product composition z 50,53 . Thus, in order to determine all pinch solutions, the pinch branch calculation is performed with the continuation method presented in Section 4.1 utilizing the pure components and all azeotropes as initial solutions. This may include those azeotropes, which are determined outside of the compositions space (cf. Section 5.4). The location of specific pinch points and the stability of these solutions is determined by the postprocessing calculations presented in Section 4.1.3.

5

Illustration

The following examples present the results of the pinch branch continuation and azeotrope calculation for mixtures of different complexity. Finally, the utilization of the method for the synthesis and design of distillation processes is illustrated. All computations have been R performed on a stand-alone PC with an an Intel⃝ Core

TM

i5 3230 M (2.60 GHz) using

R MATLAB⃝ R2011b. The pinch branch continuation is implemented in C and compiled as R mex-function, while the azeotrope computation is implemented directly in MATLAB⃝ .

5.1

Comparison with Felbab’s approach

In order to analyse the performance of the pinch branch continuation a comparison with the most recent approach by Felbab 54 is presented for the two most complex mixtures investigated in that study. This method is chosen as the benchmark, because the author claims that the suggested algorithm is simpler, faster and numerically more robust than all previous approaches. It is also the most recent approach. Since the algorithm is implemented 29 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 58

R R in Fortran, compiled as MATLAB⃝ mex-function and run on an Intel⃝ Core

TM

i5 430 M

(2.53 GHz), a direct comparison is feasible. 5.1.1

Benzene, hexafluorobenzene and n-hexane

The ternary mixture of benzene, hexafluorobenzene and n-hexane has a special topology, due to the double azeotrope for benzene and hexafluorobenzene, which was already introduced in Section 2.1 (cf. Figure 1). The results of the azeotrope computation for the ternary mixture, assuming pure n-hexane as net product are illustrated in Figure 8. Remember that the choice of the net product composition determines the pinch branches and consequently the univolitility curves that are calculated (cf. Section 3).

Figure 8: Pinch branches and azeotropes calculated for the benzene/hexafluorobenzene/nhexane (BHnH) mixture at 1 atm for a pure n-hexane product. Interestingly, but as expected, the two binary azeotropes (AZ (BH)1 and AZ (BH)2 ) are located on two separate univolatility curves for benzene and hexafluorobenzene (KB = KH ). Note, that the topology of this mixture is inconsistent, according to Eq. (44). As already described in Section 2.1 and also noted by Felbab 54 , ASPEN PLUS

TM

”azeotrope search”

fails in the attempt to calculate these two azeotropes (tested with ASPEN Plus 30 ACS Paragon Plus Environment

TM

V7.2 and

Page 31 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


V8.8). However, the cource of the residue curves in a ternary map indicate the existence of the two azeotropes on the binary benzene-hexafluorobenzene edge. Felbab 54 reports a CPU time of 0.329 sec, which is about four times larger than the 0.066 sec required with the current method. Note, that all CPU times reported for the current method represent the mean value of 50 consecutive evaluations.

5.1.2

Acetone, chloroform, methanol and benzene

The second example considers a quaternary mixture of acetone, chloroform, methanol and benzene at atmospheric pressure, which exhibits six azeotropes. The results computed with the NRTL model are listed in Table 1, providing information on composition, number of stable EV, type and temperature for each SP. Table 1: Acetone, chloroform, methanol and benzene azeotropes (NRTL) at 1atm. SP 1 2 3 4 5 6 7 8 9 10

acetone chloroform methanol 1.000 0.778 0.315 0.340 0.345

1.000 0.227 0.658 0.232 0.655

1.000 0.222 0.606 0.437 0.342 0.428 -

benzene

EV

type

T[◦ C]

1.000 0.394 0.021 -

1 1 3 3 0 2 1 0 2 2

S S SN SN UN S S UN S S

56.13 61.09 64.53 80.12 55.23 57.96 57.12 53.73 57.13 64.14

Figure 9 presents an exemplary illustration of the pinch branches for a continuation utilizing pure benzene as product specification. However, any pure component product suffices to detect all azeotropes. As illustrated, higher-order azeotropes (e.g. Az(ACMB)) are located on higher-order univolatility curves (KA = KC = KM ), which bifurcate from lower-order univolatility curves (KC = KM ). Felbab 54 reports a computational time of 0.117 sec for the computation of the pinch branches and azeotropes. The current method requires 0.265 sec.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 9: Pinch branches and azeotropes calculated for the quaternary acetone/chloroform/methanol/benzene mixture at 1 atm for pure benzene as the product. The comparison for the two complex mixtures shows that the current implementation is at least competitive in comparison to the method of Felbab 54 , providing an efficient and reliable computation of the pinch branches and azeotropes in homogeneous mixtures.

5.2

Application to heterogeneous mixtures

In order to demonstrate the application of the suggested computational method to heterogeneous mixtures, for which no alternative method for pinch branch computation has been reported, we investigate the most complex mixtures treated in the context of azeotrope computation.

5.2.1

N-propanol, water and toluene

The calculation of the azeotropes for the ternary mixture of n-propanol, water and toluene was investigated by Wasylkiewicz et al. 36 , who showed the importance of the correct determination of phase stability. In this case, the application of a homogeneous model results in a topologically consistent characterization with three binary azeotropes. The use of a 32 ACS Paragon Plus Environment

Page 32 of 58

Page 33 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


heterogeneous model employing the correct determination of phase stability also results in a topologically consistent system, but revealing two binary VLE azeotropes, as well as one binary and one ternary VLLE azeotrope in the type I miscibility gap 64 .

Figure 10: Pinch branches and azeotropes calculated for the propanol/water/toluene mixture at 1 atm for a pure n-propanol net product.

ternary

n-

Figure 10 illustrates the pinch branches and azeotropes determined by the method suggested in this work using n-propanol as net product for the computation. All azeotropes and the phase stability are determined correctly by the method in a CPU time of 0.035 sec, using the NRTL model. 5.2.2

Ethanol, water, cyclohexane and benzene

The quaternary mixture of ethanol, water, cyclohexane and benzene is the most complex mixture investigated in prior publications on the calculation of heterogeneous azeotropes 36,48 . The mixture exhibits five homogeneous and six heterogeneous azeotropes, including a quaternary heterogeneous azeotrope at atmospheric pressure. Table 2 lists all singular points, as well as their characterization determined by the current method for the NRTL model. 33 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 58

Table 2: Ethanol, water, cyclohexane and benzene azeotropes (NRTL) at 1atm. SP 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

ethanol water cyclohexane 1.000 0.895 0.445 0.281 0.441 0.423 0.321 0.301 0.277 0.286

1.000 0.105 0.191 0.299 0.152 0.162

1.000 0.558 0.461 0.450 0.527 0.699 0.337 0.355

benzene

EV

type

stability

T[◦ C]

1.000 0.554 0.527 0.701 0.116 0.550 0.386 0.197

3 3 3 3 2 2 1 2 2 1 2 1 2 1 2

SN SN SN SN S S S S S S S S S S UN

VLE VLE VLE VLE VLE VLE VLLE VLLE VLE VLE VLE VLLE VLLE VLLE VLLE

78.30 61.09 100.01 80.12 78.14 67.70 64.02 69.34 64.91 64.76 77.53 62.38 69.48 67.60 61.90

The resulting pinch branches and azeotropes from a computation using pure ethanol as product composition are illustrated in Figure 11.

Figure 11: Pinch lines and azeotropes for the heterogeneous ethanol/water/cyclohexane/benzene (EWCB) mixture and ethanol as net product.


Page 35 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The method determines all azeotropes and their stability correctly in 0.126 sec. Even though a fair comparison with the results of Wasylkiewicz et al. 36 and Harding and Floudas 48 can only be made, if the original codes were run on the same machine, a comparison of the CPU times reported in the original publications, taking into account the benchmarking results given in Appendix C, suggests that the current method is at least competitive. The computation is slightly faster then the updated 0.157 sec (= 40sec ) for the method of Wasylkiewicz 254 et al. 36 and significantly faster then the updated 1.47 sec (=

221sec ) 150

for the approach of

Harding and Floudas 48 .

5.3

Multicomponent mixtures

Since the current method is completely algorithmic it can be applied to mixtures with five or more components. Case studies with a five and seven component mixture are presented in the following.

5.3.1

Acetone, chloroform, methanol, ethanol and benzene

The quinary mixture of acetone, chloroform, methanol ethanol and benzene represents the most complex homogeneous mixture for which results on the computation of azeotropes are reported 34,41,43,46,49 . It exhibits six binary, two ternary and one quaternary azeotrope at atmospheric pressure. For a complete comparison with previously reported results, the current method is applied for three different GE models. The necessary CPU time varies with model complexity and is 0.735 sec for the NRTL model, 0.470 sec for the UNIQUAC model and 0.070 for the Wilson model. Obviously the calculations with the simpler Wilson model are much faster than with the more complex UNIQUAC and NRTL models. For all three models the compositions and temperatures of the computed azeotropes were validated with ASPEN PLUS

TM

”azeotrope search”, using the same parameters. For com-

parison with previous results CPU times are again updated according to the speed up ratios described in Appendix C. The results are summarized in Table 3. 35 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 58

Table 3: Comparison of CPU times for different methods (reported = values from the article, updated = updated values, based on benchmarking values (cf. Table 7). article

Wilson reported updated

Maier et al. 41 (T -dependent) 105.7 sec Maier et al. 41 (Tref ) 14.1 sec 46 Harding et al. 16.0 sec Salomone and Espinosa 43 58.1 sec 49 Bonilla-Petriciolet et al.

0.43 0.06 0.11 0.61

NRTL reported updated

sec 608.7 sec sec 273.0 sec sec 47.8 sec sec 1.01 sec

2.46 sec 1.10 sec 0.32 sec 0.06 sec

The computational results with the current method are generally competitive and agree well with the updated values for the approach of Maier et al. 41 , evaluating the activity coefficients for a fixed temperature (Tref ). While this simplification results in a significant reduction in CPU time, it may also produce topologically wrong results 41 . Even after the optimistic update of the reported CPU times, only the CPU time reported by Bonilla-Petriciolet et al. 49 for a combination of simulated-annealing and a Newton solver is superior to our method.

5.3.2

A complex mixture of seven components

Finally a mixture of seven components, consisting of n-hexane (n-H), ethanol (E), n-butyraldehyde (n-BA), water (W), n-butanol (n-B), isobutanol (i-B) and o-xylene (o-X) is considered. The thermodynamical equilibrium of this test system is modeled by the NRTL model in conjunction with the Redlich-Kwong equation of state for the description of non-idealities of the R vapor phase. Property parameters are taken from ASPEN Plus⃝ (Version 24.0), including

the estimation of missing property parameters by means of the UNIFAC model. For this system 13 VLLE and 8 VLE azeotropes are determined by means of the current method. The topological consistency of the classification of the system is confirmed by means of Eq. (44). All singular points, the stability and their characterization are listed in Table 4.


Page 37 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 4: Singular points of the complex mixture of seven components at 1atm. SP 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28

n-H

E

1.000 1.000 0.965 0.930 0.614 0.662 0.338 0.280 0.790 0.633 0.257 0.552 0.158 0.041 0.525 0.757 0.784 0.895

n-BA

W

n-B

1.000 0.386 0.720 0.155 0.698 0.733 0.304 -

1.000 0.210 0.110 0.135 0.261 0.267 0.171 0.207 0.686 0.638 0.210 0.753 0.687 0.794 0.105

1.000 0.035 0.818 0.006 0.247 0.172 -

i-B

o-X

EV

type

stability

T[◦ C]

5 3 3 6 6 5 6 4 5 3 4 2 1 2 4 1 0 1 2 1 2 4 3 3 5 4 5 2

S S S SN SN S SN S S S S S S S S S UN S S S S S S S S S S S

VLE VLE VLE VLE VLE VLE VLE VLE VLE VLE VLE VLE VLE VLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLLE VLE

68.73 78.31 74.88 100.02 117.75 107.68 144.29 68.59 116.94 68.17 107.64 62.68 58.35 72.24 61.43 56.32 55.77 67.38 67.49 57.34 61.26 90.35 88.32 61.43 92.64 90.11 93.67 78.15

1.000 1.000 0.182 0.070 0.965 0.035 0.036 0.314 0.245 0.117 0.141 0.206 -

R ASPEN Plus⃝ ”azeotrope search” identifies 20 of the 21 azeotropes (tested with ASPEN

Plus

TM

V7.2). Only SP 18, the ternary ethanol-water-n-butyraldehyde azeotrope is not

confirmed. The univolatility curves and azeotropes for this ternary mixture are illustrated in Figure 12. As expected, the ternary azeotrope is located at the intersection of the three univolatility curves which each pass through one of the binary azeotropes. While this ternary R azeotrope is not detected by ASPEN Plus⃝ ”azeotrope search” the remaining singular points

are correctly characterized as saddles and stable nodes, such that an unstable node is absent,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

which already indicates that the calculations are corrupted. The current method correctly determines all 21 azeotropes in 1.29 sec.

Figure 12: Univolatility curves and azeotropes for the ternary subsystem of ethanol, water and n-butyraldehyde at 1 atm.

5.4

Computations outside the composition space

While the calculation of azeotropes outside of the composition space (for mole fractions > 1 and < 0) has no physical meaning for the separation of a given mixture at a fixed pressure, such azeotropes can traverse inside the composition space at modified pressure. Depending on the type of thermodynamic model (cf. Appendix B), the continuation approach can extend the calculation outside of the composition space in order to determine ”inactive” azeotropes (having no impact on the topology of the mixture at the current pressure). Figure 13 shows the univolatility curves determined by the continuation approach for all three pure net products of a ternary mixture of acetone, ethanol and water, as well as the active and inactive azeotropes outside of the composition space, at different pressures. 38 ACS Paragon Plus Environment

Page 38 of 58

Page 39 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 13: Univolatility curves and azeotropes (active (AZ) and inactive (AZIA ) for different pressures of 1 atm (left) and 5 bar (right)). While four azeotropes are determined for each pressure, three of the azeotropes are inactive at 1 atm. All four azeotropes are active at 5bar, resulting in a significant change in the topology of the mixture, connected with the development of a third distillation region, which can be seen from the trajectory of the simple distillation boundary (SDB) at the different pressure values. The inactive azeotropes move inside of the composition space as the intersections of the two univolatility curves (KA = KE and KA = KW ) and the two binary edges for acetone and ethanol, as well as acetone and water. Therefore, analysis of the univolatility curves for varying pressure also provides information on the movement of azeotropes. Figure 14 illustrates the trajectories of the different azeotropes with increasing pressure, determined by means of a simple sensitivity analysis. Obviously the ethanol-water azeotrope (AZ(EW )) is rather insensitive to pressure, while the other three azeotropes move significantly with pressure. While this simple sensitivity analysis can produce valuable insight, a continuation approach with the pressure as continuation parameter, as presented by Wasylkiewicz et al. 78 presents a more efficient means to perform such calculations. The current 39 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

method can provide the necessary initial solutions for such an efficient continuation approach. In addition, inactive azeotropes are also valid solutions for the determination of pinch points. Consequently, they can be used for the continuation of pinch branches for arbitrary product compositions and may facilitate the determination of additional pinch branches.

Figure 14: Trajectories of the different azeotropes (active and inactive) for varying pressures between 1 atm and 20 bar.

5.5

Process synthesis and design

Finally the application of the method for synthesis and design is illustrated by means of a case study on the separation of the quaternary mixture of acetone, chloroform, benzene and toluene. The mixture exhibits only one azeotrope that separates the composition space into two distillation regions. However, several authors have shown that the simple distillation boundary (SDB) can be crossed to some extent at finite reflux 15,53,79 . Different sequences for the separation of an equimolar mixture have been proposed e.g. by Thong and Jobson 80 and Kraemer et al. 81 , exploiting the special curvature of the distillation boundary. The


Page 40 of 58

Page 41 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


configuration with the lowest energy requirement reported by Kraemer et al. 81 is illustrated in Figure 15. In order to achieve the necessary purity of 99 mol% for each component a toluene recycle is used. The necessary recycle flow rate is determined by an optimization problem minimizing the overall heat duty, subject to a feasibility constraint, which is derived from the pitchfork bifurcation distillation boundary (PDB) 82 and the rectification body method (RBM) 18 as described by Br¨ uggemann and Marquardt 83 .

Figure 15: Process configuration proposed by Kraemer et al. 81 for the separation of an equimolar mixture of acetone, chloroform, benzene, and toluene. By means of the current method the necessary recycle flow rate can be determined based on a simple bisection, while the computation of the pinch branches and pinch points provides the essential information required to determine the minimum energy demand based on the RBM 18 . In a first step, the equimolar feed stream is mixed with pure toluene. Both products, the acetone-rich distillate and the lean bottoms product are located in the same distillation region, if the pinch lines for both products connect the same stable and unstable nodes. Thus, the minimum recycle flow rate is determined such that the pinch line of the bottoms product for the mixed feed stream connects the acetone and toluene vertex. An illustration of the shift of the pinch line from one into the other distillation region for two different recycle flow ratios is enclosed in the supporting information.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 42 of 58

Figure 16: Pinch branches and rectification bodies for the first column in the process configuration (cf. Figure 15). The resulting recycle flow ratio of R/F = 0.52 is in excellent agreement with the optimized value reported by Kraemer et al. 81 . Further, using the pinch information determined from the current method, the minimum energy requirement of each column is determined based on the RBM criteria. The pinch branches and the rectification bodies for the rectifying section (RBRS ) and the stripping section (RBSS ) of the first column are illustrated in Figure 16, while the product streams and heat duties for each column are listed in Table 5. The result from this sequential computation is in very good agreement with the optimization-based results presented by Kraemer et al. 81 . Table 5: Product streams and required heat duties for the three columns in the process flowsheet shown in Figure 15. distillate bottoms product heat duties [mol/s] [mol/mol] [mol/s] [mol/mol] QB [kW] QC [kW] column 1 column 2 column 3 process

2.5 [1.0 0.0 0.0 0.0] 5.0 [0.0 0.5 0.5 0.0] 2.5 [0.0 1.0 0.0 0.0]

12.7 7.7 2.5

[0.0 0.197 0.197 0.606] [0.0 0.0 0.0 1.0] [0.0 0.0 1.0 0.0]


475 406 266 1147

462 389 265 1116

Page 43 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


6

Conclusions

The current article presents a unifying approach for the computation of all azeotropes and pinch points in homogeneous and heterogeneous mixtures. It is based on the combination of a continuation approach for the computation of the pinch branches and the equivalency of the pinch branches and univolatility curves for pure components as net product specifications. Numerical robustness is improved by means of a reformulation of the pinch equation system, which can be solved inside as well as outside the composition simplex, depending on the thermodynamic model, and which does not become singular at the singular points. By embedding a phase stability test into the continuation approach, it is also the first reported approach for the computation of pinch branches in heterogeneous mixtures. Several case studies demonstrate the reliability and competitiveness of the unifying method with previously reported specialized methods for the computation of pinch branches and azeotropes. Besides this, the current method presents a general tool which provides all the necessary information required for the application of topological synthesis methods and feasibility analysis 10,84,85 as well as pinch-based shortcut methods 16,18 . The application to synthesis and design is illustrated for the separation of an azeotropic quaternary mixture. Calculation and examination of the univolatility curves outside of the composition space also facilitates the determination of ”inactive” azeotropes, which can traverse inside of the composition space for varying pressure. These azeotropes present important information concerning the feasibility of pressure-swing distillation processes and represent additional starting points for pinch branch continuation. While it is generally assumed that the existence of azeotropes always entails the existence of the according univolatility curves 9,61 , a proof is missing. The derivation of such proof could build up on the work of Doherty and Perkins 86 , who showed that ternary azeotropes are either connected to binary azeotropes or a pure component vertex by a univolatility curve. Future work will also focus on the extension to reactive mixtures.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 44 of 58

Supporting Information Available The supporting information contains information on the thermodynamic property models and parameters used for the calculations in the current article, as well additional illustrations for Section 5.5. This material is available free of charge via the Internet at http://pubs. acs.org/.

A

Pseudo-homogeneous K-values

The definition of the pseudo-homogeneous K-values given in Eq. (23) is derived, starting with Eq. (21)

0 = xi − (1 − ϕ) · xIi − ϕ · xII i

, i = 1, . . . , nc ,

and introducing the general VLE conditions for each of the liquid phases (Eqs. (19)-(20)), 0 = yi∗ − Ki (xI , y ∗ , T, p) · xIi

, i = 1, . . . , nc ,

0 = yi∗ − Ki (xII , y ∗ , T, p) · xII i

, i = 1, . . . , nc ,

resulting in yi∗ yi∗ − ϕ · Ki (xI , y ∗ , T, p) Ki (xII , y ∗ , T, p) (1 − ϕ) · Ki (xII , y ∗ , T, p) + ϕ · Ki (xI , y ∗ , T, p) ⇔ 0 = xi − yi∗ Ki (xI , y ∗ , T, p) · Ki (xII , y ∗ , T, p) Ki (xI , y ∗ , T, p) · Ki (xII , y ∗ , T, p) ⇔ 0 = yi∗ − xi (1 − ϕ) · Ki (xII , y ∗ , T, p) + ϕ · Ki (xI , y ∗ , T, p) 0 = xi − (1 − ϕ) ·

, i = 1, . . . , nc , , i = 1, . . . , nc , , i = 1, . . . , nc .

The calculation requires the correct determination of the composition of the separate liquid phases (xI and xII ), which is performed by the homotopy-continuation approach of Bausa and Marquardt 57 . Table 6 lists several compositions and K-values for the bifurcated pinch


Page 45 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


branch in Figure 4 (right). While the K-values for isopropanol and water for each of the liquid phases deviate significantly, the pseudo-homogeneous K-values for these components are equal throughout the pinch branch, illustrating the equivalence to the univolatility curve. The compositions and the K-values for the ternary and binary azeotrope on this univolatility curve are highlighted in bold letters. Table 6: Exemplary compositions and K-values for the bifurcated pinch branch for the pure cyclohexane product - Figure 4 (right).

B

xP

xW

yP

yP

ϕ

xI P

xI W

xII P

xII P

I KP

I KW

II KP

II KW

Kp,P

Kp,W

0.009 0.025 0.034 0.052 0.065 0.090 0.111 0.129 0.152 0.185 0.236 0.256 0.268 0.347 0.426 0.503 0.521 0.521 0.524 0.547 0.614 0.671 0.691

0.051 0.058 0.063 0.072 0.079 0.094 0.107 0.119 0.136 0.162 0.200 0.216 0.224 0.282 0.339 0.394 0.407 0.407 0.405 0.393 0.355 0.321 0.309

0.048 0.112 0.138 0.176 0.196 0.221 0.234 0.242 0.247 0.250 0.255 0.256 0.257 0.261 0.264 0.266 0.267 0.267 0.271 0.305 0.438 0.611 0.691

0.286 0.264 0.256 0.243 0.237 0.229 0.225 0.222 0.220 0.219 0.217 0.216 0.215 0.213 0.210 0.209 0.209 0.209 0.210 0.219 0.253 0.293 0.309

0.049 0.057 0.061 0.071 0.079 0.096 0.113 0.130 0.164 0.225 0.324 0.366 0.390 0.563 0.753 0.953 1.000 -

0.009 0.025 0.034 0.053 0.067 0.092 0.112 0.127 0.141 0.151 0.166 0.172 0.175 0.196 0.213 0.229 0.232 -

0.003 0.003 0.003 0.004 0.005 0.006 0.007 0.008 0.009 0.009 0.011 0.011 0.011 0.013 0.015 0.016 0.017 -

0.006 0.019 0.026 0.041 0.052 0.077 0.105 0.139 0.210 0.302 0.380 0.402 0.412 0.464 0.496 0.517 0.521 -

0.993 0.981 0.974 0.959 0.947 0.922 0.893 0.859 0.785 0.686 0.597 0.571 0.558 0.491 0.446 0.413 0.407 -

5.51 4.46 4.02 3.31 2.94 2.41 2.09 1.90 1.75 1.66 1.53 1.49 1.46 1.34 1.24 1.17 1.15 -

105.96 82.56 75.21 59.27 51.48 39.47 32.55 28.47 25.32 23.28 20.62 19.60 19.21 16.47 14.41 12.98 12.64 -

7.60 6.04 5.40 4.33 3.75 2.86 2.23 1.74 1.17 0.83 0.67 0.64 0.62 0.56 0.53 0.52 0.51 -

0.29 0.27 0.26 0.25 0.25 0.25 0.25 0.26 0.28 0.32 0.36 0.38 0.39 0.43 0.47 0.51 0.51 -

5,58 4.52 4.08 3.37 2.99 2.45 2.10 1.87 1.62 1.35 1.08 1.00 0.96 0.75 0.62 0.53 0.51 0.51 0.52 0.56 0.71 0.91 1.00

5,58 4.52 4.08 3.37 2.99 2.45 2.10 1.87 1.62 1.35 1.08 1.00 0.96 0.75 0.62 0.53 0.51 0.51 0.52 0.56 0.71 0.91 1.00

Evaluation for negative compositions

Calculation of pinch branches and azeotropes outside of the composition space (for mole fractions > 1 and < 0) is possible, since the thermodynamic model for the description of the equilibrium relationship is nothing but a mathematical function 53 . Extending the calculations outside the composition space is interesting for two reasons: 1. Azeotropes with negative mole fractions can deviate inside of composition space at modified pressure levels, as illustrated in Section 5.4. These azeotropes can further be used as additional initial solutions for the calculation of pinch branches and pinch points.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2. Pinch points outside the composition space are also fixed points of the tray-to-tray equation system and may influence the trajectory of tray-to-tray profiles similar to pinch points inside the composition space 54 . The extension of the calculations outside the composition space can therefore provide valuable information for process synthesis and design. However, there are several reasons why care needs to be taken, if these solutions outside the composition space are to be taken into account. Although Beneke et al. 55 try to show the thermodynamic consistency of negative compositions by testing the Gibbs-Duhem relationship in combination with the NRTL model, the physical significance is at least questionable, since the results are purely theoretical and cannot be verified. In addition, these authors report local regions of discontinuity of the logarithm of the activity coefficient according to the NRTL model. Other thermodynamic models are even more restricted. While calculations with the Wilson model can at least slightly be extended outside the composition space, the UNIQUAC model is completely unsuited for such an evaluation, as indicated by Beneke et al. 55 and Felbab 54 . Consequently, even the possibility of performing calculations outside the composition space depends strongly on the choice of thermodynamical model. It is also interesting to note that regions of complex eigenvalues may exist outside of the composition space even for ideal mixtures, as previously described by Holland et al. 74 . The information about stable, unstable or saddle focus pinch points can be utilized within an analytical approach like the column profile maps 73,74 . Pinch-based methods, that try to utilize the complete information on fixed points, like the RBM 18,23 or the continuous distillation region method 21 need to be revised in order to include these complex fixed points appropriately. However, as pointed out by Köhler et al. 13 , the product specifications for the use of a pinchbased shortcut method should always reflect the assumption of an infinite number of trays (infinite column height) and therefore represent some kind of sharp split. In that case all relevant pinch points are located inside the physical composition space, at least for the recti46 ACS Paragon Plus Environment

Page 46 of 58

Page 47 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


fying and stripping section. Therefore inaccurate results of pinch-based shortcut methods for non-sharp splits, as e.g. reported by Felbab et al. 56 or Beneke et al. 55 , are not the result of an insufficiency of the pinch-based shortcut method, but rather of an inappropriate product specification.

C

Comparison of CPU time

In order to compare the results obtained for the implementation of the current method, R performed on a stand-alone PC with an Intel⃝ Core

TM

i5 3230 M (2.60 MHz) with those

reported in earlier articles we relate the CPU time to reported benchmarking results. A completely fair comparison would require to execute exactly the same code on the different machines, to take into account the different kinds of operations (floating point, if-then-else, assignments, exp, sqrt, etc.). Table 7: Benchmarking values for the machines used in different articles. Note: Since no MIPS values for the HP9000 J-2240 were available, the values for the slower HP9000 J-210 machine - P7200 are included. publication

machine

MIPS

speed-up rate

this work Felbab 54 Bonilla-Petriciolet et al. 49 Salomone and Espinosa 43 Harding and Floudas 48 Maier et al. 41 Wasylkiewicz et al. 36

R Intel⃝ Core i5 3230 M (2.60 GHz) TM R ⃝ Intel Core i5 430 M (2.53 GHz) TM R Intel⃝ Pentium M (1.73 GHz) TM R Intel⃝ Pentium II (333 MHz) HP9000 J-2240 - P8200 (240 MHz) Ultra 1/140 - UltraSparc (143 MHz) Alpha Station 250 - A21064A (266 MHz)

53263 29299 3377 555 336 215 198

1 1.75 15.77 95.94 149.75 247.73 254.14

TM

Since this is not possible and there is very little data reported for the old machines used by Wasylkiewicz et al. 36 , Maier et al. 41 and Harding and Floudas 48 , the results are related to the reported MIPS (millions of instructions per second) values from Dhrystone and Whetstone benchmarking tests 87,88 . The machines used for the computations in the different articles, the reported MIPS values and the speed-up rate in relation to the machine used for the current work are presented in Table 7. Even with these (ideal) speed-up ratios the current 47 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

method is capable of producing the correct results in at least competitive CPU time.

References 1. Humphrey, J. L.; Siebert, A. F. Separation technologies: an opportunity for energy savings. Chemical Engineering Progress 1992, 88, 32–41. 2. Sorensen, E. Distillation; Elsevier, 2014; pp 145–185. 3. Wankat, P. C. Separation process engineering: Includes mass transfer analysis, 3rd ed.; Prentice Hall: Upper Saddle River and NJ, 2012. 4. Harvey, L. D. D. Energy efficiency and the demand for energy services: Energy and the new reality; Earthscan: London and Washington and DC, 2010. 5. Koeijer, G. M. d.; Rivero, R. Entropy production and exergy loss in experimental distillation columns. Chemical Engineering Science 2003, 58, 1587–1597. 6. Widagdo, S.; Seider, W. D. Azeotropic distillation. AIChE Journal 1996, 42, 96–130. 7. Harding, S. T.; Floudas, C. A. European Symposium on Computer Aided Process Engineering - 11, 34th European Symposium of the Working Party on Computer Aided Process Engineering; Computer Aided Chemical Engineering; Elsevier, 2001; Vol. 9; pp 153–158. 8. Swietoslawski, W. Azeotropy and Polyazeotropy.; Pergamon Press: New York, 1963. 9. Kiva, V. N.; Hilmen, E. K.; Skogestad, S. Azeotropic phase equilibrium diagrams: a survey. Chemical Engineering Science 2003, 58, 1903–1953. 10. Rooks, R. E.; Julka, V.; Doherty, M. F.; Malone, M. F. Structure of distillation regions for multicomponent azeotropic mixtures. AIChE Journal 1998, 44, 1382–1391.


Page 48 of 58

Page 49 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


11. Thong, D. Y. C.; Jobson, M. Multicomponent homogeneous azeotropic distillation 1. Assessing Product Feasibility. Chemical Engineering Science 2001, 56, 4369–4391. 12. Ryll, O.; Blagov, S.; Hasse, H. ∞/∞-Analysis of homogeneous distillation processes. Chemical Engineering Science 2012, 84, 315–332. 13. Köhler, J.; Pöllmann, P.; Blass, E. A review on minimum energy calculations for ideal and nonideal distillations. Industrial & Engineering Chemistry Research 1995, 34, 1003– 1020. 14. van Dongen, D. B.; Doherty, M. F. Design and synthesis of homogeneous azeotropic distillations. 1. Problem formulation for a single column. Industrial & Engineering Chemistry Fundamentals 1985, 24, 454–463. 15. Wahnschafft, O. M.; Koehler, J. W.; Blass, E.; Westerberg, A. W. The product composition regions of single-feed azeotropic distillation-columns. Industrial & Engineering Chemistry Research 1992, 31, 2345–2362. 16. Julka, V.; Doherty, M. F. Geometric behavior and minimum flows for nonideal multicomponent distillation. Chemical Engineering Science 1990, 45, 1801–1822. 17. Köhler, J.; Aguirre, P.; Blass, E. Minimum reflux calculations for nonideal mixtures using the reversible distillation model. Chemical Engineering Science 1991, 46, 3007–3021. 18. Bausa, J.; Watzdorf, R. v.; Marquardt, W. Minimum energy demand for nonideal multicomponent distillations in complex columns. Computers & Chemical Engineering 1996, 20, 55–60. 19. Thong, D. Y. C.; Jobson, M. Multicomponent homogeneous azeotropic distillation 2. Column design. Chemical Engineering Science 2001, 56, 4393–4416. 20. Lucia, A.; Amale, A.; Taylor, R. Distillation pinch points and more. Computers & Chemical Engineering 2008, 32, 1342–1364. 49 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

21. Urdaneta, R. Y.; Bausa, J.; Br¨ uggemann, S.; Marquardt, W. Analysis and conceptual design of ternary heterogeneous azeotropic distillation processes. Industrial & Engineering Chemistry Research 2002, 41, 3849–3866. 22. Kraemer, K.; Harwardt, A.; Skiborowski, M.; Mitra, S.; Marquardt, W. Shortcut-based design of multicomponent heteroazeotropic distillation. Chemical Engineering Research and Design 2011, 89, 1168–1189. 23. Br¨ uggemann, S.; Marquardt, W. Shortcut methods for nonideal multicomponent distillation: 3. Extractive distillation columns. AIChE Journal 2004, 50, 1129–1149. 24. Avami, A. Conceptual Design of Double-Feed Reactive Distillation Columns. Chemical Engineering & Technology 2013, 36, 186–191. 25. Wahnschafft, O. M.; Lerudulier, J. P.; Westerberg, A. W. A problem decomposition approach for the synthesis of complex separation processes with recycles. Industrial & Engineering Chemistry Research 1993, 32, 1121–1141. 26. Thong, D. Y.-C.; Liu, G.; Jobson, M.; Smith, R. Synthesis of distillation sequences for separating multicomponent azeotropic mixtures. Chemical Engineering and Processing 2004, 43, 239–250. 27. Petlyuk, F. B.; Danilov, R.; Skouras, S.; Skogestad, S. Identification and analysis of possible splits for azeotropic mixtures. 2. Method for simple columns. Chemical Engineering Science 2012, 69, 159–169. 28. Petlyuk, F. B.; Danilov, R.; Skouras, S.; Skogestad, S. Identification and analysis of possible splits for azeotropic mixtures—1. Method for column sections. Chemical Engineering Science 2011, 66, 2512–2522. 29. Petlyuk, F.; Danilov, R.; Burger, J. A novel method for the search and identification


Page 50 of 58

Page 51 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of feasible splits of extractive distillations in ternary mixtures. Chemical Engineering Research and Design 2015, 99, 132–148. 30. Jaubert, J.-N.; Privat, R. Possible Existence of a Negative (Positive) Homogeneous Azeotrope When the Binary Mixture Exhibits Positive (Negative) Deviations from Ideal Solution Behavior (That is, When gE is Positive (Negative)). Industrial & Engineering Chemistry Research 2006, 45, 8217–8222. ˆ 31. Lecat, M. Sur l’azeAótropisme, particulie‘rement des syste‘mes binaires a‘ constituants chimiquement voisins: Azeotropism of binary systems of chemically similar constituents. C. R. Acad. Sci. Paris 1926, 880–882. 32. Wahnschafft, O. M.; Le Rudulier, J. P.; Blania, P.; Westerberg, A. W. Split: II. Authomated synthesis of hybrid liquid separation systems. Computers & Chemical Engineering 1992, 16, 305–312. 33. Missen, R. W. On Criteria for Occurrence of Azeotropes in Isothermal and Isobaric Binary Systems. The Canadian Journal of Chemical Engineering 2005, 83, 667–674. 34. Fidkowski, Z. T.; Malone, M. F.; Doherty, M. F. Computing azeotropes in multicomponent mixtures. Computers & Chemical Engineering 1993, 17, 1141–1155. 35. Eckert, E.; Kubicek, M. Computing heterogeneous azeotropes in multicomponent mixtures. Computers & Chemical Engineering 1997, 21, 347–350. 36. Wasylkiewicz, S. K.; Doherty, M. F.; Malone, M. F. Computing all homogeneous and heterogeneous azeotropes in multicomponent mixtures. Industrial & Engineering Chemistry Research 1999, 38, 4901–4912. 37. Tolsma, J. E.; Barton, P. I. Computation of heteroazeotropes. Part I: Theory. Chemical Engineering Science 2000, 55, 3817–3834.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

38. Tolsma, J. E.; Barton, P. I. Computation of heteroazeotropes. Part II: efficient calculation of changes in phase equilibrium structure. Chemical Engineering Science 2000, 55, 3835– 3853. 39. Qi, Z. W.; Sundmacher, K. Geometrically locating azeotropes in ternary systems. Industrial & Engineering Chemistry Research 2005, 44, 3709–3719. 40. Barbosa, D.; Doherty, M. F. The simple distillation of homogeneous reactive mixtures. Chemical Engineering Science 1988, 43, 541–550. 41. Maier, R. W.; Brennecke, J. F.; Stadtherr, M. A. Reliable computation of homogeneous azeotropes. AIChE Journal 1998, 44, 1745–1755. 42. Hua, J. Z.; Maier, R. W.; Tessier, S. R.; Brennecke, J. F.; Stadtherr, M. A. Interval analysis for thermodynamic calculations in process design: a novel and completely reliable approach. Fluid Phase Equilibria 1999, 160, 607–615. 43. Salomone, E.; Espinosa, J. Prediction of homogeneous azeotropes with interval analysis techniques exploiting topological considerations. Industrial & Engineering Chemistry Research 2001, 40, 1580–1588. 44. Zharov, W.; Serafimov, L. A. Physicochemical Fundamentals of Distilaltions and Rectifications (in Russian); Khimiya: Leningrad, 1975. 45. Pöllmann, P.; Bauer, M. H.; Blass, E. Investigation of vapour—liquid equilibrium of non-ideal multicomponent systems. Gas Separation & Purification 1996, 10, 225–241. 46. Harding, S. T.; Maranas, C. D.; McDonald, C. M.; Floudas, C. A. Locating All Homogeneous Azeotropes in Multicomponent Mixtures. Industrial & Engineering Chemistry Research 1997, 36, 160–178. 47. Maranas, C. D.; Floudas, C. A. Finding all solutions of nonlinearly constrained systems of equations. Journal of Global Optimization 1995, 7, 143–182. 52 ACS Paragon Plus Environment

Page 52 of 58

Page 53 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


48. Harding, S. T.; Floudas, C. A. Locating all heterogeneous and reactive azeotropes in multicomponent mixtures. Industrial & Engineering Chemistry Research 2000, 39, 1576– 1595. 49. Bonilla-Petriciolet, A.; Iglesias-Silva, G. A.; Hall, K. R. Calculation of homogeneous azeotropes in reactive and non-reactive mixtures using a stochastic optimization approach. Fluid Phase Equilibria 2009, 281, 22–31. 50. Fidkowski, Z. T.; Malone, M. F.; Doherty, M. F. Nonideal multicomponent distillation: Use of bifurcation theory for design. AIChE Journal 1991, 37, 1761–1779. 51. Julka, V.; Doherty, M. F. Geometric nonlinear analysis of multicomponent nonideal distillation: a simple computer-aided design procedure. Chemical Engineering Science 1993, 48, 1367–1391. 52. Keller, H. B. Numerical solution of bifurcation and nonlinear eigenvalue problems. Applications of Bifurcation Theory. 1977; pp 359–384. 53. Pöllmann, P.; Blass, E. Best products of homogeneous azeotropic distillations. Gas Separation & Purification 1994, 8, 194–229. 54. Felbab, N. An Efficient Method of Constructing Pinch Point Curves and Locating Azeotropes in Nonideal Distillation Systems. Industrial & Engineering Chemistry Research 2012, 51, 7035–7055. 55. Beneke, D. A.; Kim, S. B.; Linninger, A. A. Pinch Point Calculations and Its Implications on Robust Distillation Design. Chinese Journal of Chemical Engineering 2011, 19, 911– 925. 56. Felbab, N.; Hildebrandt, D.; Glasser, D. A new method of locating all pinch points in nonideal distillation systems, and its application to pinch point loci and distillation boundaries. Computers & Chemical Engineering 2011, 35, 1072–1087. 53 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

57. Bausa, J.; Marquardt, W. Quick and reliable phase stability test in VLLE flash calculations by homotopy continuation. Computers & Chemical Engineering 2000, 24, 2447– 2456. 58. Westerberg, A. W.; Wahnschafft, O. M. Synthesis of distillation-based separation systems. Advances in chemical engineering 1996, 23, 63–170. 59. Laroche, L. Homogeneous azeotropic distillation: entrainer selection; California Institute of Technology, PhD thesis: California, 1991. 60. Hilmen, E. K. Separation of Azeotropic Mixtures: Tools for analysis and studies on batch distillation operation; Norwegian University of Science and Technology, PhD-thesis: Norway, 2000. 61. Petlyuk, F. B. Distillation theory and its application to optimal design of separation units; Cambridge University Press: Cambridge, 2004. 62. Br¨ uggemann, S. Rapid screening of conceptual design alternatives for distillation processes; Fortschrittberichte VDI, VDI Verlag, Reihe 3, Nr.841: D¨ usseldorf, 2005. 63. Skiborowski, M.; Harwardt, A.; Marquardt, W. Efficient optimization-based design for the separation of heterogeneous azeotropic mixtures. Computers & Chemical Engineering 2015, 72, 34–51. 64. Treybal, R. E. Liquid Extraction, 2nd ed.; McGraw-Hill: New-York, 1963. 65. Guddat, J.; Guerra Vazques, F.; Jongen, H. T. Parametric optimization: Singularities, pathfollowing and jumbs; Teubner: Stuttgart, 1990. 66. Br¨ uggemann, S.; Oldenburg, J.; Zhang, P.; Marquardt, W. Robust dynamic simulation of three-phase reactive batch distillation columns. Industrial & Engineering Chemistry Research 2004, 43, 3672–3684.


Page 54 of 58

Page 55 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


67. Kangas, J.; Malinen, I.; Tanskanen, J. Modified bounded homotopies in the solving of phase stability problems for liquid–liquid phase-splitting calculations. Industrial & Engineering Chemistry Research 2011, 50, 7003–7018. 68. Zhang, H.; Bonilla-Petriciolet, A.; Rangaiah, G. P. A review on global optimization methods for phase equilibrium modeling and calculations. The Open Thermodynamics Journal 2011, 5, 71–92. 69. Seydel, R.; Hlavacek, V. Role of continuation in engineering analysis. Chemical Engineering Science 1987, 42, 1281–1295. 70. Fowler, M. UML distilled: A brief guide to the standard object modeling language, 3rd ed.; The Addison-Wesley object technology series; Addison-Wesley: Boston, 2004. 71. Deuflhard, P.; Hohmann, A. Eine algorithmisch orientierte Einf¨ uhrung, 3rd ed.; Numerische Mathematik; De Gruyter: Berlin, 2002; Vol. 1. 72. Choi, S. H.; Harney, D. A.; Book, N. L. A robust path tracking algorithm for homotopy continuation. Computers & Chemical Engineering 1996, 20, 647–655. 73. Tapp, M.; Holland, S. T.; Hildebrandt, D.; Glasser, D. Column profile maps. 1. Derivation and interpretation. Industrial & Engineering Chemistry Research 2004, 43, 364– 374. 74. Holland, S. T.; Tapp, M.; Hildebrandt, D.; Glasser, D. Column Profile Maps. 2. Singular Points and Phase Diagram Behaviour in Ideal and Nonideal Systems. Industrial & Engineering Chemistry Research 2004, 43, 3590–3603. 75. Doherty, M. F. Properties of liquid vapor composition surfaces for multicomponent mixtures with constant latent-heat. Chemical Engineering Science 1985, 40, 1979–1980. 76. Bausa, J.; Watzdorf, R. v.; Marquardt, W. Shortcut methods for nonideal multicomponent distillation: 1. Simple columns. AIChE Journal 1998, 44, 2181–2198. 55 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

77. Serafimov, L. A. Thermodynamic Topological Analysis and the Separation of Multicomponent Polyazeotropic Mixtures. Theoretical Foundations of Chemical Engineering 1987, 21, 44–54. 78. Wasylkiewicz, S. K.; Kobylka, L. C.; Castillo, F. J. L. Pressure Sensitivity Analysis of Azeotropes. Industrial & Engineering Chemistry Research 2003, 42, 207–213. 79. Krolikowski, L. J. Determination of distillation regions for non-ideal ternary mixtures. AIChE Journal 2006, 52, 532–544. 80. Thong, D. Y. C.; Jobson, M. Multicomponent homogeneous azeotropic distillation 3. Column sequence synthesis. Chemical Engineering Science 2001, 56, 4417–4432. 81. Kraemer, K.; Kossack, S.; Marquardt, W. Efficient optimization-based design of distillation processes for homogeneous azeotropic mixtures. Industrial & Engineering Chemistry Research 2009, 48, 6749–6764. 82. Br¨ uggemann, S.; Marquardt, W. Conceptual design of distillation processes for mixtures with distillation boundaries: I. Computational assessment of split feasibility. AIChE Journal 2011, 57, 1526–1539. 83. Br¨ uggemann, S.; Marquardt, W. Conceptual design of distillation processes for mixtures with distillation boundaries. II. Optimization of recycle policies. AIChE Journal 2011, 57, 1540–1556. 84. Knight, J. R.; Doherty, M. In Foundations of computer-aided process design; Siirola, J. J., Grossmann, I. E., Stephanopoulos, G., Eds.; Elsevier: Amsterdam and New York, 1990; p 419. 85. Wasylkiewicz, S. K.; Kobylka, L. C.; Castillo, F. J. L. Optimal design of complex azeotropic distillation columns. Chemical Engineering Journal 2000, 79, 219–227.


Page 56 of 58

Page 57 of 58

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


86. Doherty, M. F.; Perkins, J. D. Dynamics of distillation processes. 3. Topological structure of ternary residue curve maps. Chemical Engineering Science 1979, 34, 1401–1414. 87. Longbottom, R. Historic Data and Older Benchmarks. 2016;

http://www.

roylongbottom.org.uk/mips.htm. 88. Notebookcheck, Benchmarking results for Intel core i5 processors. 2016; http://www. notebookcheck.net/.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 17: Table of Content (TOC) graphic.


Page 58 of 58

A Unifying Approach for the Calculation of ... - ACS Publications

Recommend Documents