Semi-smooth Newton methods for nonlinear complementarity formulation of compositional two-phase flow in porous media

https://doi.org/10.1016/j.jcp.2019.109163Get rights and content

Highlights

  • Compositional two-phase flow with phase transitions can be formulated as a nonlinear complementarity problem (NCP).

  • Semi-smooth Newton method using the minimum function for NCP formulation is not robust for highly nonlinear problems and problems with heterogeneous media.

  • Use of the Fischer-Burmeister function accelerates convergence.

  • We develop a new Jacobian smoothing method that is efficient and scalable for a wide range of problems.

Abstract

Simulating compositional multiphase flow in porous media is a challenging task, especially when phase transition is taken into account. The main problem with phase transition stems from the inconsistency of the primary variables such as phase pressure and phase saturation, i.e. they become ill-defined when a phase appears or disappears. Recently, a new approach for handling phase transition has been developed, whereby the system is formulated as a nonlinear complementarity problem (NCP). Unlike the widely used primary variable switching (PVS) method which requires a drastic reduction of the time step size when a phase appears or disappears, this approach is more robust and allows for larger time steps. One way to solve an NCP system is to reformulate the inequality constraints as a non-smooth equation using a complementary function (C-function). Because of the non-smoothness of the constraint equations, a semi-smooth Newton method needs to be developed. In this work, we consider two methods for solving NCP systems used to model multiphase flow: (1) a semi-smooth Newton method for two C-functions: the minimum and the Fischer-Burmeister functions, and (2) a new inexact Newton method based on the Jacobian smoothing method for a smooth version of the Fischer-Burmeister function. We show that the new method is robust and efficient for standard benchmark problems as well as for realistic examples with highly heterogeneous media such as the SPE10 benchmark.

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ϕ(ρlwSl)t+(ρlwqljlh)=0,ϕ(ρlhSl+ρghSg)t+(ρlhql+ρghqg+jlh)=0, where the

Nonlinear complementarity problem

In its simplest form, a nonlinear complementarity problem with respect to a smooth function f:RNRN is to find a vector uRN such thatu0,f(u)0,uTf(u)=0. A slightly more general form of the last equation in (3.1) readsg(u)Tf(u)=0, where g:RNRN 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)

  • A. Abadpour et al.

    Method of negative saturations for modeling two-phase compositional flow with oversaturated zones

    Transp. Porous Media

    (2008)
  • G. Acs et al.

    General purpose compositional model

    Soc. Pet. Eng. J.

    (1985)
  • M. Aganagić

    Newton's method for linear complementarity problems

    Math. Program.

    (1984)
  • ASCEM
  • I. Ben Gharbia et al.

    Nonconvergence of the plain Newton-min algorithm for linear complementarity problems with a P-matrix

    Math. Program.

    (2012)
  • I. Ben Gharbia et al.

    An algorithmic characterization of P-matricity

    SIAM J. Matrix Anal. Appl.

    (2013)
  • I. Ben Gharbia et al.

    An algorithmic characterization of P-matricity II: corrections, refinements, and validation

    SIAM J. Matrix Anal. Appl.

    (2018)
  • A. Bourgeat et al.

    Two-phase partially miscible flow and transport modeling in porous media; application to gas migration in a nuclear waste repository

    Comput. Geosci.

    (2008)
  • A. Bourgeat et al.

    Compositional two-phase flow in saturated-unsaturated porous media: benchmarks for phase appearance/disappearance

  • B. Chen et al.

    A non-interior-point continuation method for linear complementarity problems

    SIAM J. Matrix Anal. Appl.

    (1993)
  • X. Chen et al.

    Smoothing methods and semismooth methods for nondifferentiable operator equations

    SIAM J. Numer. Anal.

    (Jan 2001)
  • 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.

    View full text