State Function‐Based Flash Specifications for Open Systems in the Absence or Presence of Chemical Reactions

Diverse engineering ﬁelds request ﬂash calculations like isothermal ﬂash, isenthalpic ﬂash, and isentropic ﬂash. They can be cast as minimization of a thermodynamic state function and solved by Michelsen’s Q-function approach. Flash calculations for open systems, i


Introduction
Phase equilibrium calculation is essential to process simulation in many engineering disciplines and scientific areas, such as chemical engineering, petroleum engineering, mechanical engineering, and geochemistry.The calculation determines the partitioning of species between equilibrium phases, and it is often coupled with mass-and energybalances for process analysis and design.In a broad sense, phase equilibrium calculation also covers reaction equilibrium, although the term usually refers to the situations with phase equilibrium only.
Phase equilibrium calculation covers a wide range of topics, like flash, saturation point, phase diagrams, critical loci, and even composition gradients in a petroleum reservoir.It is impossible to make a full list but a classification of these calculations, like the one given by Michelsen in 1993 1 , helps us to differentiate their different characteristics.
Michelsen showed that for a mixture with a given feed composition, its single-stage equilibrium calculation requires just two specifications, which is somewhat similar to the Duhem theorem.He classified the equilibrium calculation problems into three types 1 : flash specifications or specifications leading to thermodynamic state function minimiza-

Accepted Article
This article is protected by copyright.All rights reserved.
tion, phase fraction specifications, and other specifications.The first type refers to isothermal (TP) flash, isenthalpic (PH) flash, isentropic (PS) flash, etc.The second type is represented by saturation point calculation, but the specified phase fraction can also be a value between 0 and 1.The last type includes somewhat exotic specifications like (P ,V ) and (T , H ) and indirect specifications like critical point calculation.Among the three types, only the first type always corresponds to a unique minimum of a thermodynamic state function, and thus can be solved in a robust manner with the help of stability analysis 2 .Among all the flash specifications, TP flash is the most widely used and has well-developed algorithms especially for two-phase calculation.The modern algorithms for isothermal flash are characterized by utilizing stability analysis for robustness and a second-order method for efficiency.Since its solution corresponds to a minimum in the Gibbs energy, it is advantageous to design the algorithm as Gibbs energy minimization.As an example of the modern flash algorithms, Michelsen's algorithm 3 consists of two important steps: stability analysis 2 and phase split 4 .In the phase split part, mole numbers in a certain phase are selected as independent variables in a second-order minimization of the Gibbs energy.A noteworthy variation of the classical TP flash algorithms is the reduced variables methods or reduction methods 5;6;7;8;9;10;11 .The methods consist in using a different set of independent variables (reduced variables) to reduce significantly the size of the Jacobian matrix in certain situations.Most of these algorithms are not formulated as minimization, but Michelsen et al. (2013)  12 showed that this is possible for a special reduced variables formulation 11 .According to some recent studies 12;13 , the speed gain of the reduced variables methods is rather modest over the classical methods.The findings of Petitfrere and Nichita 14 go in the same direction and show that reduced methods are more efficient only for systems with many components.Another markedly different and promising approach for phase split calculation is the RAND-based methods 15;16;17;18 .These methods are designed to calculate chemical and phase equilibrium simultaneously.Even for situations with only phase equilibrium, the RAND-based methods also show advantages for multiphase flash due to their smaller and structured matrices.These methods can also utilize the Gibbs energy minimization to guide the convergence.
For other state function-based flash specifications, Michelsen in 1999 19 provided a systematic approach to handle these problems.These flash specifications include PH, PS, VT, UV and VS.Michelsen pointed out that all these specifications, just as the TP specification, can be formulated as minimization of a certain thermodynamic state func-

Accepted Article
This article is protected by copyright.All rights reserved.
tion, thus corresponding to a unique solution.He proposed different Q -functions for these problems i .Since most of these Q -functions correspond to a saddle point, he suggested using nested loops for a robust solution.In addition, he showed how to establish a common Jacobian to solve all the specifications using an efficient Newton method.
Although his discussion was only for two-phase flash examples with the conventional selection of the independent variables (mole numbers), the approach is not limited by the number of phases or the selection of independent variables.For instance, Paterson et al. illustrated how to use his approach for multiphase PH flash 20 and how to solve other flash specifications under the RAND framework 21 .
In contrast to Michelsen's Q -function approach, another way to handle other flash specifications is to perform unconstrained minimization of the state function by use of its natural variables as the independent variables in the flash formulation.This is actually practiced in the classical TP flash where (T , P , n)-based thermodynamics is used.
For (T , P , n)-based thermodynamics, the thermodynamic properties are evaluated at given (T , P , n) as we normally do for an equation of state (EoS).In the classical TP flash, T and P are fixed and n is the set of independent variables in the flash formulation.The resulting procedure is an unconstrained Gibbs energy minimization.Similarly, the other five specifications become unconstrained minimization if their flash formulations are based on their respective natural variables.Except for VT flash using a (T ,V , n)-based formulation, the others are seldom investigated in the literature.A reason for this is that estimating properties at (T ,V , n) is natural for most equations of state, whereas estimating properties, especially derivatives, at (P , H , n), (P , S, n), (U ,V , n), and (V , S, n) is possible but cumbersome.Nevertheless, the approach was tried by Brantferger et al. 22 and Sun et al. 23 for PH flash and Castier 24 for UV flash.Recently, Paterson et al. generalized the approach for all specifications under the classical framework 17 and the RAND framework 21;25 .
The flash specifications discussed above correspond to extrema in thermodynamic state functions in non-open systems.It should be noted that thermodynamic state functions or potentials can also be introduced for open systems through the Legendre transform 26;27;28 .In these new potentials ii , mole numbers are replaced by their conjugate chemical potentials as the new natural variables.Depending on how many mole numbers are left in the natural i Michelsen used objective functions denoted by Q in his other algorithms as well, e.g., his multiphase flash algorithm and his chemical equilibrium algorithm, we use Q -functions specifically for the functions that he introduced in his 1999 paper 19 .
ii In Michelsen's discussion on state function-based flash specifications, the state functions are actually various thermodynamic potentials.

Accepted Article
This article is protected by copyright.All rights reserved.
variables for the new potentials, the corresponding system can be semi-open to completely open.We propose here to utilize these new potentials to extend the flash-type calculations to open systems.The new equilibrium calculations are still called "flash" although the original flash concept refers to the separation in a flash drum with a given feed.
Nevertheless, like the classical flash with a given feed, the new flash for open systems is also under the first type of Michelsen's classification because all the new flash calculations can be formulated as thermodynamic state-function minimization.
Equilibrium in an open system is in fact quite common in thermodynamics-related subjects.A well-known example in statistical mechanics is the grand canonical ensemble where chemical potentials (µ), total volume (V ) and temperature (T ) are held constant.The macroscopic system corresponding to the grand canonical ensemble needs to minimize the grand potential, also called the Landau potential 29 .Alberty gave a comprehensive discussion on using the Legendre transform in chemical thermodynamics 27;28 and proposed several non-classical thermodynamic potentials useful for many scientific and engineering problems.Some of Alberty's examples are for open systems, such as reactions in a vessel connected to a reservoir of a given component by a semi-permeable membrane and biochemical reactions held at constant chemical potentials.Another area that often requires flash in open systems is geochemistry, where equilibrium calculations with complex speciation and mineral precipitations are often performed at constant partial pressure of carbon dioxide or at constant pH 30 .
In this work, we illustrate how to perform these extended "flash" calculations for open systems in the absence or presence of chemical reactions by extending Michelsen's Q -function approach.We first briefly review the Legendre

Accepted Article
This article is protected by copyright.All rights reserved.

Legendre Transform in Thermodynamics
The Legendre transform is an important mathematical manipulation used in classical mechanics, statistical mechanics and thermodynamics.The Legendre transform g = g (x , y ) of a strictly convex and continuous function f = f (x, y ) is defined as 31 It can be shown that x = ∂f ∂x y and x = − ∂g ∂x y . If f = f (x, y ) has a unique minimum, so does g = g (x , y ) .
In thermodynamics, the Legendre transform is used to generate more convenient state functions.In the Gibbs formalism, the fundamental relation of thermodynamics is written for the internal energy U as a function of natural extensive variables only, including entropy S , volume V and number of moles of the system n.The integral and differential forms for U (S,V , n) are given in Table 1.Using the Legendre transform, we can create from U (S,V , n) three new thermodynamic state functions, i.e., A(T ,V , n), H (S, P , n) and G (T , P , n), with new sets of natural variables.
Each Legendre transform replaces an extensive variable with an intensive one.A single transform is needed to obtain A and H , and two transforms to obtain G .These new functions retain the same amount of information as U (S,V , n) 26 , e.g., the equilibrium condition can be formulated as the minima in these functions under different constraints.In addition to the Gibbs formalism, the less common Massieu formalism is useful to our later discussion on state functionbased flash calculations 32 .The Massieu formalism starts with entropy S (U ,V , n) as a function of U , V and n, with its differential and integral expressions also given in Table 1.From S , we can create A/T with one Legendre transform, and G /T with two.We here introduce a modification of the Massieu formalism by assuming the original expression for S is in (H , P , n) instead of (U ,V , n).Under the modified Massieu formalism, S is converted to G /T and A/T after one and two transforms, respectively.It should be noted that for the transform to A/T , we actually replace an intensive variable (P ) with an extensive one (V /T ), which is the opposite to what we usually perform for thermodynamic functions.Table 1 lists the aforementioned potentials where the mole numbers of the system are kept constant.In section "General Construction Procedure for the Q -function", we utilize these potentials and create new potentials for open systems based on them.For a more comprehensive discussion of other useful potentials introduced by the Legendre transform,

Accepted Article
This article is protected by copyright.All rights reserved.
the readers are referred to Alberty's papers 27;28  and (U ,V ).Except for (V , S ), the other five flash specifications have important practical applications.All these flash specifications correspond to the minimization of a thermodynamic state function.Michelsen proposed two approaches to handle the problem.He showed that each flash specification corresponds to a Q -function.The solution to the flash problem is at the stationary point of the Q -function, and the stationary point is usually a saddle point, except for (T , P ) specification in (T , P )-based thermodynamics and (T ,V ) and (T , P ) specifications in (T ,V )-based thermodynamics.For the specifications corresponding to a saddle problem, Michelsen suggested using a nested loop, e.g.solving a (T , P )-flash in the inner loop and a maximization of the Q -function in the outer loop.This Q -function approach provides robustness, but it is not very efficient.Michelsen suggested using a Newton approach for an efficient solution, and he showed that all the six flash specifications share a common Jacobian matrix.
The Q -functions suggested by Michelsen can be introduced through the method of Lagrange multipliers.Let us take the (U ,V )-flash as an example.If (T , P )-based thermodynamics is used, the minimization formulation is the following: min subject to

Accepted Article
This article is protected by copyright.All rights reserved.
The problem consists in a minimization with respect to (N C × N P + 2) variables.The linear constraints (Equation 3c) are used to remove N C dependent phase mole numbers, leaving just two non-linear constraints: Equations 3a and 3b.
A Lagrangian function can be constructed to convert the problem into finding the stationary points of the following function: At the solution, the two Lagrange multipliers should be given by The above derivatives can be obtained more easily from Massieu's formulation.Also note that Equations 5a and 5b are valid only if the derivatives, as well as the values of P and T , are taken at the solution.However, if we use P and T as normal variables to express the two multipliers and substitute them into Equation 4, we obtain the expression for the Q -function: The Q -function is not equivalent to the Lagrangian function, since Equations 5a and 5b are valid only at the solution.Nonetheless, we can show that the solution is at the saddle point of this Q -function, meaning that if we perform a (P , T )-flash in the inner loop and solve for T and P in the outer loop, the outer loop is a maximization

Accepted Article
This article is protected by copyright.All rights reserved.
of Q .It is also interesting to notice that in Equation 4, we have essentially applied the Legendre transform first to change from (−S ) to (G /T ), by the terms (-λ U U ) and (-λ V V ), and from (G /T ) to Q , by the terms λ U U spec and λ V V spec .
Hence, the nature of the problem is not changed through this manipulation.Both transforms applied here are based on Massieu's formulation of thermodynamic functions.For the transform from (−S ) to (G /T ), it involves The transform from (G /T ) to Q is the reverse process, but it is a rigorous transform only at the solution.In other words, the Q -function has the same value as the entropy at the solution.The other Q -functions in Michelsen's work, no matter if (T , P )-based or (T ,V )-based, can be obtained in a similar manner.

| General Construction Procedure for the Q -function
We can generalize the construction of a Q -function through the method of Lagrange multipliers for a flash problem in which a given thermodynamic function Υ needs to be minimized.For simplicity, we first discuss (T , P )-based thermodynamics where the solution of the (T , P )-flash can be formulated as an unconstrained minimization and forms the inner loop of the final Q -function-based approach.The state function Υ has a different set of natural variables from the set (T , P , n) for the Gibbs energy, but Υ and G are linked through a series of Legendre transforms.
We denote the different natural variables of Υ by q k , and Υ and G should be linked by Legendre transforms to ensure both functions contain essentially the same information: In the above equation, ∂Υ ∂q k are conjugate variables of q k , i.e., some variables in the set (T , P , n).N var is the number of variables being transformed.The variables q k appear as non-linear constraints in the minimization using the method of Lagrange multipliers:

Accepted Article
This article is protected by copyright.All rights reserved.
At the solution, we have The Q -function is, then, constructed as: Actually, Equation 9becomes slightly more complex if q k is an intensive variable, e.g., chemical potential (detailed discussion can be found in section "Open System Flash without Reactions").In such a case, the conjugate variable for , is an extensive variable, e.g., mole number.For any thermodynamic model, the thermodynamic properties are always calculated for each phase.Therefore, we need to distinguish ∂Υ m ∂q m k and q l m for different phases and we add a superscript m to both Υ and q k to denote this.For a system with N P phases, this leads to constraints on q m k : Since all q m k will be equal to q spec k at the solution, we can denote it by q k , their corresponding terms in the Lagrange function can be combined in the construction of the Q -function, yielding If q k is an intensive variable, its conjugate, ∂Υ ∂q k , is a homogeneous Euler function of degree one, such that Finally, Equation 12becomes

Accepted Article
This article is protected by copyright.All rights reserved.
The term ∂Υ ∂q k in Equation 14represents the total extensive property for the system.We can see that the final expression for the Q -function 11b is actually preserved.The above complication can be illustrated by the example of a two-phase (T , P )-flash using (T ,V )-based thermodynamics.The constraints on pressure need to be written for both gas and liquid phases: The conjugate variables for P v and P l are V v and V l , respectively.The Q -function is constructed as follows: In summary, for a flash problem corresponding to minimization of a thermodynamic state function Υ, if we choose to use thermodynamics based on the natural variables of another state function Ξ, e.g., G for (T , P )-based thermodynamics and A for (T ,V )-based thermodynamics, in the solution process, and Υ and Ξ are related by Legendre transforms: Then, we have the Q -function in the form In the solution of this flash problem, the Legendre transform is used to find out which thermodynamic function should be minimized for the specification, and it is also used in the construction of the Q -function.When working

Accepted Article
This article is protected by copyright.All rights reserved.
with the common flash specifications, one does not need to be aware of the underlying Legendre transform reasoning, since the involved thermodynamic relations are familiar.But for the generalization of the Q -function approach, it is important to understand how the Q -function is constructed through Legendre transforms.In Table 2, we list the Q -functions introduced by Michelsen for six flash specifications together with their Υ and Ξ functions, and the corresponding Legendre transforms.
It is interesting to note that among Michelsen's Q-functions, except for three situations, including TP flash using (T , P )-based thermodynamics, TV flash and TP flash using (T ,V )-based thermodynamics, all the others give a saddle point problem.We attempt to interpret this using the following form of Equation 20: where p k are the conjugate variables for q k and the natural variables for Ξ.We first look at the situation where p consists of extensive variables only.For a two-phase flash with an arbitrary specification, we use r to denote all the independent variables except those fixed by the specifications, and p is a subset of r.If r consists of extensive variables only, the Hessian for Q assumes the following form: The superscripts v, l and z stand for properties related to the vapor phase, the liquid phase and to the total system, respectively.In Equation 21, the independent variables are the ones related to the vapor phase (superscript v) and the ones related to the total system (superscript z).Since the derivatives are given for each phase, we utilize p z = p v + p l and r z = r v + r l to reach the expression in the right-hand side of Equation 21.For a state function obtained from the Legendre transform, it is positive semi-definite (PSD) at the equilibrium in the space or subspace spanned by its extensive variables but negative semi-definite at the equilibrium in the space or subspace spanned by its intensive variables 33 .The PSD for all the extensive variables is required also by the local phase stability.The submatrices A and C are therefore both PSD.If we designed nested loops where the inner one iterates over r v and the outer one over p l in a decoupled manner, both the inner and the outer are minimizations.If p and r consist of both extensive and intensive variables, the submatrices are no longer positive semi-definite at the solution and will be a saddle problem

Accepted Article
This article is protected by copyright.All rights reserved.
by nature.For the TP flash using (V , T )-based thermodynamics, the independent variables (vapor phase moles, vapor volume and liquid volume) are all extensive and the flash can be formulated as minimization of the Q -function.

Open System Flash without Reactions
The or we can write it in its differential form Here, we used the superscript "(k )" to denote that k transforms have been applied to the original state function (A in this case).We can make N C transforms for A and obtain the so-called Landau potential or grand-potential (Φ): The macroscopic system with such a potential corresponds to the grand canonical ensemble in statistical mechanics, where the system is completely open, with chemical potentials specified, and none of the mole numbers are conserved.

Accepted Article
This article is protected by copyright.All rights reserved.
The corresponding flash problem will be a (µ,V , T ) one.The state functions obtained by transforming G are or in its differential form However, if N C transforms of this kind are performed, we get a zero potential as a result: Its differential form is actually the Gibbs-Duhem equation: It indicates that, at constant T and P , we can specify the maximum of N C − 1 chemical potentials or, in other words, we need to specify at least one mole number for the system.
For the flash specifications based on the new state functions, we can construct the corresponding Q -functions again by using Equation 18.For any flash specification in Table 2, it can be extended to a corresponding flash for an open system where k chemical potentials are kept constant.The new function Q -function is: Here, we add the superscript "closed" to denote the Q -functions listed in Table 2.We discuss below the situation where N S components have constant chemical potentials but varying total mole numbers.The updated Q -functions are given in Table 3.Note that the superscript (N S ) denotes that the Legendre transform µ m n m is applied N S times.
For example, G (N S ) (T , P , n , µ ) means that the Υ function is obtained by applying the Legendre transform (−µ m n m ) to G , and G (N S ) needs to be minimized at constant T , P , n † , µ † , where n † and µ † stand for all the mole numbers and chemical potentials, respectively, except for m.These Q -functions can be used in a similar manner to those in Michelsen's work.

Accepted Article
This article is protected by copyright.All rights reserved.
For efficient solution of these problems, Newton's method is needed.The stationary point conditions of a Qfunction are essentially the same as the equal fugacity criteria plus various constraints such as ln f m − ln f spec m = 0.
We limit our discussion to the six (T , P )-based flash specifications in Table 2 and 3, using (T , P )-based thermodynamics, and assume that there are two phases at most.At the stationary point, the derivatives of the Q -function w.r.t. the vapor-phase mole numbers, lnT , ln P and the overall mole number of component m should be equal to zero.
The gradient w.r.t. the same independent variables can have expressions with different scaling factors, depending on whether the Q -function is from G or G /T .They will be rescaled to give the same general form, and further differentiation of the rescaled gradient expressions provides the Jacobian for the Newton method.It turns out that all the six specifications share a common expression of the Jacobian with

Accepted Article
This article is protected by copyright.All rights reserved.
Different flash specifications have different r T , r P and r N , which are listed in Table 4.In Equations 31, we have One should note that but F is not necessarily equal to 1 for the open system.
Here we discuss specifically the (T , P , n † , µ † ) situation.For such a situation, the independent variables for the resulting formulation of Equations 29 and 30 are all extensive, including the vapor phase mole numbers and the total feed moles for those components with specified chemical potentials.

Accepted Article
This article is protected by copyright.All rights reserved.
The problem can be cast as a Q -function minimization.For the two-phase situation with the chemical potential of component 1 specified, we suggest the following procedure to initialize the calculation: at the given T and P , a set of Wilson K -factors are calculated, and the Rachford-Rice equation is solved to yield initial phase fractions (β l and β v ) and initial vapor and liquid phase compositions.Then, the mole fractions of the vapor and liquid phases are used to update the fugacity coefficients using the equation of state.After that, the liquid and vapor mole fractions of component 1 (x 1 and y 1 ) and the total mole number of component 1 (z 1 ) are updated by The new compositions are used to update the K -values.They are calculated as the ratio of fugacity coefficients from an EoS, and then used as input for solving the Rachford-Rice equation in the next round.After five iterations of the above procedure, we switch to the second-order approach using Equation 30.We remark that if the composition derivatives in Equation 30 are dropped, the iteration scheme becomes a first order one.

Open System Flash with Reactions
The new Q -functions defined in section "Open System Flash without Reactions" can be used for systems with both phase equilibrium and reaction equilibrium.The extension to open systems follows the same formula (Equation 29), although reaction equilibrium brings extra complexity to the solution procedure.The methods for solving chemical and phase equilibrium fall into two main categories, stoichiometric and non-stoichiometric approaches 34 .In our recent studies, RAND-based non-stoichiometric methods have been developed and tested.We chose to construct the open system flash solution based on the modified RAND approach discussed by Paterson et al. 15 .and Tsanas et al. 16;18 .
The modified RAND method differs from other CPE methods in that it uses elemental chemical potentials ( λ and λ

Accepted Article
This article is protected by copyright.All rights reserved.
for dimensional and dimensionless version of it, respectively) as independent variables instead of phase mole numbers.
This greatly reduces the size of the system and leads to a more structured matrix for multiphase systems.The phase mole numbers needed for evaluating thermodynamic properties are expressed from the elemental chemical potentials, which involves a slightly more complicated derivation.We follow the derivation of the modified RAND formulation of Paterson et al. 25 and extend it to open systems, where the chemical potentials of N S components are fixed and these N S components are allowed to enter or leave the system.Paterson et al.'s derivation leads to a system of size N E + N P + 2 for six specifications in a closed system, in which N E and N P are the number of elements and of phases, respectively.Our extension adds another N S equations in the system.Since the derivation presented here shares a lot of similarities with the derivation of Paterson et al. 25 , we focus on the most critical steps and emphasize the differences.
For non-stoichiometric methods in a closed system at constant T and P , we choose to work with an augmented objective function, formed by the sum of the Gibbs energy and an element balance constraint: Equation 35 introduces two new variables: the formula matrix (A), and the mole numbers of the elements (b).Elements are abstract entities that are conserved throughout the flash calculation, and they can be understood as the building blocks of the actual components.The system is fully characterized by the mole numbers of the elements in the system (b), and its full specification is (T , P , b).The addition of the second term of the right-hand side of Equation 35to the Gibbs energy does not change the nature of the potential.
A Q -function for a reactive system can be constructed by adding the second term of the right-hand side of Equation 35 to the Q -functions presented in the previous sections.Here, we take the Q -function for the S,V , b † , µ †specification as an example: The stationary conditions with respect to the number of moles in each phase (n i ,j ) and the elemental chemical potential

Accepted Article
This article is protected by copyright.All rights reserved.
( λ) require Equations 37a and 37b also have a clear physical meaning.The first equation represents that the chemical potential of component i in all phases is equal to its chemical potential calculated from the elemental chemical potentials.The second one represents the elemental conservation equation.For an open system, Equation 37a does not change, but Equation 37b needs to be modified and is addressed later in this section.
The chemical potential term in Equation 37a can be linearized in terms of mole numbers, temperature and pressure.After some manipulations 15;25 , we get with e j = 1 RT ∂µ j ∂T , and Multiplying Equation 38 by 1 T and using the relations 1 T ∆n j = ∆β j , 1 T x j = 1 and 1 T M j = x T j , we get: The equation relates the dimensionless equilibrium elemental potentials (λ) to the dimensionless molar Gibbs energy of each phase: Linearization of Equation 37b is less straightforward because the system is not closed anymore and the number of

Accepted Article
This article is protected by copyright.All rights reserved.
moles of elements is no longer constant throughout the whole flash calculation.We express them, as follows: In this equation, n tot is a vector of size N C that collects the total mole numbers of all components in the whole system.
We, then, can rewrite Equation 37b as In the linearization of Equation 43, we need to consider both ∆n i ,j and ∆n tot i , where ∆n i ,j involves all the components in all the phases, and ∆n tot i involves the components whose chemical potentials are specified: − n i ,j + ∆n i ,j + n tot i + ∆n tot i = 0 (44a) Here, A spec is the formula matrix only for the N S components, and it is of the size N E × N S and ∆n spec is the change in the overall mole numbers for the N S components.The introduction of matrix A spec and vector ∆n spec simplifies the notation, as we do not have to carry on with the total mole numbers of species that cannot move through the boundaries of the system.Substitution of Equation 38 into Equation 44b gives If T and P are constant and the system is closed, Equations 41 and 45 provide the N E +N P equations for the multiphase reaction equilibrium.
Regarding the specifications for entropy and pressure, the stationary conditions can be obtained by differentiation of Equation 36 with respect to T and P , respectively.Divided by RT , they are

Accepted Article
This article is protected by copyright.All rights reserved.
Linearization of them with respect to the variables (T , P , n i ,j and n spec m ) yields e T j ∆n j = 0 (47a) For the enthalpy and internal energy specifications, a similar approach can be used, and we would get the following as a result of the linearization: The linearizations presented in Equations 47a, 47b, 48a and 48b are the same as in Paterson et al. ( 2018) 15 , and those equations do not change for open systems.The relation between ξ j and e j is given by Furthermore, it can be shown that Unlike in a closed system, we do not have A N P j =1 ∆n j = 0 throughout the whole process, but the equality is valid at the solution.Therefore, we can also replace ξ j with e j as in Paterson et al. 25 .With other simplifications that Paterson et al. 25 performed, we can rewrite the four equations in two general equations:

Accepted Article
This article is protected by copyright.All rights reserved.
In short, the final forms of these two equations do not change for the open system.The final expressions for these two equations can be obtained by substituting Equation 38 into them and reformulating the equations.
Finally, we differentiate Equation 36with respect to all n spec m .Since vector b is related to n through Equation 42, we have the following stationary conditions: Equation 52 also has a clear physical meaning.It shows that the chemical potential in any phase of the component whose chemical potential is held constant (component m), should be equal the specified chemical potential, through the elemental chemical potentials (λ) of the elements that are the building blocks of component m.The assemblage of Equations 45, 41, 51a, 51b and 52 produces a system of N P + N E + N S + 2 equations: In this system of equations, we have

Accepted Article
This article is protected by copyright.All rights reserved.
β j M j e j (54e) β j e T j M j e j (54k) Here, r T and r P have slightly different expressions than the ones presented in Table 4 because the independent variables in that case were lnT and ln .The expressions for this case are shown in Table 5.
In this work, we have only implemented the (T , P , b † , µ † ) specification for a system modelled by (T , P )-based thermodynamics.We use the initial estimate procedure developed earlier for state function minimization problems proposed by Tsanas et al. 16 .It consists of the minimization of a given function with respect to λ at constant total mole numbers of the system.After a step in λ, the mole fractions of each species in each phase are updated until the stationary conditions of the function are met.

Accepted Article
This article is protected by copyright.All rights reserved.

| Results for the Open System Flash without Reactions
In this section, we present the results of a (T , P )-based implementation of the algorithm outlined in Section 4. A fivecomponent hydrocarbon mixture of methane, ethane, propane, n-octance and n-nonane is used with their critical properties given in Table 6.The SRK EoS is used in all the calculations.
For the reference calculation, two-phase flash at the (T , P , n)-specification is performed using the second-order algorithm suggested by Michelsen 19 .Before switching to the second-order calculation, five successive substitutions are performed without acceleration (Equation 34).The results are summarized in Figure 1. Figure 1 (a) shows the vapor phase fraction in this typical phase envelope for a hydrocarbon mixture.For the open system flash with a specified chemical potential for methane, the specification is (T , P , n † , µ C 1 ).The system has essentially a semi-permeable membrane allowing the methane exchange with a reservoir having a constant fugacity of methane.We actually use a subset of the Equations 30 and 31 because T and P are already specified here.
It should be noted that we only use the fugacity of methane from the (T , P , n)-flash to specify the problem but not its results in assisting the (T , P , n † , µ C 1 )-flash calculation.The calculation algorithm used for the (T , P , n † , µ C 1 )flash is as described in section "Open System Flash without Reactions".The algorithm is a second-order Q -function minimization procedure, but we can reduce it to a first-order one by dropping the composition derivatives of the fugacity coefficients in the Jacobian.The results for the first-order and the second-order calculations are presented in Figures 2 and 3, respectively, for comparison.Both the first-order and the second-order calculations converge and give essentially the same vapor fraction values as in Figure 1.However the first-order one uses much more iterations, especially close to the critical point.
The number of iterations for the second-order (T , P , n † , µ C 1 )-flash is comparable to the corresponding (T , P , n)flash, and the number of iterations is counted after five initial iterations with the initialization method for an open system (see section "Open System Flash without Reactions").The number of iterations does not exceed 18 even

Accepted Article
This article is protected by copyright.All rights reserved.
close to the critical point.Table 7 presents a calculation example at 336.5 K and 6 MPa.We first solved the (T , P , n)flash for z 1 = 0.4.Then we performed calculations at specified fugacity of methane (ln(f 1 /MPa) = 1.538798) starting from z 1 = 0.3 with different algorithms.In addition to the first and second order flash algorithms using Q-function minimization, we tested two other simple algorithms using the PT-flash as their inner loop.The last two simple algorithms are actually one-dimensional root-finding, using the secant method and the Newton method, respectively.
The inner loop at a certain z k 1 gives a residual ln The adjustment is made using either the secant method or the Newton method for z k +1 1 in the next iteration.For the Newton method, its iteration scheme is given by with the derivative Derivative can be obtained from the sensitivity calculation after the convergence of the TP flash in the inner loop (see the supplementary material to this paper for more details).For the secant method, the derivative is approximated by the slope of the secant line.Table 7 shows that all the methods at the specified fugacity converged to the same vapor phase fraction and the total methane amount.The first order Q -function minimization converged after 14 iterations, whereas the other algorithms converged faster.The second order Q -function minimization requires the smallest number (4) of iterations.Although the Newton and the secant root-finding algorithms only needed 5 and 6 outer loop iterations, respectively, each of their outer loop iteration involves several TP flash iterations in the inner loop.The total numbers of iterations calculated from the TP flash iterations in the inner loops are 21 and 24, respectively.It shows that the simple root-finding methods can give quadratic or nearly quadratic convergence in the outer loop, but the total numbers of iterations are much larger than that for the second-order Q -function minimization.We provide more detailed results with different initial z 1 in the supplementary material to this paper, and they all support the same finding.Compared to the root-finding methods, the Q -function minimization handles the problem more systematically.The Q -function can be checked during the iteration to safeguard the convergence.More importantly, the formulation is intended for the general situation with more than one chemical potential being fixed.In this situation,

Accepted Article
This article is protected by copyright.All rights reserved.
the root-finding algorithms will be more complicated.The convergence behavior of the two Q -function methods (first and second-order) and the two root-finding algorithms (the secant method and Newton's method) are shown in Fig- ure 4. The two Q -function methods presented first-order and the second-order convergence, as expected.Newton's method also converged quadratically, but the secant method did not.Figure 4 (c) also shows that both Q -function methods converged to a minimum of the Q -function.
In a second example for the open flash at a constant fugacity of methane, we use the same mixture, but instead of using the methane fugacity at each T and P from Figure 1 (c), we simply specify the methane fugacity to a single fixed value.This test also generates the two-phase region under a constant methane fugacity condition, which is different from our classical TP phase diagram.In the open flash, we fix the moles of all the species except for methane.The system has 0.1, 0.1, 0.2, 0.2 moles of ethane, propane, n-octane and n-nonane, respectively.Methane is free to enter or leave the system, thus making the overall mole fraction of the system vary in the finally presented TP phase diagram for this open system.Figure 5 shows the phase envelope for ln(f 1 /MPa) = 0.6.The pure methane saturation line is also presented in the same figure.The PT phase diagram for the open system is very different from the normal PT phase diagram in Figure 1, with the two-phase region smaller for the open system.The most distinctive feature of the new phase envelope is that it stretches all the way to the pure methane saturation line.Figure 5 (a) shows the vapor phase fraction at different temperatures and pressures.Its value varies from 0 and 1, and there is a critical point around 550 K and 6 MPa to where various vapor phase fraction lines tend to converge.The number of iterations, depicted in Figure 5 (b), shows that the most difficult part to converge is close to the tail of the phase envelope.In this region, the system behavior changes abruptly with small changes in T and P , giving challenges to initialization in this area.Figure 5 (c) shows the total mole fraction of methane in the system.It can be seen that, in the area close to the saturation line, the system consists of nearly pure methane.This means that at those T and P , methane enters the system to such a large extent that the other components at fixed mole numbers become traces.Finally, the dotted and the dash-dotted lines show that the phase envelope meets the saturation line at the saturation temperature and pressure at the specified fugacity of methane.

Accepted Article
This article is protected by copyright.All rights reserved.

| Results for the Open System Flash with Reactions
For the open flash with reactions, we exemplify it with the system containing propene, water, 1-propanol and propane.
Similar systems have been studied previously 16;35;36;37 .There is only one equilibrium reaction in the system: Its equilibrium constant (K ) is assumed to be temperature independent and equal to 23 (K = 23) 37 .The pressure of the system is fixed at 1 bar.The SRK EoS is used in the calculation, and the model parameters are provided in Table 8.
As a reference calculation, the (T , P , b)-flash using the modified RAND approach 16 is performed for the initial composition of 1 mole propene, 1 mole water and 0.1 mole propane.Propane is an inert and its mole number is conserved in all calculations.The results at different temperatures are presented in Figure 6, showing a single-phase liquid region below 270.9 K, two-phase region between 270.9 and 344.6 K, and a single-phase vapor region above 344.6K.The liquid phase below 270.9K is rich in propanol, and the other components are present in small molar fractions (the second most abundant component has its mole fraction lower than 10%).The two-phase region consists of a propanol-rich liquid phase and a vapor phase of varying composition.The vapor phase above 344.6K contains mostly propanol, but the other components are more abundant in this phase than in the liquid below 270.9 K.Besides the variation of phase fractions, liquid mole fractions, and vapor mole fractions, Figure 6 also shows the chemical potentials that can be used in the specifying the open system flash calculations.
In the open system flash, we fix the water chemical potential at the values illustrated in Figure 6 (d).The initial amounts of species are 1 mole propene, 2 moles water and 0.1 mole propane.Water is allowed to leave the system in order to keep the same element amounts as in the (T , P , b) example.Figure 7 presents the results for the open flash with reaction at 1 bar and the same range of T as before.The phase fractions and equilibrium phase compositions are essentially the same as those in Figure 6. Figure 7 demonstrates that the same results from Figure 6 were successfully achieved by using the open system algorithm.
Figure 8 shows the convergence behavior of the (T , P , b † , µ water )-flash using the second-order Q -function minimization at three temperatures (259.9K, 300.5 K, and 350.1K).At 259.9 K and 350.1 K, the system is in single-phase

Accepted Article
This article is protected by copyright.All rights reserved.
regions.There is only a liquid phase and a vapor phase, at 259.9 K and 350.1 K, respectively.At 300.5 K, the system is in the two-phase region.The error in Figure 8 is defined as where q is an iteration number.We remark that this type of error calculation may give artificial increases in the error value.At all the three temperatures, the algorithm converge quadratically.

Conclusions and Future Work
In this work, we analyzed the Q -functions suggested by Michelsen for state function-based flash specifications.We showed that a Q -function relates the thermodynamic state function to be minimized for a flash specification and the state function whose natural variables are used in the flash formulation -Michelsen's paper essentially describes how to perform state function minimization in (T , P , n)-based thermodynamics or (T ,V , n)-based thermodynamics.The two state functions are related by Legendre transforms.The Q -function has the same physical meaning at the solution as the state function to be minimized, and it is not equivalent to the Lagrangian function except at the solution.
Based on the analysis, we extended the Q -function approach to open systems.For an open system, it is possible to define a state function to be minimized with Legendre transforms.We categorized the corresponding equilibrium calculation through minimizing the state function as the flash-type calculation.New Q -functions can be introduced to solve the flash problems for open systems.For open systems without reactions, we suggested Q -functions and a formulation using the classical framework, with the mole numbers as independent variables; for those with reactions, we suggested Q -functions and a formulation using the recently proposed modified RAND framework.For both situations, a common Jacobian can be obtained to cover various specifications.
Flash for open systems at constant T and P leads to Q -function minimization.We discussed specifically how to solve this type of flash for two-phase non-reactive and reactive situations.We validated our algorithms using the Q -function minimization for a non-reactive hydrocarbon mixture and a four-component reactive mixture by a comparison with the classical flash results at the same conditions.For the hydrocarbon mixture where the fugacity of methane is kept constant, a distinctively different two-phase region is obtained for this system allowing the exchange

Accepted Article
This article is protected by copyright.All rights reserved.
of methane with an external reservoir.The algorithms give quadratic convergence for the tested problems.
The general formulations provided in this work include the counterparts for classical PH flash, PS flash, VT flash, and UV flash, although our examples are only for the two-phase open systems at constant T and P .It should also be noted that the formulations are not limited by the number of phases either.In principle, the formulations lay out the basis for robust and efficient algorithms for other types of flash, although it requires more detailed study for the specific problems.Even for isothermal processes only, the current work has made it possible to develop more minimization-based algorithms for many open system problems, such as, compositional grading in a petroleum reservoir, membrane separation processes, and geochemical calculations at constant partial pressure or constant pH.
The last one requires additional considerations on electrolyte models, electroneutrality, and mineral precipitations which will be discussed in detail in our future paper.
TA B L E 1 Thermodynamic potentials and their natural variables for a closed system.

Accepted Article
This article is protected by copyright.All rights reserved.♦ The original variables are (1/T ) and V .

Accepted Article
This article is protected by copyright.All rights reserved.

Accepted Article
This article is protected by copyright.All rights reserved.

Accepted Article
This article is protected by copyright.All rights reserved.F I G U R E 1 Phase envelope calculated using the two-phase second order (T , P , n) flash: (a) vapor phase fraction in the two-phase region; (b) the number of iterations needed in the second-order step; (c) the natural logarithm of the fugacity of methane, which is used to specify the calculation at the fixed methane fugacity.

Accepted Article
This article is protected by copyright.All rights reserved.

Accepted Article
This article is protected by copyright.All rights reserved.

Accepted Article
This article is protected by copyright.All rights reserved.

Accepted Article
This article is protected by copyright.All rights reserved.
transform and the construction of thermodynamic state functions needed in later discussion.Then, we analyze the general theory underlying Michelsen's Q -function construction and the role of the Legendre transform.After that, we show how Michelsen's Q -functions can be extended to open systems.With the extended Q -functions, we present the detailed algorithms for open-system flash without and with reactions.For open systems without reactions, we adopt the classical framework using mole numbers as independent variables; for open systems with reactions, we adopt the RAND framework instead.Finally, we present examples for open systems without and with reactions.

Figure 1 (
b) presents the number of iterations used in the second-order step.Figure1 (c)shows the natural logarithm of the fugacity of methane.Those fugacity values are to be used in the subsequent calculation to specify the chemical potential in an open system.

2
Phase envelope calculated using the two-phase first order (T , P , n † , µ C 1 )-flash: (a) vapor phase fraction in the two-phase region; (b) the number of iterations needed in the first-order step.

3
Phase envelope calculated using the two-phase second order (T , P , n † , µ C 1 )-flash: (a) vapor phase fraction in the two-phase region; (b) the number of iterations needed in the second-order step.
z 1 / z i F I G U R E 5 Phase envelope calculated using the two-phase second order (T , P , n † , µ C 1 )-flash at ln(f 1 /MPa) = 0.6: (a) vapor phase fraction in the two-phase region; (b) the number of iterations needed in the second-order step.(c) total mole fraction of methane in the system (z 1 / z i ).

General Expression for the Q -function | State Function-based Flash Specifications Studied by Michelsen
. Apart from flash for open systems discussed here, Alberty also discussed processes with gravitational and centrifugal work, with surface work, and with electrical work.
classical state function-based flash specifications were developed for closed systems with constant overall mole number.Since Legendre transforms can be used to create other state functions for open systems, the "flash" concept in a general sense, i.e., determination of the component partitioning between phases through minimization, can be extended to open systems.
We construct these new state functions from a basis function, e.g., the Helmholtz energy A or the Gibbs energy G .Each transform forms a new set of natural variables, where the number of moles of a certain component in the old set is replaced by its conjugate potential.The mole numbers that are not replaced are still conserved in the flash calculation.However, those replaced by chemical potentials are no longer conserved, and the system as a whole is no longer a closed one.We start by analyzing the state functions generated from the Helmholtz energy.The new state function from A after k transforms assumes TA B L E 3 Q -functions for open systems.
TA B L E 6 SRK parameters for the components used in the analysis of the open system case without chemical reaction * .
* Binary interaction parameters were all set to 0.
TA B L E 7 Solution for the phase equilibrium problem at 336.5 K and 6 MPa.The number of iterations is counted from the initial value of z 1 = 0.3.B L E 8 SRK parameters for the components used in the analysis of the open system case with chemical reaction * .
†The size of vector n is NC − 1Accepted ArticleThis article is protected by copyright.All rights reserved.TA * Binary interaction parameters were all set to 0.