Ordering-based nonlinear solver for fully-implicit simulation of three-phase flow

The Fully Implicit method (FIM) is often the method of choice for the temporal discretization of the partial differential equations governing multiphase flow in porous media. The FIM involves solving large coupled systems of nonlinear algebraic equations. Newton-based methods, which are employed to solve the nonlinear systems, can suffer from convergence problems—this is especially true for large time steps in the presence of highly nonlinear flow physics. To overcome such convergence problems, the time step is usually reduced, and the Newton steps are restarted from the solution of the previous (converged) time step. Recently, potential ordering and the reduced-Newton method were used to solve immiscible three-phase flow in the presence of buoyancy and capillary effects (e.g., Kwok and Tchelepi, J. Comput. Phys. 227(1), 706–727 2007). Here, we improve the robustness of the potential-based ordering method in the presence of gravity. Furthermore, we also extend this nonlinear approach to interphase mass transfer. Our algorithm deals effectively with mass transfer between the liquid and gas phases, including phase disappearance (e.g., gas going back in solution) and reappearance (e.g., gas coming out of solution and forming a separate phase), as a function of pressure and composition. Detailed comparisons of the robustness and efficiency of the potential-based solver with state-of-the-art nonlinear/linear solvers are presented for immiscible two-phase (Dead-Oil), Black-Oil, and compositional problems using heterogeneous models. The results show that for large time steps, our nonlinear ordering-based solver reduces the number of nonlinear iterations significantly, which leads to gains in the overall computational cost.     

A parallel, implicit, cell‐centered method for two‐phase flow with a preconditioned Newton–Krylov solver

A new parallel solution technique is developed for the fully implicit three‐dimensional two‐phase flow model. An expandedcell‐centered finite difference scheme which allows for a full permeability tensor is employed for the spatial discretization, and backwardEuler is used for the time discretization. The discrete systems are solved using a novel inexact Newton method that reuses the Krylov information generated by the GMRES linear iterative solver. Fast nonlinear convergence can be achieved by composing inexact Newton steps with quasi‐Newton steps restricted to the underlying Krylov subspace. Furthermore, robustness and efficiency are achieved with a line‐search backtracking globalization strategy for the nonlinear systems and a preconditioner for each coupled linear system to be solved. This inexact Newton method also makes use of forcing terms suggested by Eisenstat and Walker which prevent oversolving of the Jacobian systems. The preconditioner is a new two‐stage method which involves a decoupling strategy plus the separate solutions of both nonwetting‐phase pressure and saturation equations. Numerical results show that these nonlinear and linear solvers are very effective.     

Nonlinear solver for three-phase transport problems based on approximate trust regions

Implicit transport solvers used in reservoir simulation can take longer time steps than explicit solvers, but for long time steps, the commonly used Newton-Raphson’s method will often fail to converge. The convergence issues may manifest themselves as oscillating residuals even though the implicit discretization itself is stable. This behavior occurs because the fractional flow-type flux functions often change between convex and concave during long time steps, resulting in multiple contraction regions for the Newton-Raphson solver. The common strategy to overcome this is to set limits on the saturation changes during the nonlinear iteration, but such a limit has to be determined on a case by case basis, excess iterations may be required, and practical convergence is not guaranteed for a given problem. Previous work on this problem by multiple authors has resulted in solvers based on trust regions, where unconditional convergence can be obtained for incompressible two-phase flow provided a priori analytical knowledge of the flux function exists. The goal of our work is to extend this methodology to a solver where inflection points demarking the different contraction regions do not need to be explicitly known. Instead, these values are estimated during the solution process, giving improved convergence by a local computation for each interface in the simulation model. By systematically reducing updates over regions known to produce convergence issues, it is possible to greatly reduce the computational expense, making the same formulation suitable for an arbitrary number of components. We present a series of numerical results, including arbitrary time-step lengths for two and three-phase gravity segregation, as well as three-dimensional gas and water injection problems with wells and a mixture of both viscous and gravity-dominated flow regimes. The test cases are a systematic validation on a wide variety of both analytical and tabulated relative permeability curves.     

Accelerating the solution of linear systems appearing in two-phase reservoir simulation by the use of POD-based deflation methods

We explore and develop a Proper Orthogonal Decomposition (POD)-based deflation method for the solution of ill-conditioned linear systems, appearing in simulations of two-phase flow through highly heterogeneous porous media. We accelerate the convergence of a Preconditioned Conjugate Gradient (PCG) method achieving speed-ups of factors up to five. The up-front extra computational cost of the proposed method depends on the number of deflation vectors. The POD-based deflation method is tested for a particular problem and linear solver; nevertheless, it can be applied to various transient problems, and combined with multiple solvers, e.g., Krylov subspace and multigrid methods.

     

A one-dimensional local discontinuous Galerkin Richards’ equation solution with dual-time stepping

We present a compact, high-order Richards’ equation solver using a local discontinuous Galerkin finite element method in space and a dual-time stepping method in time. Dual-time stepping methods convert a transient problem to a steady state problem, enabling direct evaluation of residual terms and resolve implicit equations in a step-wise manner keeping the method compact and amenable to parallel computing. Verification of our solver against an analytical solution shows high-order error convergence and demonstrates the solvers ability to maintain high accuracy using low spatial resolution; the method is robust and accurately resolves numerical solutions with time steps that are much larger than what is normally required for lower-order implicit schemes. Resilience of our solver (in terms of nonlinear convergence) is demonstrated in ponded infiltration into homogeneous and layered soils, for which HYDRUS-1D solutions are used as qualitative references to gauge performance of two slope limiting schemes.

     

Multiscale finite-volume formulation for multiphase flow in porous media: black oil formulation of compressible,three-phase flow with gravity

Most practical reservoir simulation studies are performed using the so-called black oil model, in which the phase behavior is represented using solubilities and formation volume factors. We extend the multiscale finite-volume (MSFV) method to deal with nonlinear immiscible three-phase compressible flow in the presence of gravity and capillary forces (i.e., black oil model). Consistent with the MSFV framework, flow and transport are treated separately and differently using a sequential implicit algorithm. A multiscale operator splitting strategy is used to solve the overall mass balance (i.e., the pressure equation). The black-oil pressure equation, which is nonlinear and parabolic, is decomposed into three parts. The first is a homo geneous elliptic equation, for which the original MSFV method is used to compute the dual basis functions and the coarse-scale transmissibilities. The second equation accounts for gravity and capillary effects; the third equation accounts for mass accumulation and sources/ sinks (wells). With the basis functions of the elliptic part, the coarse-scale operator can be assembled. The gravity/capillary pressure part is made up of an elliptic part and a correction term, which is computed using solutions of gravity-driven local problems. A particular solution represents accumulation and wells. The reconstructed fine-scale pressure is used to compute the fine-scale phase fluxes, which are then used to solve the nonlinear saturation equations. For this purpose, a Schwarz iterative scheme is used on the primal coarse grid. The framework is demonstrated using challenging black-oil examples of nonlinear compressible multiphase flow in strongly heterogeneous formations.     

Hybrid discretization of multi-phase flow in porous media in the presence of viscous,gravitational, and capillary forces

Multi-phase flow in porous media in the presence of viscous, gravitational, and capillary forces is described by advection diffusion equations with nonlinear parameters of relative permeability and capillary pressures. The conventional numerical method employs a fully implicit finite volume formulation. The phase-potential-based upwind direction is commonly used in computing the transport terms between two adjacent cells. The numerical method, however, often experiences non-convergence in a nonlinear iterative solution due to the discontinuity of transmissibilities, especially in transition between co-current and counter-current flows. Recently, Lee et al. (Adv. Wat. Res. 82, 27–38, 2015) proposed a hybrid upwinding method for the two-phase transport equation that comprises viscous and gravitational fluxes. The viscous part is a co-current flow with a one-point upwinding based on the total velocity and the buoyancy part is modeled by a counter-current flow with zero total velocity. The hybrid scheme yields C¹-continuous discretization for the transport equation and improves numerical convergence in the Newton nonlinear solver. Lee and Efendiev (Adv. Wat. Res. 96, 209–224, 2016) extended the hybrid upwind method to three-phase flow in the presence of gravity. In this paper, we present the hybrid-upwind formula in a generalized form that describes two- and three-phase flows with viscous, gravity, and capillary forces. In the derivation of the hybrid scheme for capillarity, we note that there is a strong similarity in mathematical formulation between gravity and capillarity. We thus greatly utilize the previous derivation of the hybrid upwind scheme for gravitational force in deriving that for capillary force. Furthermore, we also discuss some mathematical issues related to heterogeneous capillary domains and propose a simple discretization model by adapting multi-valued capillary pressures at the end points of capillary pressure curves. We demonstrate this new model always admits a consistent solution that is within the discretization error. This new generalized hybrid scheme yields a discretization method that improves numerical stability in reservoir simulation.     

Shock trains on a planar beach: quasi-analytical and fully numerical solutions

This study, part of the Special Issue dedicated to the 70th anniversary of Professor Efim Pelinovsky, focuses on a topic that has been central in Professor Pelinovsky’s research, i.e. the analytical and numerical modelling of shallow water waves. We specifically focus on the evolution of trains of shock waves on a planar beach. Antuono (J Fluid Mech 658:166–187, 2011) has, for the first time, proposed a quasi-analytical solution for a train of shock waves forced by a constant Riemann invariant. The present contribution clarifies the validity of such solution and its value for benchmarking nonlinear shallow water equation solvers. Hence, the same tests of Antuono (J Fluid Mech 658:166–187, 2011) have been run by means of the solver of Brocchini et al. (Coast Eng 43(2):105–129, 2001) revealing surprisingly and reassuring good agreements. This provides significant support to the mentioned analytical solution and allows to critically analyse the eventual discrepancies, due to the practicalities of running numerical shallow water solutions (e.g. influence of the boundary conditions, of the numerical resolution, etc.).     

Coupled hydromechanical‐fracture simulations of nonplanar three‐dimensional hydraulic fracture propagation

This paper presents an algorithm and a fully coupled hydromechanical‐fracture formulation for the simulation of three‐dimensional nonplanar hydraulic fracture propagation. The propagation algorithm automatically estimates the magnitude of time steps such that a regularized form of Irwin's criterion is satisfied along the predicted 3‐D fracture front at every fracture propagation step. A generalized finite element method is used for the discretization of elasticity equations governing the deformation of the rock, and a finite element method is adopted for the solution of the fluid flow equation on the basis of Poiseuille's cubic law. Adaptive mesh refinement is used for discretization error control, leading to significantly fewer degrees of freedom than available nonadaptive methods. An efficient computational scheme to handle nonlinear time‐dependent problems with adaptive mesh refinement is presented. Explicit fracture surface representations are used to avoid mapping of 3‐D solutions between generalized finite element method meshes. Examples demonstrating the accuracy, robustness, and computational efficiency of the proposed formulation, regularized Irwin's criterion, and propagation algorithm are presented.     

Analytical and numerical solutions for carbonated waterflooding

We develop a Riemann solver for transport problems including geochemistry related to oil recovery. The example considered here concerns one-dimensional incompressible flow in porous media and the transport for several chemical components, namely H₂O, H⁺, OH^?, CO₂, \(\text {CO}_{3}^{2-}\), \(\text {HCO}_{3}^{-}\), and decane; they are in chemical equilibrium in the aqueous and oleic phases, leading to mass transfer of CO₂ between the oleic and aqueous phases. In our ionic model, we employ equations with zero diffusion coefficients. We do so because it is well known that for upscaled equations, the convection terms dominate the diffusion terms. The Riemann solution for this model can therefore be applied for upscaled transport processes in enhanced oil recovery involving geochemical aspects. In our example, we formulate the conservation equations of hydrogen, oxygen, hydrogen, and decane, in which we substitute regression expressions that are obtained by geochemical software. This can be readily done because Gibbs phase rule together with charge balance shows that all compositions can be rewritten in terms of a single composition, which we choose to be the hydrogen ion concentration (p H). In our example, we use the initial and boundary conditions for the carbonated aqueous phase injection in an oil reservoir containing connate water with some carbon dioxide. We compare the Riemann solution with a numerical solution, which includes capillary and diffusion effects. The significant new contribution is the effective Riemann solver we developed to obtain solutions for oil recovery problems including geochemistry and a variable total Darcy velocity, a situation in which fractional flow theory does not readily apply. We thus obtain an accurate solution for a carbonated waterflood, which elucidates some mechanisms of low salinity carbonated waterflooding.     

Preconditioned IDR(s) iterative solver for non‐symmetric linear system associated with FEM analysis of shallow foundation

Non‐associated flow rule is essential when the popular Mohr–Coulomb model is used to model nonlinear behavior of soil. The global tangent stiffness matrix in nonlinear finite element analysis becomes non‐symmetric when this non‐associated flow rule is applied. Efficient solution of this large‐scale non‐symmetric linear system is of practical importance. The standard Krylov solver for a non‐symmetric solver is Bi‐CGSTAB. The Induced Dimension Reduction [IDR(s)] solver was proposed in the scientific computing literature relatively recently. Numerical studies of a drained strip footing problem on homogenous soil layer show that IDR(s = 6) is more efficient than Bi‐CGSTAB when the preconditioner is the incomplete factorization with zero fill‐in of global stiffness matrix K_ep (ILU(0)‐K_ep). Iteration time is reduced by 40% by using IDR(s = 6) with ILU(0)‐K_ep. To further reduce computational cost, the global stiffness matrix K_ep is divided into two parts. The first part is the linear elastic stiffness matrix K_e, which is formed only once at the beginning of solution step. The second part is a low‐rank matrix Δ, which is re‐formed at each Newton–Raphson iteration. Numerical studies show that IDR(s = 6) with this ILU(0)‐K_e preconditioner is more time effective than IDR(s = 6) with ILU(0)‐K_ep when the percentage of yielded Gauss points in the mesh is less than 15%. The total computation time is reduced by 60% when all the recommended optimizing methods are used. Copyright © 2013 John Wiley & Sons, Ltd.     

A finite volume / discontinuous Galerkin method for the advective Cahn–Hilliard equation with degenerate mobility on porous domains stemming from micro-CT imaging

A numerical method is formulated for the solution of the advective Cahn–Hilliard (CH) equation with constant and degenerate mobility in three-dimensional porous media with non-vanishing velocity on the exterior boundary. The CH equation describes phase separation of an immiscible binary mixture at constant temperature in the presence of a conservation constraint and dissipation of free energy. Porous media / pore-scale problems specifically entail images of rocks in which the solid matrix and pore spaces are fully resolved. The interior penalty discontinuous Galerkin method is used for the spatial discretization of the CH equation in mixed form, while a semi-implicit convex–concave splitting is utilized for temporal discretization. The spatial approximation order is arbitrary, while it reduces to a finite volume scheme for the choice of element-wise constants. The resulting nonlinear systems of equations are reduced using the Schur complement and solved via inexact Newton’s method. The numerical scheme is first validated using numerical convergence tests and then applied to a number of fundamental problems for validation and numerical experimentation purposes including the case of degenerate mobility. First-order physical applicability and robustness of the numerical method are shown in a breakthrough scenario on a voxel set obtained from a micro-CT scan of a real sandstone rock sample.     

A fully-coupled discontinuous Galerkin method for two-phase flow in porous media with discontinuous capillary pressure

In this paper, we formulate and test numerically a fully-coupled discontinuous Galerkin (DG) method for incompressible two-phase flow with discontinuous capillary pressure. The spatial discretization uses the symmetric interior penalty DG formulation with weighted averages and is based on a wetting-phase potential/capillary potential formulation of the two-phase flow system. After discretizing in time with diagonally implicit Runge-Kutta schemes, the resulting systems of nonlinear algebraic equations are solved with Newton’s method and the arising systems of linear equations are solved efficiently and in parallel with an algebraic multigrid method. The new scheme is investigated for various test problems from the literature and is also compared to a cell-centered finite volume scheme in terms of accuracy and time to solution. We find that the method is accurate, robust, and efficient. In particular, no postprocessing of the DG velocity field is necessary in contrast to results reported by several authors for decoupled schemes. Moreover, the solver scales well in parallel and three-dimensional problems with up to nearly 100 million degrees of freedom per time step have been computed on 1,000 processors.     

Finite volume approximation for an immiscible two-phase flow in porous media with discontinuous capillary pressure

We consider an immiscible incompressible two-phase flow in a porous medium composed of two different rocks so that the capillary pressure field is discontinuous at the interface between the rocks. This leads us to apply a concept of multivalued phase pressures and a notion of weak solution for the flow which have been introduced in Cancès and Pierre (SIAM J Math Anal 44(2):966–992, 2012). We discretize the problem by means of a numerical algorithm which reduces to a standard finite volume scheme in each rock and prove the convergence of the approximate solution to a weak solution of the two-phase flow problem. The numerical experiments show in particular that this scheme permits to reproduce the oil-trapping phenomenon.     

Cu (II) removal intensification using Fe<Subscript>3</Subscript>O<Subscript>4</Subscript> nanoparticles under inert gas and magnetic field in a microchannel

In the present study, Cu (II) ions removal from aqueous solution was intensified by exciting magnetic nanoparticles under inert gas, magnetic field and combination of these two mixing methods in a T-type microchannel. The flow patterns and liquid–liquid two-phase mass transfer were studied in three different magnet distances from mixing channel (3, 6 and 10 mm) and also in the presence of different inert gas flow rates (1, 3 and 5 mL/min). Depending on the mixing method and the flow rate of both phases, several distinct flow patterns were observed including slugs, droplet, parallel and dispersed flows. The performances of mixing techniques for mass transfer enhancement based on relative removal efficiency ratio (λ) and mass transfer coefficient ratio (γ) were compared with simple layout (without nanoparticles, magnetic field and inert gas). The results showed that simultaneous using of inert gas and magnetic field can drive the nanoparticles as mixer. Liquid–liquid mass transfer with 27–62% enhancement in E and 235–285% in K _L a compared with plain one was observed.     

Finite element modelling of three-phase flow in deforming saturated oil reservoirs

A numerical simulation is presented for three-dimensional three-phase fluid flow in a deforming saturated oil reservoir. The mathematical formulation describes a fully coupled governing equation systen which consists of the equilibrium and continuity equations for three immiscible fluids flowing in a porous medium. An elastoplastic soil model, based on a Mohr–Coulomb yield surface, is used. The finite element method is applied to obtain simultaneous solutions to the governing equations where displacement and fluid pressures are the primary unknowns. The final discretized equations are solved by a direct solver using fully implicit procedures. The developed model is used to investigate the problems of three-phase fluid flow in a deforming saturated oil reservoir.     

A finite element solution of a fully coupled implicit formulation for reservoir simulation

An iterative method is presented for solving a fully coupled and implicit formulation of fluid flow in a porous medium. The mathematical model describes a set of fully coupled three-phase flow of compressible and immiscible fluids in a saturated oil reservoir. The finite element method is applied to obtain the simultaneous solution (SS) for the resulting highly non-linear partial differential equations where fluid pressures are the primary unknowns. The final discretized equations are solved iteratively by using a fully implicit numerical scheme. Several examples, illustrating the use of the present model, are described. The increased stability achieved with this scheme has permitted the use of larger time steps with smaller material balance errors.     

An adjoint-based optimization method for gas production in shale reservoirs

In this work, we consider a new model for flow in a multiporosity shale gas reservoir constructed within the framework of an upscaling procedure where hydraulic fractures are treated as (\(n-1\)) interfaces (\(n=2,3\)). Within this framework, the hydrodynamics is governed by a new pressure equation in the shale matrix which is treated as a homogenized porous medium composed of organic matter (kerogen aggregates with nanopores) and inorganic impermeable solid (clay, calcite, quartz) separated from each other by a network of interparticle pores of micrometer size. The solution of the pressure equation is strongly influenced by the constitutive response of the retardation parameter and effective hydraulic conductivity where the former incorporates gas adsorption/desorption in the nanopores of the kerogen. By focusing our analyses on this nonlinear diffusion equation in the domain occupied by the shale matrix, an optimization strategy seated on the adjoint sensitivity method is developed to minimize a cost functional related to gas production and net present value in a single hydraulic fracture. The gradient of the objective functional computed with the adjoint formulation is explored to update the controlled pressure drop aiming to optimize production in a given window of time. The combination of the direct approach and gradient-based optimization using the adjoint formulation leads to the construction of optimal production scenarios under controlled pressure decline in the well. Numerical simulations illustrate the potential of the methodology proposed herein in optimizing gas production.     

Numerical Homogenization of Two-Phase Flow in Porous Media

Homogenization has proved its effectiveness as a method of upscaling for linear problems, as they occur in single-phase porous media flow for arbitrary heterogeneous rocks. Here we extend the classical homogenization approach to nonlinear problems by considering incompressible, immiscible two-phase porous media flow. The extensions have been based on the principle of preservation of form, stating that the mathematical form of the fine-scale equations should be preserved as much as possible on the coarse scale. This principle leads to the required extensions, while making the physics underlying homogenization transparent. The method is process-independent in a way that coarse-scale results obtained for a particular reservoir can be used in any simulation, irrespective of the scenario that is simulated. Homogenization is based on steady-state flow equations with periodic boundary conditions for the capillary pressure. The resulting equations are solved numerically by two complementary finite element methods. This makes it possible to assess a posteriori error bounds.