Simulation and Design Methods for Multiphase Multistream Heat Exchangers*

11th 11th IFAC IFAC Symposium Symposium on on Dynamics Dynamics and and Control Control of of Process Systems, including Biosystems 11th IFAC IFACSystems, Symposium on Dynamics Dynamics and Control Control of of 11th Symposium on and Process including Biosystems Available online at www.sciencedirect.com June 2016. Trondheim, Norway Process including Biosystems Process Systems, including Biosystems June 6-8, 6-8,Systems, 2016. NTNU, NTNU, Trondheim, Norway June June 6-8, 6-8, 2016. 2016. NTNU, NTNU, Trondheim, Trondheim, Norway Norway

ScienceDirect

IFAC-PapersOnLine 49-7 (2016) 839–844

Simulation and Design Methods for Simulation Simulation and and Design Design Methods Methods for for ⋆⋆ Multiphase Multistream Heat Exchangers Multiphase Multistream Heat Exchangers Multiphase Multistream Heat Exchangers ⋆

* ** Harry Harry A. A. J. J. Watson Watson ** Paul Paul I. I. Barton Barton ** ** * Paul I. Barton ** Harry A. J. Watson Harry A. J. Watson Paul I. Barton * * Massachusetts Institute of Technology , Cambridge, MA 02139 USA Massachusetts Institute of Technology , Cambridge, MA 02139 USA * * Massachusetts Institute of ,, Cambridge, (e-mail: [email protected]). Massachusetts Institute of Technology Technology Cambridge, MA MA 02139 02139 USA USA (e-mail: [email protected]). ** (e-mail: [email protected]). ** Massachusetts Institute of Technology , Cambridge, MA 02139 USA (e-mail: [email protected]). Massachusetts Institute of Technology , Cambridge, MA 02139 USA ** ** Massachusetts Institute Institute of Technology , Cambridge, MA 02139 USA (e-mail: [email protected]) Massachusetts of Technology , Cambridge, MA 02139 USA (e-mail: [email protected]) (e-mail: (e-mail: [email protected]) [email protected])

Abstract: Multistream Multistream heat heat exchangers exchangers (MHEXs) (MHEXs) are are found found in in many many energy energy intensive and Abstract: intensive and Abstract: Multistream heat exchangers (MHEXs) are found in many energy intensive and industrially relevant cryogenic processes. However, design and optimization of such processes is Abstract: Multistream heat exchangers (MHEXs) are found in many energy intensive and industrially relevant cryogenic processes. However, design and optimization of such processes is industrially relevant processes. However, design and of processes rendered difficult difficult by cryogenic the inability inability to MHEXs robustly, as most most current current flowsheet-level industrially relevant cryogenic processes. However, design and optimization optimization of such such processes is is rendered by the to simulate simulate MHEXs robustly, as flowsheet-level rendered difficult by to MHEXs robustly, current flowsheet-level MHEX models models solve onlyinability an energy energy balance with with no constraint constraint onmost second law feasibility. feasibility. The rendered difficult by the the inability to simulate simulate MHEXs robustly, as as most current flowsheet-level MHEX solve only an balance no on second law The MHEX models solve only an energy balance with no constraint on second law feasibility. The new model described herein combines an extension of the classical pinch analysis algorithm MHEX models solve only an energy balance with no constraint on second law feasibility. The new model described herein combines an extension of the classical pinch analysis algorithm new herein combines an the pinch algorithm with model explicitdescribed dependence on the the heat exchange exchange area to toof formulate nonsmooth equation system new model described herein combines an extension extension offormulate the classical classical pinch analysis analysis algorithm with explicit dependence on heat area aa nonsmooth equation system with explicit dependence on the heat exchange area to formulate a nonsmooth equation system which can be solved for up to three unknown variables in an MHEX. The resulting model is with on to thethree heat exchange to formulate a nonsmooth equationmodel system whichexplicit can bedependence solved for up unknown area variables in an MHEX. The resulting is which can be solved for up to three unknown variables in an MHEX. The resulting model is further augmented to simulate realistic thermodynamic phenomena, such as phase change and which be solvedtofor up to realistic three unknown variables phenomena, in an MHEX. The model is furthercan augmented simulate thermodynamic such as resulting phase change and further augmented to simulate simulate realistic thermodynamic phenomena, suchequations. as phase phase change change and and equilibrium, which are are also naturally naturally described by additional additional nonsmooth equations. further augmented to realistic thermodynamic phenomena, such as equilibrium, which also described by nonsmooth equilibrium, equilibrium, which which are are also also naturally naturally described described by by additional additional nonsmooth nonsmooth equations. equations. © 2016, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. Keywords: Heat exchangers, Process models, Process simulators, Numerical methods, Energy Energy Keywords: Heat exchangers, Process models, Process simulators, Numerical methods, Keywords: Heat Heat exchangers, exchangers, Process Process models, models, Process Process simulators, simulators, Numerical Numerical methods, methods, Energy Energy Keywords: 1. INTRODUCTION INTRODUCTION 1. 1. INTRODUCTION 1. INTRODUCTION Multistream heat exchangers in in cryogenic cryogenic processes processes are are Multistream heat exchangers Multistream heat exchangers exchangers in cryogenic cryogenic processes are notoriously difficult difficult to model, model, simulate, simulate, andprocesses design. Such Such Multistream heat in are notoriously to and design. notoriously difficult to model, model, simulate, andintensive, design. Such Such processes are are by nature nature extremely energy intensive, and notoriously difficult to simulate, and design. processes by extremely energy and processes are by nature extremely energy intensive, and therefore stand to benefit greatly from accurate process processes are by nature extremely energy intensive, and therefore stand to benefit greatly from accurate process therefore stand to benefit benefit greatly from accurate process optimization; however, without robust and flexible process models therefore stand to greatly from accurate optimization; however, without robust and flexible models optimization; however, without robust and flexible models for heat exchange unit operations, simulation and optioptimization; however, robust and flexible for heat exchange unitwithout operations, simulation andmodels optifor heat exchange unit operations, simulation and optimization of these processes cannot be performed. Taking for heat exchange unit operations, simulation andTaking optimization of these processes cannot be performed. mization of these these gas processes cannot be performed. performed. Taking liquefied natural natural gas (LNG)cannot production as an an example, example, mization of processes be Taking liquefied (LNG) production as liquefied natural gas (LNG) production as example, the streams streams in the thegas process MHEXs are usually usually both mulliquefied natural (LNG) production as an anboth example, the in process MHEXs are multhe streams in the process MHEXs are usually both multicomponent and multiphase, which creates an even more the streams inand themultiphase, process MHEXs usually both more multicomponent which are creates an even ticomponent and multiphase, multiphase, which creates an an evenphases more challenging simulation problem, especially if the the phases ticomponent and which creates even more challenging simulation problem, especially if challenging simulation problem, especially if the phases traversed in the exchangers are not known ahead of time. challenging simulation problem, especially if the phases traversed in the exchangers are not known ahead of time. traversed in the are not ahead traversed the exchangers exchangers aresoftware not known known ahead of ofintime. time. The use use of ofinprocess process simulation is common common the The simulation software is in the The use of process simulation software is common in the literature involving processes with multistream heat exThe use of process simulation software is common in the literature involving processes with multistream heat exliterature involving with multistream heat exchangers. Commercial Commercial simulators employ proprietary modliterature involving processes processes with multistream heatmodexchangers. simulators employ proprietary changers. Commercial simulators employ proprietary models which which Commercial generally permit permit solving for aaproprietary single unknown unknown changers. simulators employ models generally solving for single els which generally permit solving for a single unknown variable, afforded by the energy balance, typically taken els which generally permit solving for a single unknown variable, afforded by the energy balance, typically taken variable, afforded by the the energy energy balance, typically taken as the the one oneafforded of the the exchanger exchanger outlet balance, temperatures. However, variable, by typically taken as of outlet temperatures. However, as the one of the exchanger outlet temperatures. However, an example, the first author’s experience using the as the of the exchanger outlet temperatures. However, an one example, the first author’s experience using the R ○ as an example, the first author’s experience using the R ○ (Aspen Technology MHEATX block from Aspen Plus as an example, author’s experience using the (Aspen Technology MHEATX block the fromfirst Aspen Plus○ R R (Aspen ○ Technology MHEATX block from Aspen Plus Inc., 2014) suggests there are rigorous checks in (Aspen Technology MHEATX from Aspen Plus Inc., 2014) block suggests there are no no rigorous checks in place place Inc., 2014) suggests there are no rigorous checks place to avoid heat exchange between two streams at very Inc., 2014) suggests there are no rigorous checks insimilar place to avoid heat exchange between two streams at veryin similar to avoid heat exchange between two streams at very similar temperatures, or prevent temperature crossovers. to avoid heat exchange between two streams at very similar temperatures, or prevent temperature crossovers. temperatures, or crossovers. temperatures, or prevent prevent temperature temperature crossovers. involves A A more more rigorous rigorous approach approach to to modeling modeling MHEXs MHEXs involves A more rigorous approach to modeling MHEXs the use of a superstructure concept (Hasan et 2009). A rigorous approach to modeling MHEXs involves themore use of a superstructure concept (Hasan et al., al.,involves 2009). the use of a superstructure concept (Hasan et al., The major disadvantage of this methodology is that simuthe use of adisadvantage superstructure concept (Hasan et al., 2009). 2009). The major of this methodology is that simuThe major disadvantage of this methodology is that simulating the MHEX involves the solution of a very complex The major disadvantage of this methodology is that simulating the MHEX involves the solution of a very complex lating the MHEX involves the solution of a very complex mixed-integer nonlinear program (MINLP) model, which lating the MHEX involves the solution of a very complex mixed-integer nonlinear program (MINLP) model, which mixed-integer nonlinear (MINLP) is to mixed-integer nonlinear program program (MINLP) model, model, which which is extremely extremely challenging challenging to find find globally. globally. is extremely challenging to find globally. is extremely challenging to find globally.   The The authors authors are are  for this research. are  The authors are forThe thisauthors research. for this research. for this research.

grateful grateful grateful grateful

to to to to

Statoil Statoil Statoil Statoil

for for for for

providing providing providing providing

financial financial financial financial

support support support support

Another Another method method for for robust robust modeling modeling of of MHEXs MHEXs borrows borrows Another method for robust modeling of borrows heavily from pinch analysis and the analysis of Another method foranalysis robust and modeling of MHEXs MHEXs borrows heavily from pinch the analysis of composite composite heavily from pinch analysis and the analysis of composite curves. Kamath et al. (2012) showed how to create heavily from pinch analysis and the analysis of composite curves. Kamath et al. (2012) showed how to create a a fully fully curves. Kamath al. how a equation-oriented (EO) modelshowed for MHEXs MHEXs bycreate considering curves. Kamath et et(EO) al. (2012) (2012) showed how to toby create a fully fully equation-oriented model for considering equation-oriented for by the as heat (HEN) equation-oriented (EO) model for MHEXs MHEXsnetwork by considering considering the unit unit operation operation(EO) as aamodel heat exchanger exchanger network (HEN) the operation as heat (HEN) that requires no Applied to Rethe unit operation as aa utilities. heat exchanger exchanger network (HEN) that unit requires no external external utilities. Appliednetwork to the the Poly Poly Rethat requires no external utilities. Applied to the Poly Refrigerant Integrated Cycle Operations (PRICO) process, that requires no external utilities. Applied to the Poly Refrigerant Integrated Cycle Operations (PRICO) process, frigerant Integrated Cycle (PRICO) process, the formulation formulation results in aaOperations moderately-sized mathematfrigerant Integrated Cycle Operations (PRICO) process, the results in moderately-sized mathematthe results in mathematical program complementarity constraints (MPCC) the formulation results in aa moderately-sized moderately-sized mathematical formulation program with with complementarity constraints (MPCC) ical program with complementarity constraints (MPCC) (3426 equations using Soave-Redlich-Kwong thermodyical program with complementarity constraints (MPCC) (3426 equations using Soave-Redlich-Kwong thermody(3426 equations using Soave-Redlich-Kwong thermodynamics) which requires completing a complicated (3426 equations using Soave-Redlich-Kwong thermodynamics) which requires completing a complicated initialinitialnamics) which completing aa complicated ization procedure procedure to obtain obtain suitable initial guess guessinitialfrom namics) which requires requires completing complicated initialization to aa suitable initial from ization obtain aa suitable guess from which solve problem, well solution method ization procedure to obtain as suitable initial guess from which to toprocedure solve the the to problem, as well as as aainitial solution method which to solve the problem, as well as a solution method suitable for MPCCs. which to solve the problem, as well as a solution method suitable for MPCCs. suitable for suitable for MPCCs. MPCCs. approaches require the solution of Both Both of of these these rigorous rigorous approaches require the solution of ofoptimization these rigorous rigorous approaches require the solutionof ofa aBoth problem to the feasibility Both these approaches require solution a hard hardof optimization problem to ensure ensure the the feasibility of of a aaMHEX. hard optimization problem to ensure the feasibility of It is also notable that the MHEX modeling literhard optimization problem to ensure the feasibility of a MHEX. It is also notable that the MHEX modeling liter-a MHEX. It also that the MHEX ature mention of of exMHEX. It is ismakes also notable notable that the dependence MHEX modeling modeling literature rarely rarely makes mention of the the dependence of heat heatliterexature makes of of exchange area the performance of the Rather, ature rarely makes mention of the the dependence of heat heat exchangerarely area on on the mention performance of dependence the operation. operation. Rather, change area on the performance of the operation. Rather, the size is usually only calculated following determination change area on the performance of the operation. Rather, the size is usually only calculated following determination the size is only calculated following determination of stream states, if In this the sizeoutput is usually usually only calculated following determination of the the output stream states, if at at all. all. In contrast, contrast, this artiartiof the output stream states, if at all. In contrast, article describes a model and solution procedure for MHEXs of output astream if at all.procedure In contrast, this articlethe describes modelstates, and solution for this MHEXs cle describes a model and solution procedure for MHEXs which avoids returning non-physical solutions, doesn’t rely cle describes a model and solution procedure MHEXs which avoids returning non-physical solutions, for doesn’t rely which avoids returning non-physical solutions, doesn’t rely on approximations or solving a hard optimization probwhich avoids returning non-physical solutions, doesn’t rely on approximations or solving a hard optimization probon approximations or solving a hard optimization problem, incorporates information about the available heat on approximations or solving a hard optimization problem, incorporates information about the available heat lem, incorporates information the heat exchange area into procedure, and lem, incorporates information about the available available heat exchange area directly directly into the theabout procedure, and handles handles exchange area directly into the procedure, and handles phase changes without needing to fix conditions a priori. exchange area without directly needing into thetoprocedure, anda handles phase changes fix conditions priori. phase changes changes without without needing needing to to fix fix conditions conditions a a priori. priori. phase 2. 2. PRELIMINARIES PRELIMINARIES 2. 2. PRELIMINARIES PRELIMINARIES This paper considers This paper considers the the modeling modeling of of countercurrent countercurrent This paper considers the modeling ofindexed countercurrent MHEXs in which a set of hot streams, by This paper considers the modeling countercurrent MHEXs in which a set of hot streams,of indexed by a a set set MHEXs in which a set of hot streams, indexed by a 𝐻𝐻, exchange heat with a set of cold streams, indexed by MHEXs in which a set of hot streams, indexed by a set set 𝐻𝐻, exchange heat with a set of cold streams, indexed by exchange aa set streams, indexed by a𝐻𝐻, 𝐶𝐶. hot stream ∈ 𝐻𝐻 aa constant molar 𝐻𝐻, exchange heat with set𝑖𝑖𝑖𝑖 of of cold streams, indexed by a set set 𝐶𝐶. Each Eachheat hotwith stream ∈ cold 𝐻𝐻 has has constant molar in aaheat set 𝐶𝐶. Each hot stream 𝑖𝑖 ∈ 𝐻𝐻 has a constant molar in , capacity flowrate 𝐹𝐹 , enters at temperature 𝑇𝑇 Cp ,i𝑖𝑖, ∈ set capacity 𝐶𝐶. Each flowrate hot stream 𝐻𝐻 has a constant molar i heat 𝐹𝐹C enters at temperature 𝑇𝑇 ,i p iin in , out in out heatexits capacity flowrate 𝐹𝐹 𝐹𝐹C , enters at temperature 𝑇𝑇 , ,i p i out in out heat capacity flowrate , enters at temperature 𝑇𝑇 , with 𝑇𝑇 ≥ 𝑇𝑇 . Similarly, and at temperature 𝑇𝑇 ,i , with 𝑇𝑇i ≥ 𝑇𝑇i i , . Similarly, and exits at temperature C𝑇𝑇piiout out , with 𝑇𝑇iiin in ≥ 𝑇𝑇iiout out . Similarly, and exits at temperature 𝑇𝑇 i and exits at temperature 𝑇𝑇i , with 𝑇𝑇i ≥ 𝑇𝑇i . Similarly,

Copyright 2016 IFAC IFAC 839 Hosting by Elsevier Ltd. All rights reserved. 2405-8963 © IFAC (International Federation of Automatic Control) Copyright © 2016, 2016 839 Copyright 2016 IFAC 839 Peer review© of International Federation of Automatic Copyright ©under 2016 responsibility IFAC 839Control. 10.1016/j.ifacol.2016.07.294

IFAC DYCOPS-CAB, 2016 840 Harry A.J. Watson et al. / IFAC-PapersOnLine 49-7 (2016) 839–844 June 6-8, 2016. NTNU, Trondheim, Norway

each cold stream 𝑗𝑗 ∈ 𝐶𝐶 has a constant molar heat capacity flowrate 𝑓𝑓Cp ,j , enters at temperature 𝑡𝑡in j , and exits at in out temperature 𝑡𝑡out , with 𝑡𝑡 ≤ 𝑡𝑡 . The overall heat transfer j j j coefficient of the MHEX is denoted by 𝑈𝑈 and the available heat transfer area is denoted by 𝐴𝐴. The product 𝑈𝑈 𝑈𝑈 is often used as the metric for the physical ability of a MHEX to exchange heat to avoid the need for accurate estimation of the value of 𝑈𝑈 in early-stage design. The energy balance around this unit operation is given by: ∑︁ ∑︁ 𝑓𝑓Cp ,j (𝑡𝑡out − 𝑡𝑡in (1) 𝐹𝐹Cp ,i (𝑇𝑇iin − 𝑇𝑇iout ) − j j ) = 0. i∈H

j∈C

However, it is not sufficient for an MHEX to satisfy only an energy balance – all heat transfer must be both feasible in the sense of the second law and also possible in terms of the available area. While it is straightforward to write equations for these conditions in the two-stream heat exchanger case, such constraints do not immediately generalize to the multistream case. However, two equations, and the necessary algorithms for their evaluation, can be formulated to enforce these requirements, as detailed in Watson et al. (2015) and summarized in the following section. The equations in this MHEX model are nonsmooth and require specialized algorithms to solve. Facchinei et al. (2014) developed a linear program (LP) based Newton method that works well for such systems. At each step, this method solves the following LP to find the next iterate: min 𝛾𝛾 γ,u ⃦ ⃦ ⃦ ⃦2 s.t. ⃦f (uk ) + G(uk )(u − uk )⃦∞ ≤ 𝛾𝛾 ⃦f (uk )⃦∞ , (2) ⃦ ⃦ ⃦ ⃦ ⃦(u − uk )⃦ ≤ 𝛾𝛾 ⃦f (uk )⃦ , ∞ ∞ u ∈ 𝑋𝑋𝑋

where G(uk ) is a generalized derivative of f at uk , 𝛾𝛾 is a variable introduced to drive convergence to the solution, and 𝑋𝑋 is a set of known bounds on the solution. The u part of the solution is then used as the next iterate for the algorithm. In practice, the infinity norms are replaced by their equivalent linear reformulations. When f is piecewisecontinuously differentiable (𝒫𝒫𝒫𝒫 1 ) in the sense of Scholtes (2012), and G(uk ) is an element of the B-subdifferential (Qi, 1993) of f at uk , the method can achieve local Q-quadratic convergence to the solution. The automatic method of Khan and Barton (2015b) can be used to evaluate B-subdifferential elements for 𝒫𝒫𝒫𝒫 1 functions using a vector forward mode of automatic differentiation. Additionally, the LP-Newton method will not necessarily fail if G(uk ) is singular, unlike methods based on linear equation solving. This makes the method particularly well suited for the application of MHEX simulation.

pinch analysis, heat transfer is considered feasible if the hot and cold composite curves corresponding to the streams in the MHEX are separated by at least some minimum temperature difference, Δ𝑇𝑇min . With the additional goal of maximizing the energy transferred in the MHEX, a feasible solution corresponds to the vertical separation between the hot and cold composite curves being exactly Δ𝑇𝑇min . This is equivalent to the more readily solved problem of reducing the horizontal distance between the composite curves at the pinch point to zero after applying a temperature shift of Δ𝑇𝑇min to the cold curve. Adding nonphysical extensions of the curves which extend to the maximum and minimum temperatures existing in the heat exchanger allows this distance to be evaluated everywhere in this range, even if both of the true composite curves were not defined at every pinch candidate temperature. With these modifications, the hot and cold composite curve enthalpy values corresponding to each pinch candidate are defined using the following expressions: ∑︁ [︀ p 𝐹𝐹Cp ,i max{0, 𝑇𝑇 p − 𝑇𝑇iout }−max{0, 𝑇𝑇 p − 𝑇𝑇iin } = 𝐸𝐸𝐸𝐸𝐸𝐸H i∈H

]︀ − max{0, 𝑇𝑇 min − 𝑇𝑇 p } + max{0, 𝑇𝑇 p − 𝑇𝑇 max } ,

(3)

∀𝑝𝑝 ∈ 𝑃𝑃𝑃 ∑︁ [︀ 𝐸𝐸𝐸𝐸𝐸𝐸Cp = 𝑓𝑓Cp ,j max{0, (𝑇𝑇 p − Δ𝑇𝑇min ) − 𝑡𝑡in j } j∈C

(4) − max{0, (𝑇𝑇 p − Δ𝑇𝑇min ) − 𝑡𝑡out j } p max + max{0, (𝑇𝑇 − Δ𝑇𝑇min ) − 𝑡𝑡 } ]︀ min p − (𝑇𝑇 − Δ𝑇𝑇min )} , ∀𝑝𝑝 ∈ 𝑃𝑃𝑃 − max{0, 𝑡𝑡

where 𝑇𝑇 max/min is the maximum/minimum hot stream temperature, 𝑡𝑡max/min is the maximum/minimum cold stream temperature, and 𝑃𝑃 = 𝐻𝐻 ∪ 𝐶𝐶 is the set of pinch point candidates with temperatures, 𝑇𝑇 p , defined by: {︂ in 𝑇𝑇i , ∀𝑝𝑝 = 𝑖𝑖 ∈ 𝐻𝐻𝐻 p 𝑇𝑇 = (5) + Δ𝑇𝑇 , 𝑡𝑡in min ∀𝑝𝑝 = 𝑗𝑗 ∈ 𝐶𝐶𝐶 j

In (3) and (4), the last two terms correspond to the nonphysical extensions of the curves.

p , ∀𝑝𝑝 ∈ 𝑃𝑃 , the miniAfter calculating 𝐸𝐸𝐸𝐸𝐸𝐸Cp − 𝐸𝐸𝐸𝐸𝐸𝐸H mum horizontal distance between the extended composite curves is found by searching through the finite set 𝑃𝑃 : p min{𝐸𝐸𝐸𝐸𝐸𝐸Cp − 𝐸𝐸𝐸𝐸𝐸𝐸H } = 0. (6) p∈P

The expression inside the min statement of (6) is negative at any pinch temperature where the hot composite curve lies to the right of the shifted cold composite curve, and positive where it lies to the left. 3.2 Formulation of MHEX area constraint

3. A NEW MODEL FOR MULTISTREAM HEAT EXCHANGERS This section summarizes the novel model formulation of Watson et al. (2015). 3.1 Formulation of MHEX feasible heat transfer constraint As shown in Kamath et al. (2012), a multistream heat exchanger can be viewed as a special case of a HEN where the external utilities are not present. Therefore, as in classical 840

An algorithm for calculating the heat exchange area for a MHEX on each iteration of an equation-solving procedure is described in this section. Let 𝐷𝐷 be the index set for the points at which the composite curves are nondifferentiable, as well as their endpoints. For 𝑑𝑑 ∈ 𝐷𝐷, let 𝑄𝑄d denote the enthalpy value at this point, and construct the triple (𝑄𝑄d , 𝑇𝑇 d , 𝑡𝑡d ), where 𝑇𝑇 d is the temperature on the hot composite curve at 𝑄𝑄d , and 𝑡𝑡d is the temperature on the cold composite curve at 𝑄𝑄d . If all such triples are sorted in

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Harry A.J. Watson et al. / IFAC-PapersOnLine 49-7 (2016) 839–844

order of nondecreasing 𝑄𝑄d value, then an adjacent pair of triples in the list demarcates an interval of the composite curves in which part of the MHEX can be modeled as a two-stream heat exchanger. The total MHEX heat transfer area can then be found by summing over all such intervals: ∑︁ Δ𝑄𝑄d 𝑈𝑈 𝑈𝑈 = , (7) d Δ𝑇𝑇LM d∈D d̸=|D|

where Δ𝑄𝑄d = 𝑄𝑄d+1 − 𝑄𝑄d is the width of enthalpy interval d 𝑑𝑑, |𝐷𝐷| = 2(|𝐻𝐻| + |𝐶𝐶|), and Δ𝑇𝑇LM is a modified version of the log-mean temperature difference across this same enthalpy interval which is defined later in this section. In standard practice, heat transfer area is calculated after heat integration is performed, and the the integrated composite curves are used to define the enthalpy intervals. However, by including (7) in the system of equations being solved simultaneously, an additional unknown for simulation is freed as the area can be specified in the problem input. However, the full sorted list of (𝑄𝑄d , 𝑇𝑇 d , 𝑡𝑡d ) data must then be calculated at each iteration, starting from an incomplete list of just the inlet and outlet temperatures which are arranged in an arbitrary order. Sorting this list into nondecreasing order based on enthalpy value sets up the intervals for (7) correctly. If the sort is performed using bubble sort, then the only operations involve taking the max or min of two values. If the for loops in bubble sort run over the entire length of the list to ensure the same number of operations are performed for a given input size (regardless of how well-sorted the input is), then this sorting operation is a composite 𝒫𝒫𝒫𝒫 1 function which maps the unsorted input to the sorted output. The missing temperatures in the triples (e.g. 𝑇𝑇 d ) are calculated as follows given the value of 𝑄𝑄d : ∑︁ )︀ (︀ 𝐹𝐹Cp ,i max{0, 𝑇𝑇 d −𝑇𝑇iout }−max{0, 𝑇𝑇 d −𝑇𝑇iin } = 0. 𝑄𝑄d − i∈H

(8) An analogous equation implicitly defines unknown cold curve temperatures. The solution value of these equations can be determined by searching and interpolating over the piecewise affine segments of the relevant composite curve. Calculating the correct generalized derivatives of the unknown temperature at the solution of (8) requires application of an implicit function result for 𝒫𝒫𝒫𝒫 1 functions from Khan and Barton (2015a). Finally, the standard definition of the log-mean temperature difference must be modified so that evaluating the function never results in undefined behavior, as follows (based on Zavala-R´ıo et al. (2005)): ⎧ 1 ⎪ if Δ𝑇𝑇 d = Δ𝑇𝑇 d+1 , ⎨ (Δ𝑇𝑇 d +Δ𝑇𝑇 d+1 ) 2 d Δ𝑇𝑇LM = (9) Δ𝑇𝑇 d+1 −Δ𝑇𝑇 d ⎪ ⎩ otherwise, d+1 d ln(Δ𝑇𝑇 )−ln(Δ𝑇𝑇 ) d where Δ𝑇𝑇 = max{Δ𝑇𝑇min , 𝑇𝑇 d − 𝑡𝑡d } is the temperature difference at the start of enthalpy interval 𝑑𝑑, and Δ𝑇𝑇 d+1 = max{Δ𝑇𝑇min , 𝑇𝑇 d+1 − 𝑡𝑡d+1 } is the temperature difference at the end of the interval. The if statement in the definition of the log-mean temperature difference is necessary to make the calculation defined for all possible inputs Δ𝑇𝑇 d > 0 and Δ𝑇𝑇 d+1 > 0. Fortunately, since this function is 841

841

continuously differentiable on the positive quadrant of R2 (Zavala-R´ıo et al., 2005), the if statement in (9) does not introduce complications for calculating correct derivatives. In summary, the set of equations describing the multistream heat exchanger now consists of (1), (6), and (7). This is a nonsmooth equation system involving three equations in three unknowns. In particular, the left-hand side of each equation is a composite 𝒫𝒫𝒫𝒫 1 function, as each can be expressed as a finite composition of continuouslydifferentiable functions and 𝒫𝒫𝒫𝒫 1 functions such as the binary max and min functions. Including bounds on the temperature variables is also recommended to avoid nonrealistic solutions where hot streams have gained heat as if they were cold streams and vice-versa. 4. MODELING PHASE CHANGES IN MULTISTREAM HEAT EXCHANGERS In the context of heat integration and MHEX simulation problems, the challenge associated with a stream changing phase is modeling the change in its heat capacity flowrate. Pinch analysis techniques assume every stream has a constant heat capacity flowrate, and therefore an affine temperature-enthalpy relationship, but na¨ıvely assuming this for a liquefying multicomponent refrigerant or natural gas stream would introduce significant error. Therefore, the temperature-enthalpy behavior of such streams must be approximated by a series of affine segments. A good choice is to approximate each of the three individual phase regions (superheated vapor, two-phase, and subcooled liquid) with one (or possibly more) affine segments. However, this is not straightforward in the general case where the process stream inlet/outlet temperatures, pressures, and compositions are possibly variables in the simulation. Kamath et al. (2012) propose a disjunctive model in which the physical streams in the process are subdivided into substreams corresponding to superheated (sup), two-phase (2p), and subcooled (sub) regions. The heat integration calculations for the MHEX are then performed using these substreams instead of the physical process streams. However, instead of using either binary variables or complementarity constraints to assign the phase substream inlet/outlet temperatures as in Kamath et al. (2012), the in/out in/out in/out values of 𝑇𝑇sup , 𝑇𝑇sub , and 𝑇𝑇2p can be expressed as continuous nonsmooth functions of the process stream temperature (𝑇𝑇 IN/OUT ), the dew point (𝑇𝑇DP ), and the bubble point (𝑇𝑇BP ) as follows: in/out = max{𝑇𝑇DP , 𝑇𝑇 IN/OUT }, 𝑇𝑇sup in/out 𝑇𝑇2p in/out 𝑇𝑇sub

= mid{𝑇𝑇DP , 𝑇𝑇

IN/OUT

IN/OUT

, 𝑇𝑇BP },

(10) (11)

= min{𝑇𝑇 , 𝑇𝑇BP }, (12) 3 where the function mid : R → R maps to the median of its three arguments and is also a 𝒫𝒫𝒫𝒫 1 function.

The possible appearance and disappearance of phases from iteration to iteration also causes problems for vapor-liquid equilibrium calculations. Consider a typical steady-state flash operation, in which a feed stream with molar flowrate 𝐹𝐹 with 𝑛𝑛c components at molar composition z separates into a liquid stream with molar flowrate 𝐿𝐿 at molar composition x and a vapor stream with molar flowrate 𝑉𝑉 at molar composition y. In the classical approach to flash

IFAC DYCOPS-CAB, 2016 842 Harry A.J. Watson et al. / IFAC-PapersOnLine 49-7 (2016) 839–844 June 6-8, 2016. NTNU, Trondheim, Norway

calculations, it is assumed that both liquid and vapor outlet streams exist and are in equilibrium, and the material balance and equilibrium relationships (expressed in terms of 𝑘𝑘i , the distribution coefficient for component 𝑖𝑖) are combined into the well-known Rachford-Rice equation (Rachford Jr. and Rice, 1952). When the heat duty of the flash, 𝑄𝑄flash , is specified (as will be the case in the MHEX simulations in this work), an energy balance is also needed. These two equations can then be solved via a Newton-type method for the vaporized fraction of the feed, VF , and the flash temperature, from which all other relevant quantities in the formulation can be recovered. However, under conditions where only one outlet stream exists, the equilibrium conditions cannot be satisfied. A new nonsmooth model which accounts for the appearance and disappearance of phases while performing flash calculations is given by solving the following equation instead of the Rachford-Rice equation: }︃ {︃ nc ∑︁ 𝑉𝑉 𝑉𝑉 𝑧𝑧i (𝑘𝑘i − 1) , − 1, − = 0. (13) mid 𝐹𝐹 𝐹𝐹 1 + VF (𝑘𝑘i − 1) i=1

hot two-phase substream, as well as the substreams in the other regions of the phase diagram. These nonsmooth model elements can be combined to simulate multiphase MHEXs in complex processes. In many cases, splitting a physical stream into just three substreams will be insufficient to capture its true temperatureenthalpy relationship. Therefore, as described in Kamath et al. (2012), the superheated, subcooled, and two-phase substreams are further discretized into 𝑛𝑛sup , 𝑛𝑛sub , and 𝑛𝑛2p affine segments, respectively, to improve the estimation of the nonlinear behavior. The substreams are subdivided into segments of equal enthalpy difference, and their inlet/outlet temperatures are therefore implicitly defined by the solution of energy balances. Given the temperature and heat load for each segment, a constant heat capacity flowrate can be determined by dividing the heat load of the segment by the difference between its inlet and outlet temperature. Note that the use of (14) and (15) will never allow this quotient to become undefined. These heat capacity flowrates are then used in evaluating (1), (6), and (7). The size of this MHEX model is given by:

Here, the three arguments in the mid function correspond to finding an all liquid outlet, an all vapor outlet and a twophase outlet, respectively. Note that the third term is just the negative of the standard Rachford-Rice residual. The working mechanism of the equation is as follows. When the outlet is all vapor, 𝑉𝑉 𝑉𝑉𝑉 = 1 and the Rachford-Rice expression is positive, so that here, the first term is equal to 1, the second term is equal to zero and the third term is negative. Thus, the mid expression picks the second term and evaluates to zero, satisfying (13) with 𝑉𝑉 𝑉𝑉𝑉 = 1. The arguments for the other phases are analogous.

𝑛𝑛var = 3 + 2𝑛𝑛str + 𝑛𝑛str [(𝑛𝑛sup −1) + (𝑛𝑛sub −1) + 2(𝑛𝑛2p −1)] , (16) where 𝑛𝑛var is the number of variables/equations needed to model the MHEX, and 𝑛𝑛str is the number of physical process streams that enter the MHEX in the flowsheet. The first term of (16) accounts for the base MHEX model consisting of (1), (6), and (7), the second term accounts for calculating the dew and bubble points of each stream involved in the heat exchanger, and the last term accounts for all the substream temperatures (and vapor fractions in the two-phase region) which must be calculated.

A final provision can be added in order to simulate streams that undergo very small (or no) temperature changes, as is the case with pure component phase changes or when phase substreams disappear based on the current value of the model variables. The inclusion of nonsmooth functions in this work allows this behavior to be modeled more precisely than in approaches which rely on smoothing approximations, while avoiding the need for binary variables or disjunctions. For instance, if it is likely that a cold substream will undergo a very small or no temperature change in the two-phase region during iteration, (11) can be replaced with the following expressions based on the discussion in Kamath et al. (2012): {︁ IN 𝑡𝑡in , 𝑇𝑇BP }, 2p =min mid{𝑇𝑇DP , 𝑇𝑇 IN mid{𝑇𝑇DP , 𝑇𝑇 , 𝑇𝑇BP } + mid{𝑇𝑇DP , 𝑇𝑇 OUT , 𝑇𝑇BP } − 𝜀𝜀 }︁ , 2 (14) {︁ OUT 𝑡𝑡out , 𝑇𝑇BP }, 2p =max mid{𝑇𝑇DP , 𝑇𝑇 IN mid{𝑇𝑇DP , 𝑇𝑇 , 𝑇𝑇BP } + mid{𝑇𝑇DP , 𝑇𝑇 OUT , 𝑇𝑇BP } + 𝜀𝜀 }︁ , 2 (15) where 𝜀𝜀 is a user-defined fictitious temperature change. Conventionally this approach is considered undesirable; however, in the specific context of nonsmooth equations, nonsmooth equation solvers will not be affected by the usual ill-conditioning caused by making such a parameter too small. Analogous equations can also be written for a

Despite the reduction in model complexity as compared to other approaches, the problem remains difficult to solve. Improved robustness of the LP-Newton algorithm was observed by replacing all the infinity norms in (2) with the 1norm, and changing the expression{︁on the right-hand side ⃦ ⃦ ⃦ ⃦2 }︁ of the first constraint to 𝛾𝛾 max ⃦f (uk )⃦1 , ⃦f (uk )⃦1 . The change of norm was motivated by the observation that the step size calculated by (2) at each iteration is bounded by unity if the largest residual value corresponds to an equation with a zero row in the generalized derivative. This occurs often with (6), as any temperature variables in the simulation will only influence the minimum distance between the composite curves over a limited range of values. The addition of the max term to the first ⃦ constraint ⃦ allows the method to take larger steps when ⃦f (uk )⃦1 > 1, which helps mitigate sensitivity to problem scaling since this method does not share the invariance to affine scaling exhibited by classical Newton methods.

842

Finally, it is important that the initial guess provided be as near to the solution as possible to aid convergence, which necessitates a robust initialization procedure. In the procedure developed by the authors, only guesses for the three MHEX model variables and the bubble/dew point temperatures of each stream are required from the user. The bubble and dew point estimates are then refined first by solving their defining equations. Initial guesses for the remaining temperatures are obtained by assuming a linear relationship between temperature and enthalpy in each phase, and then improved by solving each of the energy

IFAC DYCOPS-CAB, 2016 June 6-8, 2016. NTNU, Trondheim, Norway Harry A.J. Watson et al. / IFAC-PapersOnLine 49-7 (2016) 839–844

balances independently to generate a better estimate. Example 1 shows the method of this paper applied to simulate the PRICO process under the assumption of ideal thermodynamic behavior and Raoult’s Law. Example 1. Figure 1 shows the PRICO (Maher and Sudduth, 1975) process for producing liquefied natural gas. The PRICO process is a single-stage mixed refrigerant (MR) process. The MR stream serves as both a hot stream, high-pressure refrigerant (HPR), and as a cold stream, low-pressure refrigerant (LPR) by means of intermediate expansion and compression operations. Tables 1 and 2 give

843

Case I. In this case, the pressures and compositions in OUT the flowsheet are held fixed. Let 𝑢𝑢1 ≡ 𝑇𝑇HPR , 𝑢𝑢2 ≡ 𝑡𝑡OUT LPR , and 𝑢𝑢3 ≡ 𝑈𝑈 𝑈𝑈 be the unknown variables afforded by the base MHEX model consisting of (1), (6), and (7). The throttle valve outlet temperature, 𝑢𝑢4 ≡ 𝑡𝑡IN LPR and vapor V IN fraction, 𝑢𝑢5 ≡ F LPR , are given by performing a fixed pressure/enthalpy flash calculation around the valve. The remainder of the variable set is comprised of the unknown temperatures and vapor fractions given by the solution of the energy balances around each affine segment. Using CPLEX v12.5 (IBM, 2015) as the LP solver, the simulation converges to the solution with 𝑢𝑢*1 = 120.78 K, 𝑢𝑢*2 = 286.51 K, 𝑢𝑢*3 = 18.83 MW/K, 𝑢𝑢*4 = 116.61 K, and 𝑢𝑢*5 = 0.11 with ‖f (u* )‖∞ < 10−9 after 56 iterations (0.63 seconds) starting from 𝑢𝑢01 = 118.15 K, 𝑢𝑢02 = 275 K, 𝑢𝑢03 = 20 MW/K, 𝑢𝑢04 = 116.95 K, 𝑢𝑢05 = 0.10, and the rest of the initial guess calculated by the initialization subroutine. Figure 2 shows the composite curves at the solution. In spite of the piecewise-affine approximations, the composite curves show much of the true curvature, particularly around the natural gas stream bubble point (195.59 K). 320 300 280

Fig. 1. Flowsheet of the idealized PRICO liquefaction process for natural gas. the data for the streams involved in the MHEX for the PRICO process under three sets of simulation conditions. In each case, 𝑛𝑛sup = 3, 𝑛𝑛2p = 5, and 𝑛𝑛sub = 3 for all three process streams entering the heat exchanger. From (16), the MHEX is modeled with 45 equations and variables.

240 220 200 180 160 140 120 100

0

10000

20000

30000

40000

50000

60000

70000

80000

90000

Enthalpy (kW)

Table 1. Natural gas stream data.

Fig. 2. Composite curves for the MHEX in the PRICO process simulated under the conditions of Case I.

Property Value Flowrate (kmol/s) 1.00 Pressure (MPa) 5.500 Inlet temperature (K) 298.15 Outlet temperature (K) 118.15 Composition (mol %) N2 CH4 C2 H 6 C3 H 8 n-C4 H10 1.0 95.6 3.1 0.2 0.1

Case II. In this case, the composition of the refrigerant mixture is allowed to vary while the 𝑈𝑈 𝑈𝑈 value of the MHEX is fixed. Let 𝑢𝑢1 be the mole fraction of n-butane OUT in the refrigerant, 𝑢𝑢2 ≡ 𝑇𝑇LPR and 𝑢𝑢3 ≡ Δ𝑇𝑇min with 𝑢𝑢4 and 𝑢𝑢5 as before. The simulation converges to the solution with 𝑢𝑢*1 = 0.2595, 𝑢𝑢*2 = 285.93 K, 𝑢𝑢*3 = 0.99 K, 𝑢𝑢*4 = 114.324 K, and 𝑢𝑢*5 = 0.094 with ‖f (u* )‖∞ < 10−9 after 38 iterations (0.44 seconds) starting from 𝑢𝑢01 = 0.2562, 𝑢𝑢02 = 286.5 K, 𝑢𝑢03 = 1.2 K, 𝑢𝑢04 = 116.95 K, 𝑢𝑢05 = 0.10, and the rest of the initial guess calculated by the initialization subroutine. Figure 3 shows the composite curves at the solution, which are pinched together more closely than in Case I due to the larger available heat transfer area.

Table 2. Refrigerant stream and MHEX data. Property Flowrate (kmol/s) High pressure level (MPa) Low pressure level (MPa) HPR inlet temp. (K) HPR outlet temp. (K) LPR inlet temp. (K) LPR outlet temp. (K) Composition (mol %) Nitrogen Methane Ethane Propane n-Butane U A (MW/K) ∆Tmin (K)

Temperature (K)

260

Case I 3.47 3.695 0.170 303.15 u1 u4 u2

Case II 3.47 3.695 0.170 303.15 118.15 u4 u2

Case III 3.47 3.695 u1 303.15 118.15 u4 u2

15.32 17.79 40.85 0.41 25.62 u3 1.2

15.32 17.79 40.85 73.93 − 100u1 100u1 20 u3

15.32 17.79 40.85 0.41 25.62 10 u3

Case III. In this case, the pressure of the MR stream is allowed to vary while the 𝑈𝑈 𝑈𝑈 value of the MHEX is fixed. Let 𝑢𝑢1 be the discharge pressure of the throttle valve, 𝑢𝑢2 ≡ 𝑡𝑡OUT LPR and 𝑢𝑢3 ≡ Δ𝑇𝑇min with 𝑢𝑢4 and 𝑢𝑢5 as before. The simulation converges to the solution with 𝑢𝑢*1 = 0.117 MPa, 𝑢𝑢*2 = 286.51 K, 𝑢𝑢*3 = 2.77 K, 𝑢𝑢*4 = 112.67 K, and 𝑢𝑢*5 = 0.125 after 114 iterations (1.31 seconds) starting from 𝑢𝑢01 = 0.17 MPa, 𝑢𝑢02 = 286.5 K, 𝑢𝑢03 = 1.2 K, 𝑢𝑢04 = 116.95 K, 𝑢𝑢05 = 0.10, and the rest of the initial guess calculated by the initialization subroutine. Figure 4 shows 843

IFAC DYCOPS-CAB, 2016 844 Harry A.J. Watson et al. / IFAC-PapersOnLine 49-7 (2016) 839–844 June 6-8, 2016. NTNU, Trondheim, Norway

property calculation package; however, performing such a simulation for processes larger than PRICO is extremely challenging, and future work involves investigating both a modular approach and improvements to the nonsmooth algorithm for nonideal flash calculations.

320 300 280

Temperature (K)

260 240 220

ACKNOWLEDGEMENTS

200 180

The authors wish to thank Kamil Khan for many helpful discussions regarding the mathematical background of this work. The authors are also grateful to Truls Gundersen at the Norwegian University of Science and Technology (NTNU) for his input on the project.

160 140 120 100

0

10000

20000

30000

40000

50000

60000

70000

80000

90000

Enthalpy (kW)

Fig. 3. Composite curves for the MHEX in the PRICO process simulated under the conditions of Case II. the composite curves at the solution. Satisfying the low 𝑈𝑈 𝑈𝑈 value specification requires the pressure at the throttle valve exit to approach atmospheric pressure, leading to greater temperature change across the valve and more separation of the composite curves. 320 300 280

Temperature (K)

260 240 220 200 180 160 140 120 100

0

10000

20000

30000

40000

50000

60000

70000

80000

90000

Enthalpy (kW)

Fig. 4. Composite curves for the MHEX in the PRICO process simulated under the conditions of Case III. In each case, the model is successfully able to simulate the PRICO process under ideal thermodynamics. Cases II and III are challenging simulation problems since many of the equations in the model are sensitive to the composition and pressure of the refrigerant stream. However, the new simulation strategy converges to the solution in each case without requiring significant computational time or effort. 5. CONCLUSION A method for simulating multistream heat exchangers both with and without phase changes has been presented which differs substantially from those found in the literature to date. The final formulation is a compact nonsmooth model which is solved with equation-solving methods, in contrast to those which require solving a difficult optimization problem. Accordingly, the model is significantly less complex and allows for realistic simulation of process flowsheets involving multiphase MHEXs outside of an optimization framework. The method can be extended to incorporate cubic equations of state using the EO formulation from Kamath et al. (2010) and a physical 844

REFERENCES Aspen Technology Inc. (2014). Aspen Plus v8.4. Aspen Technology Inc., Cambridge, MA. Facchinei, F., Fischer, A., and Herrich, M. (2014). An LPNewton method: nonsmooth equations, KKT systems, and nonisolated solutions. Mathematical Programming, 146, 1–36. Hasan, M.M.F., Karimi, I.A., Alfadala, H.E., and Grootjans, H. (2009). Operational modeling of multistream heat exchangers with phase changes. AIChE Journal, 55(1), 150–171. IBM (2015). IBM ILOG CPLEX v12.5. http://www01.ibm.com/software/commerce/optimization/cplexoptimizer/index.html. Kamath, R.S., Biegler, L.T., and Grossmann, I.E. (2012). Modeling multistream heat exchangers with and without phase changes for simultaneous optimization and heat integration. AIChE Journal, 58(1), 190–204. Kamath, R.S., Biegler, L.T., and Grossmann, I.E. (2010). An equation-oriented approach for handling thermodynamics based on cubic equation of state in process optimization. Computers & Chemical Engineering, 34(12), 2085–2096. Khan, K.A. and Barton, P.I. (2015a). Generalized derivatives for hybrid systems. Submitted. Khan, K.A. and Barton, P.I. (2015b). A vector forward mode of automatic differentiation for generalized derivative evaluation. Optimization Methods & Software, 30(6), 1185–1212. Maher, J. and Sudduth, J. (1975). Method and apparatus for liquefying gases. US Patent 3,914,949. Qi, L. (1993). Convergence analysis of some algorithms for solving nonsmooth equations. Mathematics of Operations Research, 18(1), 227–244. Rachford Jr., H.H. and Rice, J.D. (1952). Procedure for use of electronic digital computers in calculating flash vaporization hydrocarbon equilibrium. Journal of Petroleum Technology, 4(10), 327–328. Scholtes, S. (2012). Introduction to Piecewise Differentiable Equations. SpringerBriefs in Optimization, New York, NY. Watson, H.A.J., Khan, K.A., and Barton, P.I. (2015). Multistream heat exchanger modeling and design. AIChE Journal, 61(10), 3390–3403. Zavala-R´ıo, A., Femat, R., and Santiesteban-Cos, R. (2005). An analytical study of the logarithmic mean temperature. Revista Mexicana de Ingenier´ıa Qu´ımca, 4, 201–212.