Smoothed particle hydrodynamics and its applications for multiphase flow and reactive transport in porous media

Smoothed particle hydrodynamics (SPH) is a Lagrangian method based on a meshless discretization of partial differential equations. In this review, we present SPH discretization of the Navier-Stokes and advection-diffusion-reaction equations, implementation of various boundary conditions, and time integration of the SPH equations, and we discuss applications of the SPH method for modeling pore-scale multiphase flows and reactive transport in porous and fractured media.     

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.     

Discontinuous Galerkin method for two-component liquid–gas porous media flows

We consider two-component (typically, water and hydrogen) compressible liquid–gas porous media flows including mass exchange between phases possibly leading to gas-phase (dis)appearance, as motivated by hydrogen production in underground repositories of radioactive waste. Following recent work by Bourgeat, Jurak, and Smaï, we formulate the governing equations in terms of liquid pressure and dissolved hydrogen density as main unknowns, leading mathematically to a nonlinear elliptic–parabolic system of partial differential equations, in which the equations degenerate when the gas phase disappears. We develop a discontinuous Galerkin method for space discretization, combined with a backward Euler scheme for time discretization and an incomplete Newton method for linearization. Numerical examples deal with gas-phase (dis)appearance, ill-prepared initial conditions, and heterogeneous problem with different rock types.     

A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics

In this paper, we present a computational framework for the simulation of coupled flow and reservoir geomechanics. The physical model is restricted to Biot’s theory of single-phase flow and linear poroelasticity, but is sufficiently general to be extended to multiphase flow problems and inelastic behavior. The distinctive technical aspects of our approach are: (1) the space discretization of the equations. The unknown variables are the pressure, the fluid velocity, and the rock displacements. We recognize that these variables are of very different nature, and need to be discretized differently. We propose a mixed finite element space discretization, which is stable, convergent, locally mass conservative, and employs a single computational grid. To ensure stability and robustness, we perform an implicit time integration of the fluid flow equations. (2) The strategies for the solution of the coupled system. We compare different solution strategies, including the fully coupled approach, the usual (conditionally stable) iteratively coupled approach, and a less common unconditionally stable sequential scheme. We show that the latter scheme corresponds to a modified block Jacobi method, which also enjoys improved convergence properties. This computational model has been implemented in an object-oriented reservoir simulator, whose modular design allows for further extensions and enhancements. We show several representative numerical simulations that illustrate the effectiveness of the approach.     

Coupled formulation and algorithms for the simulation of non‐planar three‐dimensional hydraulic fractures using the generalized finite element method

This paper presents a coupled hydro‐mechanical formulation for the simulation of non‐planar three‐dimensional hydraulic fractures. Deformation in the rock is modeled using linear elasticity, and the lubrication theory is adopted for the fluid flow in the fracture. The governing equations of the fluid flow and elasticity and the subsequent discretization are fully coupled. A Generalized/eXtended Finite Element Method (G/XFEM) is adopted for the discretization of the coupled system of equations. A Newton–Raphson method is used to solve the resulting system of nonlinear equations. A discretization strategy for the fluid flow problem on non‐planar three‐dimensional surfaces and a computationally efficient strategy for handling time integration combined with mesh adaptivity are also presented. Several three‐dimensional numerical verification examples are solved. The examples illustrate the generality and accuracy of the proposed coupled formulation and discretization strategies. Copyright © 2015 John Wiley & Sons, Ltd.     

A consistent formulation of a dilatant interface element

The paper considers a plane joint or interface element suitable for implementation into a standard non-linear finite element code. The element is intended to model discontinuities with rough contact surfaces, such as rock joints, where dilatant behaviour is present. Of particular concern is the formulation of a constitutive model which fully caters for all possible histories of opening, closing and sliding (accompained by dilation or contraction) in any direction. The non-linear incremental constitutive equations are formulated in a manner appropriate for a back-ward difference discretization in time along the path of loading. The advantage of such an approach is that no essential distinction need be drawn between opening, closing and sliding. Further, a convenient formulation of the constitutive equations is facilitated by representing the different contact conditions in relative displacement space. The state diagram in relative displacement space, however, changes from one time step to the next, and evolution equations for the updating must be formulated. These concepts are illustrated for two rock-joint models: a sawtooth asperity model and a limited dilation model. The models are based on a penalty formulation to enforce the contact constraints, and explicit equations for the tangent stiffness matrix and for the corrector step of the standard Newton–Raphson iterative algorithm are derived. These equations have been implemented as an user element into the finite element code ABAQUS⁷. Three examples are presented to illustrate the predictions of the formulation.     

Error Estimates for a Time Discretization Method for the Richards' Equation

We present a numerical analysis of a time discretization method applied to Richards' equation. Written in its saturation-based form, this nonlinear parabolic equation models water flow into unsaturated porous media. Depending on the soil parameters, the diffusion coefficient may vanish or explode, leading to degeneracy in the original parabolic equation. The numerical approach is based on an implicit Euler time discretization scheme and includes a regularization step, combined with the Kirchhoff transform. Convergence is shown by obtaining error estimates in terms of the time step and of the regularization parameter.     

Numerical analysis of thermal effects during carbon dioxide injection with enhanced gas recovery: a theoretical case study for the Altmark gas field

Prediction about reservoir temperature change during carbon dioxide injection requires consideration of all, often subtle, thermal effects. In particular, Joule?CThomson cooling (JTC) and the viscous heat dissipation (VHD) effect are factors that cause flowing fluid temperature to differ from the static formation temperature. In this work, warm-back behavior (thermal recovery after injection completed), as well as JTC and VHD effects, at a multi-layered depleted gas reservoir are demonstrated numerically. OpenGeoSys (OGS) is able to solve coupled partial differential equations for pressure, temperature and mole-fraction of each component of the mixture with a combination of monolithic and staggered approaches. The Galerkin finite element approach is adapted for space discretization of governing equations, whereas for temporal discretization, a generalized implicit single-step scheme is used. For numerical modeling of warm-back behavior, we chose a simplified test case of carbon dioxide injection. This test case is numerically solved by using OGS and FeFlow simulators independently. OGS differs from FeFlow in the capability of representing multi-componential effects on warm-back behavior. We verify both code results by showing the close comparison of shut-in temperature profiles along the injection well. As the JTC cooling rate is inversely proportional to the volumetric heat capacity of the solid matrix, the injection layers are cooled faster as compared to the non-injection layers. The shut-in temperature profiles are showing a significant change in reservoir temperature; hence it is important to account for thermal effects in injection monitoring.     

An elliptic problem with integral constraints with application to large-scale geophysical flows

We present a weak formulation of a non-standard elliptic equation whose boundary values are determined in part by integral relations. Existence and uniqueness of its solution are proved, and a finite element discretization is described, analyzed, and implemented on a test problem. The equation is a generalization of one that is solved during integration of the three-dimensional Quasigeostrophic equations, which model large-scale rotating stratified flows, where the integral constraints represent conservation of physical properties.     

A semi-discrete procedure for dynamic response analysis of saturated soils

Biot's equations of wave propagation through fluid-saturated porous elastic media are discretized spatially using the finite element method in conjunction with Galerkin's procedure. Laplace transformation of the discretized equations is used to suppress the time variable. Introducing Laplace transforms of constituent velocities at nodal points as additional variables, the quadratic set of equations in the Laplace transform parameter is reduced to a linear form. The solution in the Laplace transform space is inverted, term by term, to get the complete time history of the solid and fluid displacements and velocities. Since the solution is exact in the time domain, the error in the calculated response is entirely due to the spatial approximation. The procedure is applied to one-dimensional wave propagation in a linear elastic material and in a fluid-saturated elastic soil layer with ‘weak’, ‘strong’ as well as ‘moderate’ coupling. With refinement of the spatial mesh, convergence to the exact solution is established. The procedure can provide a useful benchmark for validation of approximate temporal discretization schemes and estimation of errors due to spatial discretization.     

带横隔板圆柱绕流特性数值模拟

为揭示尾迹区添加横隔板对圆柱绕流流场特性影响,把多步格式引入到特征线算子分裂有限元法中,建立了基于多步格式的特征线算子分裂有限元法:在每个时间步内将Navier-Stokes方程分裂成对流项和扩散项,对流项时间离散采用多步格式,在每一子时间步内沿特征线展开并显式求解。方腔流数值模拟结果表明该算法既可降低对整体时间步长的要求又可提高计算精度。对比有无横隔板圆柱绕流流场和圆柱表面压力变化表明,横隔板可以有效地抑制绕流尾迹区涡旋脱落,提高圆柱背流面压力,减少圆柱上下表面的压力差。     

A stable meshfree method for fully coupled flow-deformation analysis of saturated porous media

A fully coupled meshfree algorithm is proposed for numerical analysis of Biot’s formulation. Spatial discretization of the governing equations is presented using the Radial Point Interpolation Method (RPIM). Temporal discretization is achieved based on a novel three-point approximation technique with a variable time step, which has second order accuracy and avoids spurious ripple effects observed in the conventional two-point Crank Nicolson technique. Application of the model is demonstrated using several numerical examples with analytical or semi-analytical solutions. It is shown that the model proposed is effective in simulating the coupled flow deformation behaviour in fluid saturated porous media with good accuracy and stability irrespective of the magnitude of the time step adopted.     

Finite volume discretization with imposed flux continuity for the general tensor pressure equation

We present a new family of flux continuous, locally conservative, finite volume schemes applicable to the diagonal and full tensor pressure equations with generally discontinuous coefficients. For a uniformly constant symmetric elliptic tensor field, the full tensor discretization is second order accurate with a symmetric positive definite matrix. For a full tensor, an _M-matrix with diagonal dominance can be obtained subject to a sufficient condition for ellipticity. Positive definiteness of the discrete system is illustrated. Convergence rates for discontinuous coefficients are presented and the importance of modeling the full permeability tensor pressure equation is demonstrated.     

Isogeometric analysis for unsaturated flow problems

Unsaturated flow problems in porous media often described by Richards’ equation are of great importance in many engineering applications. In this contribution, we propose a new numerical flow approach based on isogeometric analysis (IGA) for modeling the unsaturated flow problems. The non-uniform rational B-spline (NURBS) basis is utilized for spatial discretization whereas the stable implicit backward Euler method for time discretization. The nonlinear Richards’ equation is iteratively solved with the aid of the Newton–Raphson scheme. Owing to some desirable features of an efficient numerical flow approach, major advantages of the present formulation involve: (a) numerical oscillation at the wetting front can be avoided or facilitated, simply by using either an h-refinement or a lumped mass matrix technique; (b) higher-order exactness can be obtained due to the nature of the IGA features; (c) the approach is straightforward to implement and it does not need any transformation, e.g., Kirchhoff transformation or filter algorithm; and (d) in contrast to the Picard iteration scheme, which forms linear convergences, the proposed approach can however yield quadratic convergences by using the Newton–Raphson method for solving resultant nonlinear equations. Numerical model validation is analyzed by solving a three-dimensional unsaturated flow problem in soil, and its derived results are verified against analytical solutions. Numerical applications are then studied by considering three extensive examples with simple and complex configurations to further show the accuracy and applicability of the present IGA.     

The single-cell transport problem for two-phase flow with polymer

Polymer injection is a widespread strategy in enhanced oil recovery. Polymer increases the water viscosity and creates a more favorable mobility ratio between the injected water and the displaced oil. The computational cost of simulating polymer injection can be significantly reduced if one splits the governing system of two-phase equations into a pressure equation and a set of saturation/component equations and use a Gauss–Seidel algorithm with optimal cell ordering to solve the nonlinear systems arising from an implicit discretization of the saturation/component equations. This approach relies on a robust single-cell solver that computes the saturation and polymer concentration of a cell, given the total flux and the saturation and polymer concentration of the neighboring cells. In this paper, we consider a relatively comprehensive polymer model used in an industry-standard simulator, and show that, in the case of a discretization using a two-point flux approximation, the single-cell problem always admits a solution that is also unique.     

Numerical approximation of water–air two-phase flow by the mixed finite element method

A new model for two-phase flow of water and air in soil is presented. This leads to a system of two mass balance equations and two equations representing conservation of momentum of fluid and gas, respectively. This paper is concerned with the verification of this model for the special case of a rigid soil skeleton by computational experiments. Its numerical treatment is based on the Raviart–Thomas mixed finite element method combined with an implicit Euler time discretization. The feasibility of the method is illustrated for some test examples of one- and two-dimensional two-phase flow problems.     

A time‐discontinuous Galerkin method for the dynamical analysis of porous media

We present a time‐discontinuous Galerkin method (DGT) for the dynamic analysis of fully saturated porous media. The numerical method consists of a finite element discretization in space and time. The discrete basis functions are continuous in space and discontinuous in time. The continuity across the time interval is weakly enforced by a flux function. Two applications and several numerical investigations confirm the quality of the proposed space–time finite element scheme. Copyright © 2006 John Wiley & Sons, Ltd.     

A fully coupled non-linear dynamic analysis procedure and its verification using centrifuge test results

The equations governing the dynamic behavior of saturated porous media as well as a finite element spatial discretization of these equations are summarized. A three-parameter time integration scheme called the Hilber–Hughes–Taylor α-method is used together with a predictor/multi-corrector algorithm, instead of the widely used Newmark's method, to integrate the spatially discrete finite element equations. The new time integration scheme possess quadratic accuracy and desirable numerical damping characteristics. The proposed numerical solution and bounding surface plasticity theory to describe the constitutive behaviour of soil have been implemented as the computer code DYSAC2. Predictions made by DYSAC2 code are verified using dynamic centrifuge test results for a clay embankment. Importance of initial state of a soil on its dynamic behaviour is demonstrated.