Energy landscape analysis for two-phase multi-component NVT flash systems by using ETD type high-index saddle dynamics

Highlights

• The HISD algorithm used for the first time in saddle point calculation in NVT flash.

• The Rosenbrock-Euler ETD format used to reduce the interference of system rigidity.

• At most the 1st-order saddle points in the energy landscape of the 2-component system.

• There is only one solution in the 2-component 2-phase NVT flash problem.

• All the 1st-order saddle points can converge to the equivalent local minima.

Abstract

With the increasingly accurate modeling and simulation demands and techniques, obtaining the knowledge of the multi-component phase equilibrium states holds the bottleneck challenge in the study of multi-phase fluid flow. In order to resolve the shortages in computational efficiency and stability in the conventional iterative flash calculation, a new phase equilibrium prediction method is proposed by solving the saddle point of the multi-phase system. In this paper, the HiSD (High-index saddle point dynamics) algorithm is used for the first time to calculate the saddle points on the energy landscape of the two-phase two-component NVT flash model based on the Peng-Robinson equation of state, and the up-down search algorithm of HiSD is applied to generate the solution landscape of the flash system. The Rosenbrock-Euler ETD (exponential time differencing) format is involved to reduce the interference of system rigidity to the calculation. It can be referred from the numerical analysis that there are at most the $1_{st}$-order saddle points in the energy landscape of the two-component two-phase NVT flash system, and all these saddle points are located on one straight line of the hyperplane, where the energy is equal everywhere. All these $1_{st}$-order saddle points can converge to the same or equivalent local minima, which indicates that the two-component two-phase flash system is a system with only one single solution with physical meanings. In addition, the saddle points also obey a linear relationship and the energy remains the same at different temperatures. Therefore, using the method proposed in this paper, the conventional two-step efforts of phase stability test and phase separation calculation can be simplified. The $1_{st}$-order saddle points of the system can be directly calculated, reducing the need for an initial guess. The local minima can be directly searched through the downward direction of the saddle point, which greatly reduces the calculation amount of phase equilibrium calculations. Furthermore, the minimum states at different temperatures can be calculated in batch by using one certain initial value, which significantly improves the adaptability and reliability for complex engineering problems with drastic temperature changes.

Introduction

Thanks to the rapid development in computing equipment and techniques, the past decades have borne remarkable progress in the more refined modeling and simulation of multi-phase fluid flows with multiple components. Studies on phase distributions and phase transitions have been incorporated into an increasing number of engineering applications to improve the production efficiency [1], ensure the operation safety, [2] and enhance the environmental protection [3]. As a concept defining the aggregated state of homogeneous substances with the same physical and chemical properties, phases are determined and separated by a clear physical interface, crossing which a sudden change in certain thermo properties can be captured [4]. Thus, the observations on the phase behaviors are focused on finding the change of thermodynamic properties, especially modeling the non-analytical free energy function representing phase transitions [5]. The vapor phase and liquid phase are the two major fluid phases existing in the multi-phase flow problems, and the description of the phase transitions between the two phases is usually carried out by flash calculation [6]. By calculating the compositional distributions of each phase at equilibrium, the flash calculation provides the fundamental thermodynamic knowledge for the multi-phase flow modeling and simulation, including the number of phases and the initial phase distributions with physical meanings [7]. In addition, phase equilibrium investigation can also used for independent calculation, analysis, design and optimization of engineering processes and facilities, including reducing unnecessary substances (for example, CO2 and H2O) to purify the production of hydrogen, removing acidic gases (for example, CO2 and H2S) to protect pipelines and equipment from corrosion, determining the amount of inhibitors (for example, methanol and glycol) to avoid the formation of gas hydrates, controlling asphaltene precipitation to enhance flow assurance, and many others [8]. Based on the different constraints on phase equilibrium, flash calculations are usually divided into two categories: NPT flash with constant chemical compositions, pressure, and temperature, and NVT flash with constant chemical compositions, volume, and temperature. Comparing with the NPT flash framework, NVT flash model has well-posed formulation and exhibits a unique pressure-volume relation so that the root selection procedure in PT flash framework can be avoided [9], [10]. At the same time, considering the easier and more intuitive treatment on capillarity, which is the critical mechanism in a number of popular engineering applications including shale gas recovery [11] and microfluidic fuel cell [12], NVT flash has been attracting increasing attentions in the phase equilibrium studies of modern industry.

Despite the enabled prediction on phase equilibrium states by thermodynamic modeling and numerical analysis, the conventional NVT flash calculation methods are usually limited by the necessary step of the phase stability test, which is the guarantee of thermodynamic consistency. Besides, in some of the NVT flash calculation methods, only the local minima on the energy landscape of the system is reached, while the numerical efficiency is challenged in the other methods calculating the global minima of the system [13], [14], [15], [16]. The emerging deep learning algorithm has been successfully applied in accelerating flash calculations [17] with a good performance on both accuracy and efficiency, and both the two steps of phase stability tests and phase splitting calculations have been incorporated in the integrated Thermodynamics-Informed Neural Network, also known as TINN [18]. However, the expensive and limited experimental phase equilibrium data results in the selection of iterative flash calculation results as the ground truth for training and testing in the deep neural network, which places a higher requirement on the reliability of the NVT flash scheme.

A significant breakthrough was proposed in [19], which generated a new thermodynamic analysis approach studying the minimum energy path of the total energy landscape of the two-component two-phase NVT flash system. With a preliminary understanding of the saddle point information and energy landscape of this system, the global minima can be easier to be reached. However, the minimum energy path can only inform the transition route between different local minima on the energy landscape, while the information regarding more saddle points is expected to help generate the complete energy landscape of such system. By calculating the stationary points of one given system, the stable state can be obtained by the 0-order saddle points calculation, but the algorithm stability is always a shortcoming compared with the conventional gradient method converging to the stable state. A number of methods have been proposed on resolving the unstable nature of the saddle points, including the dimer method [22], nudged elastic band method [23], [24], string method [25], [26], [27], gentlest ascent dynamic method [28], [29], minimax method [30] and optimization-based shrinking dimer method [31]. Recently, inspired by the idea of optimization-based shrinking dimer method, a high-order saddle points dynamic method (HiSD) was proposed in [33] to compute saddle points at any orders on the energy landscape of a nonlinear system. This numerical approach can help us to build the solution landscape of the object system, which is a pathway map gathered by all the stationary points on the energy landscape of the system and their connecting paths. This technique has been successfully applied to various fields, such as confined nematic liquid crystal [20], [21], period crystals and quasicrystals [35], diblock copolymers [36] and other applications [37], [38], [39]. Through these applications, the HiSD approach played a very important role to help us to find steady states (local minima) and meta-stable states ($k_{th}$ order saddle points) of the numerical model systematically and build the related solution landscape, which can be used to further understand the relevant physical laws and phenomena.

As mentioned in [19], the NVT flash systems with Peng-Robinson equation of state have ill-conditioned Hessian matrices with large condition numbers due to the strong nonlinearity of the equation of state. Normal HiSD approaches need strict step sizes to make the algorithm converge, which significantly slows down the phase equilibrium calculations and further challenges the requirement of obtaining the phase equilibrium state quickly in engineering practice. In order to build the solution landscape of the NVT flash system efficiently in a numerically-compatible manner, in this study, we need to modify the HiSD and make it possible to handle the challenges caused by the very strong stiffness terms in conventional NVT flash models. The string method can efficiently determine the energy transition pathways for complex systems with smooth energy landscape by evolving strings rather than points in the configuration space. This approach consists of two main steps, the evolution of the string by some ordinary differential equation (ODE) solvers and the reparametrization of the string by interpolation. The accuracy of this approach significantly relies on ODE solvers. Moreover, explicit numerical approaches are required because the real scale of the problem is huge. Some common ODE solvers, like Euler and Runge-Kutta methods, require a very small time step for solving ill-conditioned problems. Thus, the attentions have been focused on the exponential time differencing (ETD) scheme, which can provide a satisfactory accuracy and stability for problems with strong stiffness. The scheme involves the exact integration of the governing equation, followed by an explicit approximation of the time integration of the nonlinear term. For the nonlinear problem with strong stiffness, Hochbruck and Ostermann proposed Rosenbrock-type ETD scheme [40], [41] to treat the stiff nonlinear term without the loss of the original stability. This idea proposed a new numerical routine to establish a good linearization approach for nonlinear problems. In this paper, the idea of Rosenbrock-type ETD scheme is applied to build an ETD-type HiSD scheme to calculate the solution landscape of the two-phase two-component VT flash system.

In this paper, the high order dynamic method will be applied to calculate the saddle points in the NVT flash system, so as to generate the total energy landscape of such system to help proceed the thermodynamic analysis on phase equilibrium states, while the manual intervention in the conventional NVT flash framework is expected to be minimized. The fundamental NVT flash scheme is introduced in Section 2, followed by the construction of the high order dynamic method searching for the saddle points. Thermodynamic analysis on the two-component two-phase NVT flash system is presented in Section 3, showing the example saddle points calculations in a fluid mixture, as well as the phase equilibrium analysis based on the saddle points and energy landscape. Concluding remarks are organized in Section 4, and a simple introduction on the Peng-Robinson equation of state, a common-used model in thermodynamic analysis regarding hydrocarbons, is included in the appendix.