Construction of a minimum energy path for the VT flash model by an exponential time differencing scheme with the string method
aa r X i v : . [ phy s i c s . c o m p - ph ] A ug Construction of a minimum energy path for the VTflash model by an exponential time differencing schemewith the string method
Yuze Zhang a , Yiteng Li b , Lei Zhang a , Shuyu Sun b, ∗ a Beijing International Center for Mathematical Research, Peking University, Beijing100871, China b Physical Science and Engineering Division (PSE), King Abdullah University ofScience and Technology (KAUST), Thuwal 23955-6900, Saudi Arabia
Abstract
Phase equilibrium calculation, also known as flash calculation, plays signif-icant roles in various aspects of petroleum and chemical industries. SinceMichelsen proposed his milestone studies in 1982, through several decadesof development, the current research interest on flash calculation has beenshifted from accuracy to efficiency, but the ultimate goal remains the samefocusing on estimation of the equilibrium phase amounts and phase composi-tions under the given variable specification. However, finding the transitionroute and its related saddle points are very often helpful to study the evo-lution of phase change and partition. Motivated by this, in this study weapply the string method to find the minimum energy paths and saddle pointsinformation of a single-component VT flash model with the Peng-Robinsonequation of state. As the system has strong stiffness, common ordinary dif-ferential equation solvers have their limitations. To overcome these issues, aRosenbrock-type exponential time differencing scheme is employed to reducethe computational difficulty caused by the high stiffness of the investigatedsystem. In comparison with the published results and experimental data,1he proposed numerical algorithm not only shows good feasibility and accu-racy on phase equilibrium calculation, but also successfully calculates theminimum energy path and and saddle point of the single-component VTflash model with strong stiffness.
Keywords:
VT flash, Peng-Robinson equation of state, minimum energypath, string method, Rosenbrock-type exponential time differencingscheme
1. Introduction
The accurate knowledge of phase equilibria is of vital importance inpetroleum industry to determine the number of equilibrium phases andtheir amounts and compositions for complex reservoir fluids, which stimu-lates the development of equation-of-state-based phase equilibrium calcu-lation. As one of important applications, phase equilibrium calculation isan integral component of compositional multiphase flow simulators that areused to model recovery processes sensitive to compositional changes, suchas primary depletion of volatile oil and gas condensate reservoirs and mul-tiple contact miscible flooding [1]. Moreover, it can be used as standalonecalculation as well for analyzing, designing and optimizing processes andfacilities in order to reduce undesired species (e.g. H O) for conforming toproduct specifications, remove acid gas (e.g. CO and H S) for protectingpipelines and equipment from corrosion [2–4], determine the amount of in- ∗ Corresponding author: Shuyu Sun.
Email addresses: (Yuze Zhang), [email protected] (Yiteng Li), [email protected] (Lei Zhang), [email protected] (Shuyu Sun)
Preprint submitted to Fluid Phase Equilibria August 11, 2020 ibitor (e.g. methanol and glycols) for avoiding gas hydrate formation [5–8],control asphaltene precipitation for enhancing flow assurance [9–11], etc.The most commonly-used phase equilibrium calculation is conductedat the specified pressure ( P ), temperature ( T ) and chemical composition( N ), which is known as PT flash calculation in the petroleum industry.Usually, the Wilson correlation is used to initialize flash calculation but thisapproximation may not converge to the equilibrium solution. To overcomethis issue, stability test can be used to determine whether phase split takesplace and it will give a better initial approximation if the investigated fluidis unstable at the given condition but requires additional computationalcost. The milestone studies of Michelsen in stability testing [12] and phasesplitting calculation [13] laid a foundation for the development of PT flashcalculation using a single consistent equation of state (EOS), such as theSoave-Redlich-Kwong (SRK) EOS [14] and the Peng-Robinson (PR) EOS[15].Recently, the computation of phase equilibria under the given volume( V ), temperature ( T ) and mole composition ( N ), also known as VT flashcalculation, has become a promising alternative of the conventional PTflash. In addition to its well-posed formulation, VT flash exhibits an un-ambiguous relation between pressure and volume so that the root selectionprocedure in PT flash framework can be avoided [16]. Without the pre-processing of successive substitution iteration, the second-order VT flashmethods exhibit higher computational efficiency than the second-order PTflash methods in stability testing. Also, it can improve the overall efficiencyof phase split calculation while the performance somewhat depends on theinvestigated systems and thermodynamic models [17]. Another important3dvantage of VT flash formulation is phase behavior modeling for unconven-tional hydrocarbon mixtures in the presence of capillary pressure [18–21],since the commonly-used interfacial tension models are explicit functionof molar density rather than pressure. Numerous efforts have been madein the past few years to expand the applications of VT flash not only invarious phase equilibrium problems [22–24] but also in compositional flowsimulation [25, 26].Clearly, all the aforementioned phase equilibrium calculation focuses onthe estimation of phase amounts and phase compositions at the equilib-rium state where either the Gibbs free energy or Helmholtz free energy isminimized at the end of PT or VT flash calculation respectively. How-ever, finding transition routes between different local minima, one of whichactually represents the equilibrium state, under certain conditions is oftenhelpful to study the evolution of phase partitioning. These transition routescan be called as minimum energy paths (MEPs), where the potential forceis always parallel to the path. In addition, we can get the first-order saddlepoints’ information of the original energy landscape of the investigated sys-tem, which are bottlenecks for a large number of particular barrier-crossingevents. Furthermore, the knowledge of high-order saddle points is helpfulin finding multiple local minima [28–30], as well as the relative minimumpoint of energy stability, and even the global minimum point can be obtainedthrough enumeration. These knowledge could help phase behavior modelingof complex reservoir fluids. In this study, we focus on calculating the MEPand first-order saddle points of a single-component VT flash model. Thereare several numerical schemes proposed to calculate the MEPs [27, 31–35].Among them, the string method [27, 34] exhibits a good performance and4nables to efficiently determine energy transition pathways for complex sys-tems with smooth energy landscape by evolving strings rather than pointsin configuration space. This approach consists of two main steps, evolu-tion of the string by some ordinary differential equation (ODE) solvers andreparametrization of the string by interpolation. The accuracy of this ap-proach significantly relies on ODE solvers. It is worth mentioning thatequation-of-state-based VT flash models could have ill-conditioned Hessianmatrices with very large condition numbers (For example, in this study PREOS is employed to model a single-component system, the two eigenvaluesof the Hessian matrix differ by 10 − times). Moreover, explicit numer-ical approaches are required because the real scale of the problem is huge.Some common ODE solvers like Euler and Runge-Kutta methods cannotconverge or require a very small time step. For these problems with strongstiffness, an exponential time differencing (ETD) scheme is used here toreduce the influence of the rigid term during numerical calculation.The ETD scheme involves exact integration of the governing equationsfollowed by an explicit approximation of a temporal integral for the nonlin-ear terms. It can provide satisfactory accuracy and stability even for prob-lems with strong stiffness. This scheme was originally applied in the fieldof computational electrodynamics [36]. Then it was systematically studiedin [37] and was developed for solving stiff systems by Cox and Matthews[38]. Also, in [38], high-order multi-step ETD schemes and Rounge-Kuttaversions of ETD schemes were discussed. The linear stabilities of some ETDand modified ETD schemes were investigated by Du and Zhu [39, 40]. Forsome problems with very strong stiff terms, a bad chosen linearization ofordinary ETD schemes may cause a severe step size restriction. To avoid5hese problems, Hochbruck and Ostermann proposed Rosenbrock-type ETDscheme [41, 42] to handle the stiff nonlinear term without loss of the originalstability. Thus, in this study, we use the Rosenbrock-Euler ETD scheme tosolve the single-component VT flash model with strong stiffness. To the bestof the authors’ knowledge, this is the first time that the energy landscapeof the VT flash problem is modeled by ETD schemes. The structure of thispaper is organized as follows. In Section 2, we describe the VT flash modelusing PR EOS and give details for the single-component system. In Section3, the framework of string method is provided. In Section 4, we establishthe first-order and second-order ETD schemes for the single-component VTflash problem. In Section 5, numerical examples are presented to demon-strate the performance of the proposed ETD schemes and construct theMEP for the investigated system. At the end, we make some remarks andconclusions.
2. VT flash calculation
Let us consider a reservoir fluid mixture of M components occupiesvolume V at temperature T . If the investigated mixture is unstable andsplits into at most two phases, the total Helmholtz free energy of a two-phase system is defined by F = F ( N , V , T ) + F ( N , V , T ) , (1)where N α = [ N ,α , . . . , N M,α ] T and V α denote the mole number vector andvolume of phase α = 1 ,
2, respectively. For brevity, we will drop T in thefollowing expressions of F since T is assumed to be constant. Given that thetotal mole number of each component ( N i ) and total volume ( V ) are fixed,6e can arbitrarily select the mole composition and volume of one phase asprimary variables and compute the counterparts of the other phase by themole and volume balance equations N i, + N i, = N i , i = 1 , . . . , M (2) V + V = V . (3)The second law of thermodynamics indicates the minimum principle ofHelmholtz free energy at given N , V and T so that the phase equilibriumproblems can be formulated as the following minimization problem [43] if N and V are chosen as the primary variablesmin F ( N , V ) (4)subject to 0 ≤ N i, ≤ N i , (5a) V − M X i =1 b i N i, > , (5b)( V − V ) − M X i =1 b i ( N i − N i, ) > , (5c)where b i is the co-volume parameter of pure component. For details, theparameters of PR EOS can be found in Appendix A. In terms of the relation F = − P V + P i µ i N i , the partial derivatives of F ( N , V ) with respect to N i, and V yield ∂F∂N i, = ∂F ( N , V ) ∂N i, + ∂F ( N , V ) ∂N i, = µ i, − µ i, , (6) ∂F∂V = ∂F ( N , V ) ∂V + ∂F ( N , V ) ∂V = P − P , (7)7here µ i,α and P α are the checmical potential of component i and pressurein phase α . Equation (6) and (7) are known as the chemical equilibriumcondition and mechanical equilibrium condition, respectively. At the equi-librium state where F is minimized, µ i, = µ i, holds for all the componentsand meanwhile P = P is achieved.To efficiently solve bulk phase equilibria problems, one popular approachis making use of the local curvature of the energy function. The symmetricHessian matrix has been extensively used to design fast and robust numer-ical algorithms for VT flash calculations [16, 44–47]. With a proper linesearch scheme, the original constrained minimization problem, comprising(4) and (5a)–(5c), can be reformulated as an unconstrained minimizationproblem and it can be solved as follows H ∆ x = − g , (8)where H = B CC T D , ∆ x = ∆ N i, ...∆ N M, ∆ V , g = µ i, − µ i, ... µ M, − µ M, − P + P , (9)with B ij = ∂µ i, ∂N j, + ∂µ i, ∂N j, , (10a) C j = ∂µ i, ∂V + ∂µ i, ∂V = − ∂P ∂N i, − ∂P ∂N i, , (10b) D = − ∂P ∂V − ∂P ∂V . (10c)In order to ensure a constant decline of Helmholtz free energy over itera-tions, H is required to be positive definite as indicated in the aforementioned8iterature. Consequently, the positive definiteness of H is examined at eachiteration for the sake of safety, and then a modified Cholesky factorization[48, 49] is performed if H is found not sufficiently positive definite. Clearly,such a process requires additional operations on the Hessian matrix. Onthe other hand, it has been reported in [26] that the modified Choleskyalgorithm may hinder the convergence sometimes, although it is truly help-ful in most instances. In comparison, the proposed numerical algorithm,which combines the ETD scheme and the string method, circumvents themodified Cholesky algorithm. For a single-component fluid system, H is a2 × i that is commonly used as component index. Then the chemicalpotential and pressure of pure component reads as µ = − RT ln (cid:18) V − BN (cid:19) + RT BV − B + A √ BN ln V + (1 − √ BV + (1 + √ B ! − aNVV + 2 BV − B , (11) P = N RTV − B − AV + 2 BV − B , (12)and their derivatives with respect to mole number and volume are given asbelow (cid:18) ∂µ∂N (cid:19) V,T = RT V N ( V − B ) − aVV + 2 BV − B − aV (cid:16) V + B (cid:17)(cid:16) V + 2 BV − B (cid:17) , (13) (cid:18) ∂µ∂V (cid:19) N,T = − (cid:18) ∂P∂N (cid:19) V = − RT V ( V − B ) + 2 aN V (cid:16) V + B (cid:17)(cid:16) V + 2 BV − B (cid:17) , (14) (cid:18) ∂P∂V (cid:19) N,T = − N RT ( V − B ) + 2 A (cid:16) V + B (cid:17)(cid:16) V + 2 BV − B (cid:17) . (15)9n the above expressions, A = aN , B = bN and the universal gas constant R ≈ . · m ) / (K · mol). The energy parameter a is defined inAppendix A as well. In this study, we will calculate the minimum energypath and saddle point for a single-component VT flash model by solving ∂n g ∂t = − µ ( n g ) V g + µ ( n l ) V g ,∂V g ∂t = µ ( n g ) n g − f ( n g ) − (cid:16) µ ( n l ) n l − f ( n l ) (cid:17) , (16)where n g = N g V g , n l = N l V l and f = FV represents the Helmholtz free energydensity.
3. The string method
The main purpose of using the string method is to find the minimumenergy path (MEP) for barrier-crossing events. If the potential energy F of a system has at least two local minima m and m , a MEP is a curve γ which connects these two local minima and tangents to ∇ F ∇ F ( γ ) − ( ∇ F ( γ ) , ˆ γ )ˆ γ = 0 , (17)where ˆ γ denotes the unit tangent of γ . Basically, finding the MEP bystring is to evolve a curve connecting two local minima under some certaindynamics and the overall algorithm is an iterative application of two mainprocedures: using ordinary differential equation (ODE) solvers to evaluatethe string and reparameterizing the string by interpolation.For the first step, as shown in [27], N discrete points along the stringare evolved by ( φ i ) t = −∇ F ( φ i ) , i = 1 , . . . , N . (18)10or the reparametrization procedure, it is worth pointing out that onlythe normal component of the curve’s velocity affects the curve, while thetangential velocity only contributes to the movement of points along thecurve. Thus, it is free to choose any parametrization and here a particularparametrization of the curve Γ = { φ ( α ) : α ∈ [0 , } is used. Wheneverthe curve evolves, points on the Γ will change their positions following thedynamic (18). It is then necessary to interpolate these points and theirvalues onto a same parametrization Γ. Here we give two common ways ofparametrization. Parametrization by equal arc length
In this parametrization, we have two steps:(1) Calculate the arc length l i based on the current point set { φ i } of thecurve l i = l i − + | φ i − φ i − | , i = 1 , . . . , N , where l = 0 and then the nonuniform mesh { m i } can be obtained bynormalizing { l i } m i = l i /l N . (2) Estimate new φ ∗ i at the uniform mesh { m ∗ i = i/N } under given { m i } and { φ i } by interpolation. Parametrization by energy-weight arc length
It is evident that the saddle point in MEP has the largest energy. In orderto give finer resolution around the saddle point, another parametrizationprocess is carried out, which includes the following two steps:(1) Calculate the energy-weighted arc length l i based on the current point11et { φ i } of the curve l i = l i − + W ( F ( φ i − ) + F ( φ i )2 ) | φ i − φ i − | , i = 1 , . . . , N , where l = 0 and W ( x ) is weight function which can be some well-selectedpositive and increasing functions for x ∈ R . Then the nonuniform meshcan be generated under the same procedure as in the first parametrizationprocess.(2) Estimate new φ ∗ i at the uniform mesh { m ∗ i = i/N } under given { m i } and { φ i } .
4. Rosenbrock exponential time differencing schemes for stiff sys-tem with PR EOS
For VT flash problems with strong stiffness (16), conventional ODEsolvers do not converge and have large errors. Instead, we use Rosenbrock-Euler exponential time differencing schemes to reduce the numerical diffi-culty caused by the stiffness of the Hessian matrix. For this purpose, firsta linear term is introduced into the original system and the new system iswritten as dudt = H ( u n ) u + G ( u, u n ) , (19)where H ( u n ) = H n is the Hessian matrix of F ( u n ) , u n is the molar densityat n th time step and G ( u, u n ) = −∇ F ( u ) − H ( u n ) u . From (19), it gives ( e − tH n u ) t = e − tH n G ( u, u n ) . (20)12hen we integrate (20) on both sides from 0 to one time step τ , which yields e t n +1 H n u n +1 − e − t n H n u n = Z τ e ( t n + s ) H n G ( u ( t n + s ) , u n ) ds ,u n +1 = e τH n u n + Z τ e ( τ − s ) H n G ( u ( t n + s ) , u n ) ds . (21)We can see that, different from ordinary ODE solvers, the Rosenbrock ETDscheme gives the exact relationship between u n +1 and u n . To some extent,the influence of Hessian matrix on the convergence performance is negligibleand thus it allows us to avoid using the modified Cholesky algorithm thatis extensively used in the literature. For the term u ( t n + s ) in (21), if we usedifferent numerical methods to approximate it, we can get different ETDschemes for the system (16). For instance, by replacing u ( t n + s ) with u n ,the following first-order ETD scheme (ETD1) is established u n +1 = e τH n u n + H − n ( e τH n − I ) G ( u n , u n ) , (22)where I is the identity matrix. However, in some cases, the high-accuracysolution is needed to reduce the difficulty of numerical calculation causedby stiff system. Here we introduce a second-order ETD (ETD2) scheme inwhich u ( t n + s ) is approximated by linear interpolation˜ u n +1 = e τH n u n + H − n ( e τH n − I ) G ( u n H n , u n ) ,u n +1 = e τH n u n + Z τ e ( τ − s ) H n [(1 − sτ ) G ( u n , u n ) + sτ G ( ˜ u n +1 , u n )] ds , = e τH n u n + τ [( φ ( τ H n ) − φ ( τ H n )) G ( u n , u n ) + φ ( τ H n ) H ( u n +1 , u n )] , where φ ( α ) = ( e α − I ) α − ,φ ( α ) = ( e α − I − α )( α ) − . able 1: Compositional parameters of the investigated species in numerical examples. Component T c ,i (K) P c ,i (MPa) M w ,i (g / mol) ω i nC
5. Numerical experiment
In this section, we calculate the MEP and the saddle-point for the single-component VT flash model. We first use Rosenbrock-Euler-ETD schemeto solve the single-component VT flash problem. The numerical resultswill be compared with both published results and experimental data todemonstrate the feasibility of the Rosenbrock-Euler-ETD scheme in dealingwith the VT flash problem with strong stiffness.
The performance of the Rosenbrock-Euler-ETD scheme on the VT flashmodel using PR EOS is tested for both hydrocarbon component and non-hydrocarbon component. Compositional parameters of the investigatedspecies are shown in Table 1. To confirm the accuracy and credibility ofthe proposed algorithm, the computed results are compared with the ex-perimental data, which are obtained from NIST chemistry webbook [50].Numerical experiments are conducted in two ways: A) fixing the total molenumber N and changing the temperature T ; B) fixing the temperature T and changing the total mole number N . The total volume V in both casesare same and assumed to be unity.In the numerical experiment of hydrocarbon component, the investi-gated species is n-butane (nC ). For the experiment A, the total particle14
00 320 340 360 380 400
Temperature (K) M o l a r den s i y ( m o l / m ) ETDSemi-implicit scheme [19]Laboratory data [50]
Figure 1: Equilibrium molar density of the gas phase at T = 300 , , , ,
400 Kand the comparison with published results and experimental data (nC ) number is N = 3000 mol. Temperatures in the numerical experiment are300, 325, 350, 375 and 400 K. The time step is 10 − . For the experi-ment B, the temperature is fixed to 350 K and the total mole number N = 2000 , , , , . For the experiment A, the total particle num-ber is N = 10000 mol. Temperatures in the numerical experiment are220, 240, 260, 280 and 300 K. The time step is 10 − . For the experi-ment B, the temperature is fixed to 280 K and the total mole number N = 5000 , , , , and CO at the equilibrium state, respectively. As temperature increases, the molardensity of gas phase goes up because an increasing number of moleculesevaporate into the gas phase. To test the prediction accuracy of the pro-15
20 230 240 250 260 270 280 290 300
Temperature (K) M o l a r den s i t y ( m o l / m ) ETDSemi-implicit scheme [19]Laboratory data [50]
Figure 2: Equilibrium molar density of the gas phase at T = 220 , , , ,
300 Kand the comparison with published results and experimental data (CO ) posed ETD scheme, we compare the computed results with our previouswork [19] and experimental data. The inset figures show more details whichcannot be clearly seen from the parent figure. It can be seen that the molardensities calculated by the Rosenbrock-Euler-ETD scheme agree with theexperimental data very well. By comparing with the results of our pre-vious work, the explicit scheme we used in this paper exhibits the sameaccuracy as the semi-implicit scheme in [19]. Table 2 and Table 3 give theequilibrium mole number, volume and molar density of gaseous nC andCO at fixed temperatures, respectively. The constant molar densities canbe explained in terms of the Gibbs phase rule. Since we fix the tempera-ture, for a single-component system, the pressure is fixed as well. At thesame isothermal-isobaric condition, the densities of gas and liquid phasekeep constant so that the variation of gas volume is identically proportionalto the variation of the mole number of gas phase. As we can see, our new16cheme works better in the constant temperature cases than the previouswork [19]. Table 2: Equilibrium mole number ( N g ), volume ( V g ) and molar density ( n g ) of gaseousnC at T = 350 K with total mole number N = 2000 , , , , N N g (ETD) V g (ETD) n g (ETD) n g ([19]) laboratory data2000 326.3 0.812 402.2 402.1 402.43500 254.9 0.634 402.2 402.1 402.45000 184.2 0.458 402.2 402.1 402.46500 113.0 0.281 402.2 402.1 402.48000 41.8 0.104 402.2 402.1 402.4 Table 3: Equilibrium mole number ( N g ), volume ( V g ) and molar density ( n g ) of gaseousCO at T = 280 K with total mole number N = 5000 , , , , N N g (ETD) V g (ETD) n g (ETD) n g ([19]) laboratory data5000 2388.4 0.865 2763.3 2758.7 2766.37500 1973.7 0.715 2763.3 2758.7 2766.310000 1559.1 0.565 2763.3 2758.7 2766.312500 1144.4 0.415 2763.3 2758.7 2766.315000 729.7 0.2644 2763.3 2758.7 2766.3 In this section, an example will be presented to calculate the MEP andone of the saddle point of the VT flash model for pure nC at temperature T = 350 K. The total particle number is 3000 mol and the time step is 10 − .During the calculation process, the MEP was approximated using the string17 igure 3: Initial configuration of string in the energy landscape for nC method with 100 images and the initial string is the linear interpolationbetween (316 , .
69) and (2700 , . t | u n +1 i − u ni | < TOL , where the tolerance TOL = 10 − .The calculated MEP is shown in Figure 4. If we use the two-step proce-dure introduced in [27], We can determine the location of the saddle point(1500 , .
5) of the energy landscape by using the following climbing imagealgorithm φ ′ t = −∇ F ( φ ′ ) + 2( ∇ F ( φ ′ ) , ˆ γ ′ ) ˆ γ ′ , (23)18 igure 4: Minimum energy path of nC where φ ′ is the saddle point of ∇ F = 0 and γ ′ is the related unit tangentvector of the MEP. It can be proved that the stable solution of (23) is theunstable solution of ∇ F = 0. The initial value of the algorithm requiresthat the initial φ ′ is close to the actual saddle point. Here, we use the stringmethod to calculate a rough approximation of MEP with small number ofimages (in this work, we choose 25 images). Then φ ′ can be selected as thepoint with the highest potential energy on the MEP and γ ′ is the relatedunit tangent vector.
6. Conclusion
In this study, we have applied the string method to calculate the mini-mum energy path of the single-component VT flash model with the Peng-Robinson equation of state. Especially, in order to avoid the computationaldifficulties brought by strong stiffness of the model, the Rosenbrock-typeETD scheme, which is believed superior to other common ODE solvers, is19sed to evolve images of the string in each time step. Numerical experimentsare carried out and show good feasibility and accuracy of the proposed algo-rithm by comparing the computed results with the results of previous workand experimental data. It is worth pointing out that the local minima hasbeen connected by a MEP of the given model and consequently its saddlepoint can be accurately calculated. This work has discussed the relationbetween the local minima of a single-component VT flash system and itssaddle point for the first time. However, it is still a difficult task either toget the high-order saddle points or to find the global minimum energy pointfor a multi-component VT flash system, which is the subject of our futurework.
Acknowledgments
The authors greatly thank for the support from the National NaturalScience Foundation of China (grant number 51874262, 51936001,11802090)and the Research Funding from King Abdullah University of Science andTechnology (KAUST) through the grants BAS/1/1351-01, URF/1/4074-01,URF/1/3769-01, and REP/1/2879-01.
Appendix A. Parameters of the Peng-Robinson equation of state
The Peng-Robinson equation of state has the following form P ( N , V, T ) = N RTV − bN − aN V + 2 bN V − ( bN ) , (A.1)and accordingly the Helmholtz free energy can be defined as F ( N , V, T ) = RT M X i =1 N i (cid:18) ln N i V − (cid:19) − NRT ln (cid:18) V − BV (cid:19) + A √ B ln V + (1 − √ BV + (1 + √ B ! , (A.2)20here the coefficients A = aN and B = bN , and N = P i N i is the totalmole number. The energy parameter a and co-volume parameter b of afluid mixture are computed based on the classical Van der Waals mixingrule such that a = X i X j x i x j ( a i a j ) / (1 − k ij ) , b = X i x i b i . (A.3)In the above expression, k ij is the binary interaction coefficient. For pure-component fluid systems, x i = 1 and k ij = 0 so that a and b are reduced to a = a i = 0 . R T c P c m − r TT c !! , (A.4) b = b i = 0 . RT c P c . (A.5)Here, T c and P c represent the critical temperature and the critical pressureof a pure substance, respectively. The parameter m is a function of theacentric factor ω of pure substance m = 0 . . ω − . ω , ω ≤ . , (A.6) m = 0 . . ω − . ω + 0 . ω , ω > . . (A.7)If ω is unavailable, it can be estimated by the following correlation ω = 37 log (cid:18) P c (cid:19) T c T b − − , (A.8)where T b is the normal boiling point.21 eferences [1] Chen, Z. (2007). Reservoir simulation: mathematical techniques in oil recovery.SIAM.[2] Poormohammadian, S.J., Lashanizadegan, A., Salooki, M.K. (2015). Modelling VLEdata of CO and H S in aqueous solutions of N-methyldiethanolamine based on non-random mixing rules. International Journal of Greenhouse Gas Control, 42, 87–97.[3] Wang, T., El Ahmar, E., Coquelet, C. (2017). Alkane solubilities in aqueous alka-nolamine solutions with CPA EoS. Fluid Phase Equilibria, 434, 93–101.[4] Wang, T., El Ahmar, E., Coquelet, C., Kontogeorgis, G. M. (2018). Improvementof the PR-CPA equation of state for modelling of acid gases solubilities in aqueousalkanolamine solutions. Fluid Phase Equilibria, 471, 74–87.[5] Anderson, F.E., Prausnitz, J.M. (1986). Inhibition of gas hydrates by methanol.AIChE journal, 32(8), 1321–1333.[6] Pedersen, K.S., Michelsen, M.L., Fredheim, A.O. (1996). Phase equilibrium cal-culations for unprocessed well streams containing hydrate inhibitors. Fluid PhaseEquilibria, 126(1), 13–28.[7] Haghighi, H., Chapoy, A., Burgess, R., Tohidi, B. (2009). Experimental and ther-modynamic modelling of systems containing water and ethylene glycol: Applicationto flow assurance and gas processing. Fluid Phase Equilibria, 276(1), 24–30.[8] Burgass, R., Chapoy, A., Li, X. (2018). Gas hydrate equilibria in the presence ofmonoethylene glycol, sodium chloride and sodium bromide at pressures up to 150MPa. The Journal of Chemical Thermodynamics, 118, 193–197.[9] Sabbagh, O., Akbarzadeh, K., Badamchi-Zadeh, A., Svrcek, W.Y., Yarranton, H.W.(2006). Applying the PR-EoS to asphaltene precipitation from n-alkane dilutedheavy oils and bitumens. Energy & fuels, 20(2), 625–634.[10] Shirani, B., Nikazar, M., Mousavi-Dehghani, S.A. (2012). Prediction of asphaltenephase behavior in live oil with CPA equation of state. Fuel, 97, 89–96.[11] Jindrov´a, T., Mikyˇska, J., Firoozabadi, A. (2016). Phase Behavior Modeling ofBitumen and Light Normal Alkanes and CO2 by PR-EOS and CPA-EOS. Energy& Fuels, 30(1), 515–525.
12] Michelsen, M.L. (1982). The isothermal flash problem. Part I. Stability. Fluid phaseequilibria, 9(1), 1–19.[13] Michelsen, M.L. (1982). The isothermal flash problem. Part II. Phase-split calcula-tion. Fluid phase equilibria, 9(1), 21–40.[14] Soave, G. (1972). Equilibrium constants from a modified Redlich-Kwong equationof state. Chemical engineering science, 27(6), 1197–1203.[15] Peng, D., Robinson, D.B. (1976). A new two-constant equation of state. Industrialand Engineering Chemistry Fundamentals, 15(1), 59–64.[16] Jindrov´a, T., Mikyˇska, J. (2013). Fast and robust algorithm for calculation of two-phase equilibria at given volume, temperature, and moles. Fluid Phase Equilibria,353, 101–114.[17] Liang, X. (2018). Numerical aspects of phase equilibrium calculations with the cu-bic and association models. Industrial & Engineering Chemistry Research, 57(42),14273–14285.[18] Sandoval, D.R., Michelsen, M.L., Yan, W., Stenby, E.H. (2019). VT-based phaseenvelope and flash calculations in the presence of capillary pressure. Industrial &Engineering Chemistry Research, 58(13), 5291–5300.[19] Li, Y., Kou, J., Sun, S. (2018). Thermodynamically stable two-phase equilibriumcalculation of hydrocarbon mixtures with capillary pressure. Industrial & Engineer-ing Chemistry Research, 57(50), 17276–17288.[20] Sun, S. (2019). Darcy-scale phase equilibrium modeling with gravity and capillarity.Journal of Computational Physics, 399, 108908.[21] Nichita, D.V. (2019). Density-based phase envelope construction including capillarypressure. Fluid Phase Equilibria, 498, 33–44.[22] Travalloni, L., Castier, M., Tavares, F.W. (2014). Phase equilibrium of fluids con-fined in porous media from an extended Peng-Robinson equation of state. FluidPhase Equilibria, 362, 335–341.[23] Luo, S., Lutkenhaus, J.L., Nasrabadi, H. (2018). Multiscale fluid-phase-behaviorsimulation in shale reservoirs using a pore-size-dependent equation of state. SPEReservoir Evaluation & Engineering, 21(04), 806–820.[24] Li, Y., Qiao, Z., Sun, S., Zhang, T. (2020). Thermodynamic modeling of CO olubility in saline water using NVT flash with the cubic-plus-association equationof state. Fluid Phase Equilibria, 112657.[25] Pol´ıvka, O., Mikyˇska, J. (2014). Compositional modeling in porous media using con-stant volume flash and flux computation without the need for phase identification.Journal of Computational Physics, 272, 149–169.[26] Paterson, D., Michelsen, M.L., Stenby, E.H., Yan, W. (2018). Compositional Reser-voir Simulation Using (T, V) Variables-Based Flash Calculation. In ECMOR XVI-16th European Conference on the Mathematics of Oil Recovery (Vol. 2018, No. 1,pp. 1-15). European Association of Geoscientists & Engineers.[27] Weinan E, Ren W, Vanden-Eijnden E. (2008). Simplified and improved stringmethod for computing the minimum energy paths in barrier-crossing events. Journalof Chemical Physics, 126(16), 164103.[28] Yin, J., Yu, B., Zhang, L. (2020). Searching the solution landscape by generalizedhigh-index saddle dynamics, arXiv preprint arXiv:2002.10690.[29] Yin J., Zhang L., Zhang P. (2019). High-Index Optimization-Based Shrinking DimerMethod for Finding High-Index Saddle Points. SIAM Journal on Scientific Comput-ing, 41(6), A3576–A3595.[30] Yin, J., Wang, Y., Chen, J.Z., Zhang, P., Zhang, L. (2020). Construction of apathway map on a complicated energy landscape. Physical Review Letters, 124(9),090601.[31] Ulitsky, A., Elber, R. (1990). A new technique to calculate steepest descent pathsin flexible polyatomic systems. The Journal of chemical physics, 92(2), 1510–1511.[32] Fischer, S., Karplus, M. (1992). Conjugate peak refinement: an algorithm for find-ing reaction paths and accurate transition states in systems with many degrees offreedom. Chemical physics letters, 194(3), 252–261.[33] Henkelman, G., Uberuaga, B. P., Jnsson H. (2000). A climbing image nudged elasticband method for finding saddle points and minimum energy paths. The Journal ofChemical Physics, 113(22), 9901–9904.[34] Weinan E., Ren W., Vanden-Eijnden E. (2002). String method for the study of rareevents. Physical Review B, 66(5), 052301.[35] Berne B.J., Ciccotti G., Coker D.F. (1998). Classical and quantum dynamics in ondensed phase simulations: Proceedings of the International School of Physics.World Scientific, Singapore.[36] Taflove, A., Hagness, S.C. (2005). Computational electrodynamics: the finite-difference time-domain method. Artech house, London.[37] Beylkin, G., Keiser, J.M., Vozovoi, L. (1998). A new class of time discretizationschemes for the solution of nonlinear PDEs. Journal of Computational Physics,147(2), 362–387.[38] Cox, S.M., Matthews, P.C. (2002). Exponential time differencing for stiff systems.Journal of Computational Physics, 176(2), 430–455.[39] Du, Q., Zhu, W. X. (2004). Stability analysis and application of the exponentialtime differencing schemes. Journal of Computational Mathematics, 22, 200–209.[40] Du, Q., Zhu, W. (2005). Analysis and applications of the exponential time differ-encing schemes and their contour integration modifications. BIT Numerical Math-ematics, 45(2), 307–328.[41] Hochbruck, M., Ostermann, A., Schweitzer, J. (2006). Exponential integrators ofRosenbrock-type. Oberwolfach Reports, 3, 1107–1110.[42] Hochbruck, M., Ostermann, A., Schweitzer, J. (2009). Exponential Rosenbrock-typemethods. SIAM Journal on Numerical Analysis, 47(1), 786–803.[43] Firoozabadi, A. (2015). Thermodynamics and applications of hydrocarbon energyproduction. McGraw Hill Professional.[44] Mikyˇska, J., Firoozabadi, A. (2012). Investigation of mixture stability at given vol-ume, temperature, and number of moles. Fluid Phase Equilibria, 321, 1–9.[45] Jindrov´a, T., Mikyˇska, J. (2015). General algorithm for multiphase equilibria calcu-lation at given volume, temperature, and moles. Fluid Phase Equilibria, 393, 7–25.[46] Nichita, D.V. (2017). Fast and robust phase stability testing at isothermal-isochoricconditions. Fluid Phase Equilibria, 447, 107–124.[47] Nichita, D.V. (2018). New unconstrained minimization methods for robust flash cal-culations at temperature, volume and moles specifications. Fluid Phase Equilibria,466, 31–47.[48] Gill, P.E., Murray, W. (1974). Newton-type methods for unconstrained and linearlyconstrained optimization. Mathematical Programming, 7(1), 311–350.
49] Schnabel, R.B., Eskow, E. (1999). A revised modified Cholesky factorization algo-rithm. SIAM Journal on optimization, 9(4), 1135–1148.[50] Linstrom, P.J., Mallard, W.G. NIST Chemistry WebBook, NIST Standard Ref-erence Database Number 69, National Institute of Standards and Technology,Gaithersburg MD, 20899, https://doi.org/10.18434/T4D303, (retrieved June 12,2020)49] Schnabel, R.B., Eskow, E. (1999). A revised modified Cholesky factorization algo-rithm. SIAM Journal on optimization, 9(4), 1135–1148.[50] Linstrom, P.J., Mallard, W.G. NIST Chemistry WebBook, NIST Standard Ref-erence Database Number 69, National Institute of Standards and Technology,Gaithersburg MD, 20899, https://doi.org/10.18434/T4D303, (retrieved June 12,2020)