Numerical convergence study of iterative coupling for coupled flow and geomechanics

In this paper, we consider algorithms for modeling complex processes in porous media that include fluid and structure interactions. Numerous field applications would benefit from a better understanding and integration of porous flow and solid deformation. Important applications in environmental and petroleum engineering include carbon sequestration, surface subsidence, pore collapse, cavity generation, hydraulic fracturing, thermal fracturing, wellbore collapse, sand production, fault activation, and waste disposal, while similar issues arise in biosciences and chemical sciences as well. Here, we consider solving iteratively the coupling of flow and mechanics. We employ mixed finite element method for flow and a continuous Galerkin method for elasticity. For single-phase flow, we demonstrate the convergence and convergence rates for two widely used schemes, the undrained split and the fixed stress split. We discuss the extension of the fixed stress iterative coupling scheme to an equation of state compositional flow model coupled with elasticity and a single-phase poroelasticity model on general hexahedral grids. Computational results are presented.     

Accelerating the convergence of coupled geomechanical‐reservoir simulations

The pressure variations during the production of petroleum reservoir induce stress changes in and around the reservoir. Such changes of the stress state can induce marked deformation of geological structures for stress sensitive reservoirs as chalk or unconsolidated sand reservoirs. The compaction of those reservoirs during depletion affects the pressure field and so the reservoir productivity. Therefore, the evaluation of the geomechanical effects requires to solve in a coupling way the geomechanical problem and the reservoir multiphase fluid flow problem. In this paper, we formulate the coupled geomechanical‐reservoir problem as a non‐linear fixed point problem and improve the resolution of the coupling problem by comparing in terms of robustness and convergence different algorithms. We study two accelerated algorithms which are much more robust and faster than the conventional staggered algorithm and we conclude that they should be used for the iterative resolution of coupled reservoir‐geomechanical problem. Copyright © 2006 John Wiley & Sons, Ltd.     

A new mixed finite element method for poro‐elasticity

Development of robust numerical solutions for poro‐elasticity is an important and timely issue in modern computational geomechanics. Recently, research in this area has seen a surge in activity, not only because of increased interest in coupled problems relevant to the petroleum industry, but also due to emerging applications of poro‐elasticity for modelling problems in biomedical engineering and materials science. In this paper, an original mixed least‐squares method for solving Biot consolidation problems is developed. The solution is obtained via minimization of a least‐squares functional, based upon the equations of equilibrium, the equations of continuity and weak forms of the constitutive relationships for elasticity and Darcy flow. The formulation involves four separate categories of unknowns: displacements, stresses, fluid pressures and velocities. Each of these unknowns is approximated by linear continuous functions. The mathematical formulation is implemented in an original computer program, written from scratch and using object‐oriented logic. The performance of the method is tested on one‐ and two‐dimensional classical problems in poro‐elasticity. The numerical experiments suggest the same rates of convergence for all four types of variables, when the same interpolation spaces are used. The continuous linear triangles show the same rates of convergence for both compressible and entirely incompressible elastic solids. This mixed formulation results in non‐oscillating fluid pressures over entire domain for different moments of time. The method appears to be naturally stable, without any need of additional stabilization terms with mesh‐dependent parameters. Copyright © 2007 John Wiley & Sons, Ltd.     

Consistent tangent operators for constitutive rate equations

A general approach for obtaining the consistent tangent operator for constitutive rate equations is presented. The rate equations can be solved numerically by the user's favourite time integrator. In order to obtain reliable results, the substepping in integration should be based on a control of the local error. The main ingredient of the consistent tangent operator, namely the derivative of the stress with respect to the strain increment must be computed simultaneously with the same integrator, applied to a numerical approximation of the variational equations. This information enables finite‐element packages to assemble a consistent tangent operator and thus guarantees quadratic convergence of the equilibrium iterations. Several numerical examples with a hypoplastic constitutive law are given. As numerical integrator we used a second‐order extrapolated Euler method. Quadratic convergence of the equilibrium iteration is shown. Copyright © 2002 John Wiley & Sons, Ltd.     

Three‐dimensional modeling of problems in poro‐elasticity via a mixed least‐squares method using linear tetrahedral elements

In a previous publication we developed a new mixed least‐squares method for poro‐elasticity. The approximate solution was obtained via a minimization of a least‐squares functional, based upon the equations of equilibrium, the equations of continuity and weak forms of the constitutive relationships for elasticity and Darcy flow. The formulation involved four independent types of variables: displacements, stresses, pore pressures and velocities. All of them were approximated by linear continuous triangles. Encouraged by the computational results, obtained from the two‐dimensional implementation of the method, we extended our formulation to three dimensions. In this paper we present numerical examples for the performance of continuous linear tetrahedra within the context of the mixed least‐squares method. The initial results suggest that the method works well in the nearly and entirely incompressible limits for elasticity. For poro‐elasticity, the obtained pore pressures are stable without exhibiting the oscillations, which are observed when the standard Galerkin formulation is used. Copyright © 2010 John Wiley & Sons, Ltd.     

Convergence analysis of fixed stress split iterative scheme for anisotropic poroelasticity with tensor Biot parameter

We perform a convergence analysis of the fixed stress split iterative scheme for the Biot system modeling coupled flow and deformation in anisotropic poroelastic media with tensor Biot parameter. The fixed stress split iterative scheme solves the flow subproblem with all components of the stress tensor frozen using a multipoint flux mixed finite element method, followed by the poromechanics subproblem using a conforming Galerkin method in every coupling iteration at each time step. The coupling iterations are repeated until convergence and Backward Euler is employed for time marching. The convergence analysis is based on studying the equations satisfied by the difference of iterates to show that the fixed stress split iterative scheme for anisotropic poroelasticity with Biot tensor is contractive. We also demonstrate that the scheme is numerically convergent using the classical Mandel’s problem solution for transverse isotropy.     

A priori error estimates for the numerical solution of a coupled geomechanics and reservoir flow model with stress-dependent permeability

In this paper we consider the numerical solution of a coupled geomechanics and a stress-sensitive porous media reservoir flow model. We combine mixed finite elements for Darcy flow and Galerkin finite elements for elasticity. This work focuses on deriving convergence results for the numerical solution of this nonlinear partial differential system. We establish convergence with respect to the L ²-norm for the pressure and for the average fluid velocity and with respect to the H ¹-norm for the deformation. Estimates with respect to the L ²-norm for mean stress, which is of special importance since it is used in the computation of permeability for poro-elasticity, can be derived using the estimates in the H ¹-norm for the deformation. We start by deriving error estimates in a continuous-in-time setting. A cut-off operator is introduced in the numerical scheme in order to derive convergence. The spatial grids for the discrete approximations of the pressure and deformation do not need be the same. Theoretical convergence error estimates in a discrete-in-time setting are also derived in the scope of this investigation. A numerical example supports the convergence results.     

A single numerically efficient equation for approximating the Mohr–Coulomb and the Matsuoka–Nakai failure criteria with rounded edges and apex

This paper presents a novel formulation for defining soil failure. It plots in the principal stress space as a surface with the shape ranging between an approximation of the Matsuoka–Nakai and of the Mohr–Coulomb criteria depending on the value of a single parameter. The new function can be used as a replacement of the original equations of these well‐established criteria for implementing in a program for numerical analyses, and it is particularly effective for approximating the Matsuoka–Nakai criterion. Both the Mohr–Coulomb and the Matsuoka–Nakai failure criteria present numerical difficulties during implementation and also at run‐time. In the case of the Matsuoka–Nakai, the new formulation plots in the first octant only, whereas the original criterion plots in all octants, which causes severe convergence problems particularly for those Gauss points with low stress state, such as those on the side of a shallow footing. When the shape parameter is set to reproduce the Mohr–Coulomb failure criterion, on the other hand, the new formulation plots as a pyramid with rounded edges. Moreover, as the new function is at least of class C², the second derivatives are continuous, thus ensuring quadratic convergence of the Newton's method used within the integration scheme of the constitutive law. The proposed formulation can also provide both sharp and rounded apex of the surface at the origin of the stress space by setting accordingly one additional parameter. Copyright © 2013 John Wiley & Sons, Ltd.     

Stress‐driven integration strategies and m‐AGC tangent operator for Perzyna viscoplasticity and viscoplastic relaxation: application to geomechanical interfaces

The paper proposes a stress‐driven integration strategy for Perzyna‐type viscoplastic constitutive models, which leads also to a convenient algorithm for viscoplastic relaxation schemes. A generalized trapezoidal rule for the strain increment, combined with a linearized form of the yield function and flow rules, leads to a plasticity‐like compliance operator that can be explicitly inverted to give an algorithmic tangent stiffness tensor also denoted as the m‐AGC tangent operator. This operator is combined with the stress‐prescribed integration scheme, to obtain a natural error indicator that can be used as a convergence criterion of the intra‐step iterations (in physical viscoplasticity), or to a variable time‐step size in viscoplastic relaxation schemes based on a single linear calculation per time step. The proposed schemes have been implemented for an existing zero‐thickness interface constitutive model. Some numerical application examples are presented to illustrate the advantages of the new schemes proposed. Copyright © 2016 John Wiley & Sons, Ltd.     

u‐p semi‐Lagrangian reproducing kernel formulation for landslide modeling

This paper presents a u‐p (displacement‐pressure) semi‐Lagrangian reproducing kernel (RK) formulation to effectively analyze landslide processes. The semi‐Lagrangian RK approximation is constructed based on Lagrangian discretization points with fixed kernel supports in the current configuration. As a result, it tracks state variables at discretization points while allowing extreme deformation and material separation that is beyond the capability of Lagrangian formulations. The u‐p formulation following Biot theory is incorporated into the formulation to describe poromechanics of saturated geomaterials. In addition, a stabilized nodal integration method to ensure stability of the domain integration and kernel contact algorithms to model contact between bodies are introduced in the u‐p semi‐Lagrangian RK formulation. The proposed method is verified with several numerical examples and validated with an experimental result and the field data of an actual landslide.     

Robust iterative schemes for non-linear poromechanics

We consider a non-linear extension of Biot’s model for poromechanics, wherein both the fluid flow and mechanical deformation are allowed to be non-linear. Specifically, we study the case when the volumetric stress and the fluid density are non-linear functions satisfying certain assumptions. We perform an implicit discretization in time (backward Euler) and propose two iterative schemes for solving the non-linear problems appearing within each time step: a splitting algorithm extending the undrained split and fixed stress methods to non-linear problems, and a monolithic L-scheme. The convergence of both schemes are shown rigorously. Illustrative numerical examples are presented to confirm the applicability of the schemes and validate the theoretical results.     

A high‐resolution local RBF collocation method for steady‐state poroelasticity and hydromechanical damage analysis

In this work, we describe a meshless numerical method based on local collocation with RBFs for the solution of the poroelasticity equation. The RBF finite collocation approach forms a series of overlapping nodal stencils, over which an RBF collocation is performed. These local collocation systems enforce the governing PDE operator throughout their interior, with the intersystem communication occurring via the collocation of field variables at the stencil periphery. The method does not rely on a generalised finite differencing approach, whereby the governing partial differential operator is reconstructed at the global level to drive the solution of the PDE. Instead, the PDE governing and boundary operators are enforced directly within the local RBF collocation systems, and the sparse global assembly is formed by reconstructing the value of the field variables at the centrepoint of the local stencils. In this way, the solution of the PDE is driven entirely by the local RBF collocation, and the method more closely resembles the approach of the full‐domain RBF collocation method. By formulating the problem in this fashion, high rates of convergence may be attained without the computational cost and numerical ill‐conditioning issues that are associated with the full‐domain RBF collocation approach. An analytical solution is formulated for a 2D poroelastic fluid injection scenario and is used to verify the proposed implementation of the method. Highly accurate solutions are produced, and convergence rates in excess of sixth order are observed for each field variable (i.e. pressure and displacement) and field‐variable derivative (i.e. pressure gradients and stresses). The stress and displacement fields resulting from the solution of the poroelasticity equation are then used to describe the formation and propagation of microfractures and microfissures, which may form in the presence of large shear strain, in terms of a continuous damage variable which modifies the mechanical and hydraulic properties of the porous medium. The formation of such hydromechanical damage, and the resulting increase in hydraulic conductivity, is investigated for a pressurised injection into sandstone. Copyright © 2014 John Wiley & Sons, Ltd.     

Evolution of elastic properties in finite poroplasticity and finite element analysis

The formulation of the poroelastoplastic constitutive equations at large strains of a fully saturated material is performed focusing on the usually ignored influence of large strain plasticity on the poroelastic properties. A micromechanics approach allows to take into account the evolution of the microstructure geometry which in turn induces a coupling between elasticity and plasticity. Such a coupling results in an additional term in the macroscopic Cauchy stress rate equation derived from inclusion‐based estimates that leads to a modified Jaumann derivative. The pressure rate equation is also analysed. The finite element discretization of finite poroplasticity is then presented taking into account the elasticity–plasticity coupling. Application to the consolidation situation shows that the coupling may lead to non‐negligible effects. Copyright © 2002 John Wiley & Sons, Ltd.     

A multiscale finite element formulation for axisymmetric elastoplasticity with volumetric locking

Similar to plane strain, axisymmetric stress problem is also highly kinematics constrained. Standard displacement‐based finite element exhibits volumetric locking issue in simulating nearly/fully incompressible material or isochoric plasticity under axisymmetric loading conditions, which severely underestimates the deformation and overestimates the bearing capacity for structural/geotechnical engineering problems. The aim of this paper is to apply variational multiscale method to produce a stabilized mixed displacement–pressure formulation, which can effectively alleviate the volumetric locking issue for axisymmetric stress problem. Both nearly incompressible elasticity and isochoric J₂ elastoplasticity are investigated. First‐order 3‐node triangular and 4‐node quadrilateral elements are tested for locking issues. Several representative simulations are provided to demonstrate the performance of the linear elements, which include the convergence study and comparison with closed‐form solutions. A comparative study with pressure Laplacian stabilized formulation is also presented. Copyright © 2009 John Wiley & Sons, Ltd.     

A model for stress and plastic strain induced nonlinear,hyperelastic anisotropy in soils

The experimental evidence that cohesive and granular soils possess an elastic range in which the elasticity is both nonlinear and anisotropic—with stiffness and directional characteristics strongly dependent on stress and plastic strain (the so‐called ‘stress history’)—is given a formulation based on hyperelasticity. This is accomplished within the framework of elastoplastic coupling, through a new proposal of elastic potentials and a combined use of a plastic‐strain‐dependent fabric tensor and nonlinear elasticity. When used within a simple elastoplastic framework, the proposed model is shown to yield very accurate simulations of the evolution of elastic properties from initial directional stiffening to final isotropic degradation. Within the proposed constitutive framework, it is shown that predictions of shear band formation and evolution become closer to the existing experimental results, when compared to modelling in which elasticity does not evolve. Copyright © 2007 John Wiley & Sons, Ltd.     

Implicit integration of a mixed isotropic–kinematic hardening plasticity model for structured clays

In recent years, a number of constitutive models have been proposed to describe mathematically the mechanical response of natural clays. Some of these models are characterized by complex formulations, often leading to non‐trivial problems in their numerical integration in finite elements codes. The paper describes a fully implicit stress‐point algorithm for the numerical integration of a single‐surface mixed isotropic–kinematic hardening plasticity model for structured clays. The formulation of the model stems from a compromise between its capability of reproducing the larger number of features characterizing the behaviour of structured clays and the possibility of developing a robust integration algorithm for its implementation in a finite elements code. The model is characterized by an ellipsoid‐shaped yield function, inside which a stress‐dependent reversible stiffness is accounted for by a non‐linear hyperelastic formulation. The isotropic part of the hardening law extends the standard Cam‐Clay one to include plastic strain‐driven softening due to bond degradation, while the kinematic hardening part controls the evolution of the position of the yield surface in the stress space. The proposed algorithm allows the consistent linearization of the constitutive equations guaranteeing the quadratic rate of asymptotic convergence in the global‐level Newton–Raphson iterative procedure. The accuracy and the convergence properties of the proposed algorithm are evaluated with reference to the numerical simulations of single element tests and the analysis of a typical geotechnical boundary value problem. Copyright © 2007 John Wiley & Sons, Ltd.     

青藏高原壳幔形变数值模拟研究

现有数值模拟研究已在很大程度上较合理地给出青藏高原演化运动学和动力学过程的图像。利用连续介质快速拉格朗日分析方法,笔者进行了青藏高原壳幔形变数值模拟研究。据此得到的青藏高原三维壳幔形变特征反映纬向上主碰撞带远、近程效应的差异和经向上地壳物质“逃逸”的存在,印证了青藏高原形成过程中南北双向挤压、而且南部作用大于北部作用的可能应力场特征。青藏高原壳幔形变不仅强烈依赖于随深度变化的岩石力学性质及其距离挤压作用前锋带的远近,而且存在强烈的横向不均一性。同时,强应变(剪切)带的存在对高原岩石圈形变具有重要影响,高原形变过程中地壳尺度的耦合流及壳-幔解耦共存。但是,常规数值模拟研究尚存在很大局限性:(1)物理-力学模型单一;(2)几何模型简单;(3)边界形态与条件理想化;(4)模型内部块体划分粗糙;(5)不连续体介质处理困难。借助具有可处理大形变能力的4-D数值模拟方法,将观测资料与数值模拟相互补充是深入研究青藏高原壳幔形变的关键。     

Operator splitting methods for degenerate convection–diffusion equations II: numerical examples with emphasis on reservoir simulation and sedimentation

We present an accurate numerical method for a large class of scalar, strongly degenerate convection–diffusion equations. Important subclasses are hyperbolic conservation laws, porous medium type equations, two-phase reservoir flow equations, and strongly degenerate equations coming from the recent theory of sedimentation–consolidation processes. The method is based on splitting the convective and the diffusive terms. The nonlinear, convective part is solved using front tracking and dimensional splitting, while the nonlinear diffusion part is solved by an implicit–explicit finite difference scheme. In addition, one version of the implemented operator splitting method has a mechanism built in for detecting and correcting unphysical entropy loss, which may occur when the time step is large. This mechanism helps us gain a large time step ability for practical computations. A detailed convergence analysis of the operator splitting method was given in Part I. Here we present numerical experiments with the method for examples modelling secondary oil recovery and sedimentation–consolidation processes. We demonstrate that the splitting method resolves sharp gradients accurately, may use large time steps, has first order convergence, exhibits small grid orientation effects, has small mass balance errors, and is rather efficient.     

Numerical solution for consolidation and desiccation of soft soils

The consolidation and desiccation behaviour of soft soils can be described by two time‐dependent non‐linear partial differential equations using the finite strain theory. Analytical solutions do not exist for these governing equations. In this paper, we develop efficient numerical methods and software for finding the numerical solutions. We introduce a semi‐implicit time integration scheme, and show numerically that our method converges. In addition, the numerical solution matches well with the experimental result. A boundary refinement method is also developed to improve the convergence and stability for the case of Neumann type boundary conditions. Interface governing equations are derived to maintain the continuity of consolidation and desiccation processes. This is useful because the soil column can undergo desiccation on top and consolidation on the bottom simultaneously. The numerical algorithms has been implemented into a computer program and the results have been verified with centrifuge test results conducted in our laboratory. Copyright © 2001 John Wiley & Sons, Ltd.     

An internally consistent integration method for critical state models

A new procedure to integrate critical state models including Cam–Clay and modified Cam–Clay is proposed here. The proposed procedure makes use of the linearity of the virgin isotropic compression curve and the parallel anisotropic consolidation lines in e–ln p space which are basic features of the formulation of critical state models. Using this algorithm, a unique final stress state may be found as a function of a single unknown for elastoplastic loading. The key equations are given in this article for the Cam–Clay and modified Cam–Clay models. The use of the Newton–Raphson iterative method to minimize residuals and obtain a converged solution is described here. This new algorithm may be applied using the assumptions of linear elasticity or non‐linear elasticity within a given loading step. The new algorithm proposed here is internally consistent and has computational advantages over the current numerical integration procedures. Numerical examples are presented to show the performance of the algorithm as compared to other integration algorithms. Published in 2005 by John Wiley & Sons, Ltd.