Semi-smooth Newton methods for nonlinear complementarity formulation of compositional two-phase flow in porous media☆
Introduction
Multiphase flow is a critical process in a wide range of hydrodynamic phenomena, including carbon sequestration, reservoir simulation, and groundwater remediation. For simulation, it would be ideal to have a complete knowledge of the state and composition of the fluid phases in the flow. However, this is not an easy task given the complex physics involved. Some of the most important effects that need to be taken into consideration include capillarity, miscibility, and especially phase transitions. If not handled correctly, these effects can introduce nonphysical numerical oscillations in computational solutions of the strongly nonlinear system of partial differential equations.
Phase transitions have posed a major challenge for multiphase, multi-component models since the 1980s. If not handled correctly, they can cause numerical oscillations in solutions of these models, making such solutions physically inconsistent and unusable. There have been many attempts to address the problems with phase transitions and determine the correct local thermodynamic state for compositional multiphase flow. In general, most of these can be classified into two common classes of methods: flash calculation [2], [16], [17], [21], [34], [43] and primary variables switching (PVS) [26], [41]. Flash calculation computes the local thermodynamic state from the overall mass of the individual components. While this method is stable with regard to determining the thermodynamic state, it tends to be inefficient because it requires solution of a large nonlinear system of equations at each time step (in addition to solution of the linearized systems) to recover all the thermodynamic quantities of interest. The second class, PVS, involves adapting the primary variables to the thermodynamic constraints locally. The idea is that whenever phase transitions occur, physical variables that are physically inconsistent (indicated by negative saturation, for example) are switched to well-defined quantities. The governing equations related to those variables are also modified accordingly. Although this approach is locally more efficient than flash calculations, it suffers from irregular convergence behavior in the nonlinear solve, which is typically addressed by substantial reduction in time step size [20]. This feature is not desirable for simulations over a long period of time, usually encountered in groundwater remediation or transport of nuclides in a nuclear waste repository. In addition to flash calculations and PVS, there are other formulations to handle phase transitions such as negative saturation [1], and introduction of persistent primary variables [9], [33], [35].
Recently, a new approach has been developed for handling the phase transitions by formulating the system of equations as a nonlinear complementarity problem (NCP) [8], [30], [32]. In contrast to PVS, NCP has the advantage that the set of primary variables is consistent throughout the simulation, and no primary variable switching is needed. NCP requires a complementary function, referred to as a C-function, employed to rewrite the inequality constraints for the thermodynamic state as a non-smooth equation, which requires a semi-smooth Newton method [3], [37], [38] to solve. Most of the previous work in multiphase flow using the NCP approach employs the minimum function as the C-function due to its simplicity for implementation and the fact that it is piecewise linear with respect to the arguments. Even though the semi-smooth Newton method applied to the NCP using the minimum function as a C-function is observed to have quadratic convergence for simple problems in porous media (see [8]), we find that it exhibits poor convergence and even diverges for standard benchmark problems, as well as examples considered in this work. An alternative to the minimum function is the Fischer-Burmeister function, which has recently been employed as the C-function for NCP formulation of incompressible two-phase flow in [42]. As we will show, this choice of C-function can help mitigate the lack of robustness observed in using the Newton-min algorithm for NCP formulation of compositional two-phase flow with phase transitions. We then draw on this experience and develop a new method for the nonlinear solve based on a smooth version of the Fischer-Burmeister function. Our method can be considered a variant of the Jacobian smoothing method summarized in [25]. Compared to the non-smooth approaches that use the minimum and the Fischer-Burmeister functions, our new method is more robust and efficient for problems with highly heterogeneous media, and it also scales optimally with problem size.
We consider a two-phase, two-component system with phase transitions as our model problem. We describe this model in detail in section 2, and in section 3, we describe the NCP formulation for it. We briefly review the semi-smooth Newton framework and introduce our new algorithm in section 4. In section 5, several numerical tests are presented that demonstrate the robustness and scalability of the new algorithm. Some concluding remarks as well as future work are presented in section 6.
Section snippets
Governing equations
In this work, we consider a simplified two-phase two-component model with phase transitions. The phases consist of liquid and gas, and the components are water and hydrogen. We also make the following assumptions: (1) water does not vaporize so the gas phase contains only hydrogen, and (2) the amount of hydrogen dissolved in the liquid phase is small. For the two components, the mass conservation equations read where the
Nonlinear complementarity problem
In its simplest form, a nonlinear complementarity problem with respect to a smooth function is to find a vector such that A slightly more general form of the last equation in (3.1) reads where is another smooth function. As we have mentioned in section 2, for the solution of (2.1) and (2.2) to be valid, the pressure, saturation, and hydrogen concentration in the liquid phase must satisfy the constraints in (2.6) and (2.7). These conditions
Solution algorithm
We consider solving the coupled system consisting of (2.1), (2.2), and (3.3) fully implicitly. We use a cell-centered finite volume method for spatial discretization, as it is a natural way to preserve the mass conservation property of the balance (2.1) and (2.2). In addition, it can deal with the case of discontinuous permeability coefficients, and it is relatively straightforward to implement. For the time domain, we employ the backward Euler method to avoid a CFL stability restriction on the
Numerical results
In this section, we describe the results of numerical experiments for solving the NCP systems using both the semi-smooth Newton's approach using the minimum and the Fischer-Burmeister functions, and the Jacobian smoothing method with the smooth Fischer-Burmeister function for the NCP formulation. All of these methods are implemented in Amanzi, a parallel open-source multi-physics C++ code developed as a part of the ASCEM project [4]. Although Amanzi was first designed for simulation of
Conclusions
In this work, we have developed a new Jacobian smoothing method based on the smooth Fischer-Burmeister function to solve the discrete nonlinear systems resulting from the fully implicit discretization of the NCP formulation for compositional multiphase flow in porous media with phase transitions. Additionally, we performed various numerical experiments to compare our method with a semi-smooth Newton approach for two choices of C-function: the minimum and the Fischer-Burmeister functions. The
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
References (43)
- et al.
Gas phase appearance and disappearance as a problem with complementarity constraints
Math. Comput. Simul.
(2014) - et al.
Algebraic multigrid preconditioners for two-phase flow in porous media with phase transitions
Adv. Water Resour.
(2018) - et al.
Comparison of various formulations of three-phase flow in porous media
J. Comput. Phys.
(1997) - et al.
Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media
Adv. Water Resour.
(1995) - et al.
A new approach for phase transitions in miscible multi-phase flow in porous media
Adv. Water Resour.
(2011) - et al.
Results of the MoMaS benchmark for gas phase appearance and disappearance using generalized MHFE
Adv. Water Resour.
(2014) The isothermal flash problem. Part II. Phase-split calculation
Fluid Phase Equilib.
(1982)- et al.
Uncertainties in practical simulation of CO2 storage
Int. J. Greenh. Gas Control
(2012) - et al.
On the selection of primary variables in numerical formulation for modeling multiphase flow in porous media
J. Contam. Hydrol.
(2001) - et al.
Nonlinearly preconditioned semismooth Newton methods for variational inequality solution of two-phase flow in porous media
J. Comput. Phys.
(2017)
Method of negative saturations for modeling two-phase compositional flow with oversaturated zones
Transp. Porous Media
General purpose compositional model
Soc. Pet. Eng. J.
Newton's method for linear complementarity problems
Math. Program.
Nonconvergence of the plain Newton-min algorithm for linear complementarity problems with a P-matrix
Math. Program.
An algorithmic characterization of P-matricity
SIAM J. Matrix Anal. Appl.
An algorithmic characterization of P-matricity II: corrections, refinements, and validation
SIAM J. Matrix Anal. Appl.
Two-phase partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository
Comput. Geosci.
Compositional two-phase flow in saturated-unsaturated porous media: benchmarks for phase appearance/disappearance
A non-interior-point continuation method for linear complementarity problems
SIAM J. Matrix Anal. Appl.
Smoothing methods and semismooth methods for nondifferentiable operator equations
SIAM J. Numer. Anal.
Cited by (0)
- ☆
This work is funded by the U.S. Department of Energy Office of Advanced Scientific Computing Research, Applied Mathematics program, under Award Number DE-SC0009301.