In this work, a new method for PHN phase equilibrium calculations is presented. The problem is formulated as an unconstrained maximization of the entropy with respect to mole numbers. Enthalpy is a dependent variable and its dependence on mole numbers is taken into account by solving a nonlinear equation (the enthalpy balance equation) for temperature at each iteration level. The calculation framework is similar to that in PT flash. A combined successive substitution - Newton method is used, with mole numbers or natural logarithms of equilibrium constants as independent variables. The iterations in the PHN flash are “PT-like” (for a sequence of temperatures converging to the final equilibrium temperature).
A new successive substitution iteration (SSI) method is proposed for PHN flash calculations. The iterative sequence consists of two inner loops, in which the Rachford-Rice equation is solved first (using convex transformations and narrow solution windows) for the vapor phase fraction and then the enthalpy constraint equation is solved (using Brent’s method) for temperature. The equilibrium constants are updated in the outer loop. An ascent direction is guaranteed in the SSI method and an efficient line search procedure ensures an increase in entropy at each iteration. If non-feasible conditions are encountered at the initial guess or during iterations, a transition temperature is calculated, and the sub-phase method is used; one of the equilibrium phases is artificially split into two sub-phases having the same composition, in such a way that the enthalpy constraint is honored. As a result, a pseudo-surface of the entropy is constructed in the non-feasible domain, allowing the iterative method to advance and eventually reach the feasible domain closer to the solution. The new SSI method is much simpler than various versions of the so-called direct substitution (DS) method. However, the SSI method in PHN phase equilibrium is not as robust as its PT counterpart (this feature was recognized to be common to other major specifications, such as VTN or UVN) and its limitations are discussed. Several unfavorable convergence behavior patterns are identified; i.e., the SSI method is not able to reduce the gradient norm below a certain (usually low) value to the desired accuracy or exhibit slow oscillating convergence. In fact, fortunately, SSI does not need to be effectively converged in practice; the purpose is to bridge between a potentially poor initial guess and the domain near the solution where Newton iterations exhibit unproblematic quadratic convergence.
In the second-order Newton method, a modified Cholesky factorization is used (with a diagonal correction of the Hessian matrix when it becomes indefinite), with a line search procedure to guarantee an increase in entropy at each iteration. Analyses of the block structure of the Hessian matrices for different choices of independent variables (natural variables, mole numbers and temperature, mole numbers in PHN and PT flashes) reveal the relation between all these methods and suggest how to easily modify existing codes to incorporate the new procedures. The elements of the Hessian matrix in the proposed method are related to those of the Hessian in the conventional PT flash by a rank-one term, containing first-order partial derivatives of enthalpy with respect to mole numbers and temperature. Criteria adapted to PHN flash are used for switching from SSI to Newton iterations.
The proposed method is tested for several mixtures of various complexities, including binary mixtures, synthetic mixtures, and reservoir fluids, with emphasis on near-critical, narrow boiling cases and low-temperature conditions and proved to be rapid and robust. A detailed convergence analysis is provided for both SSI and Newton methods. In the binary case, contour plots of entropy in the plane defined by mole numbers allow examination of the convergence paths of different methods. Several domains in the P-H plane are identified where the DS method and/or starting directly with Newton method (or after some partial Newton iterations) struggle to converge or diverge, while the proposed SSI-Newton converges, with only a few Newton iterations required after switching. In this work, it is the first time that i) successive substitution iterations are formulated in PHN flash calculations, similar to the widely used SSI method in conventional PT flashes and ii) the direct maximization of entropy using a Newton method with mole numbers as independent variables is presented and analyzed in detail. The proposed PHN phase equilibrium calculation method is currently being implemented in a compositional reservoir simulator for thermal enhanced oil recovery and carbon dioxide storage.