(126h) A Hybrid Newton-Trust-Region Method in Multiphase Isenthalpic Flash Calculations for CO2 Storage Simulations

Carbon dioxide storage in depleted reservoirs/saline aquifers is considered as one of the most effective solutions for carbon capture and storage. A successful carbon storage implementation requires robust and accurate compositional flow/seepage simulation for the fluids containing a high percentage of CO₂. In such simulations, conventional isobaric-isothermal (PT) flash calculation may encounter computational difficulties in simulating the so-called narrow-boiling behavior. In other words, in a narrow temperature interval, the total enthalpy and phase fractions change significantly, for both CO₂-hydrocarbon mixtures and pure CO₂. In these extreme cases, a PT flash calculation is not capable of capturing the vapor-liquid equilibria of pure CO₂ since such equilibria are implicitly represented through a two-phase line in the PT plane. Thus, in compositional simulation, PH flash calculations may be applied in some cases instead of PT flashes to conduct CO₂-storage simulations in the correct manner.

PH flash calculations employ pressure and enthalpy as the specification variables. Due to the reason that the total enthalpy is an extensive variable implicitly related to an equation of state, the solution of the PH flash is always found by satisfying in a certain manner the enthalpy constraint (Michelsen, 1987). In the literature, the partial Newton and full Newton methods solve the enthalpy equality together with the Rachford-Rice equations and fugacity equalities. Direct maximization of the entropy is naturally using a second-order optimization method with mole numbers and phase enthalpies as independent variables. However, the solution temperature is indirectly obtained by a one-dimensional search to honor the enthalpy constraint. Another practical solution method is the nested approach, which treats the enthalpy equality separately in the outer loop, while a well-designed PT flash calculation algorithm is executed in the inner loop to minimize the Gibbs free energy. In this way, the isenthalpic flash problem is reduced to a univariate search in terms of temperature, which can be formulated as either a root-finding or a maximization problem. The convergence is guaranteed in the nested approach, provided a robust PT solver is available. Several authors used the nested approach as the backup solution method in the rare case when their derivative-based methods fail. On the other hand, the nested approach is relatively inefficient due to multiple calls of the PT solver.

A second-order method is of great importance in solving flash calculation problems. In the case of the PT flash, the successive substitution (SSI) method is applied first, and we switch to a second-order method for a final convergence. In terms of the PH flash, a series of direct substitutions are used before applying some second-order methods. Michelsen’s second-order method aims to minimize a penalty function (i.e., the so-called augmented Q function) involving both the enthalpy equality and Gibbs free energy. The saddle-point of the Q-function corresponds to a minimum of the augmented Q function. The latter is minimized with respect to temperature and mole numbers. However, neither Michelsen (1987), nor later research on the PH flash covered detailed formulations and testing for this method. All partial derivatives required in Michelsen’s second-order method, as well as details on the implementation and results of the application of the corresponding numerical algorithm to solve two- and three-phase isenthalpic flash problems are given here for the first time. The proposed hybrid algorithm consists in Newton iterations (with line search to ensure a decrease of the objective function) combined with a trust-region method (Petitfrere and Nichita, 2015) when the Hessian is not positive definite (that is, the Newton fails to provide a descent direction).

We conduct extensive convergence tests using several representative fluids, including examples showing a narrow-boiling behavior. The initial guesses of the independent variables (temperature and mole numbers) are obtained from ideal K-values for the two-phase cases and from stability testing for the three-phase cases. The proposed calculation algorithm proved to be robust and efficient. Convergence is obtained in the quasi-totality of test points; some outliers require the use of the extremely robust, but slower nested approach. The results demonstrate that the proposed algorithm is an ideal candidate for robust and fast isenthalpic flash routines in complex compositional simulations, considering two- and three-phase equilibria and including narrow-boiling behavior.