A fully implicit and consistent finite element framework for modeling reservoir compaction with large deformation and nonlinear flow model. Part I: theory and formulation

Reservoir depletion results in rock failure, wellbore instability, hydrocarbon production loss, oil sand production, and ground surface subsidence. Specifically, the compaction of carbonate reservoirs with soft rocks often induces large plastic deformation due to rock pore collapse. On the other hand, following the compaction of reservoirs and failure of rock formations, the porosity and permeability of formations will, in general, decrease. These bring a challenge for reservoir simulations because of high nonlinearity of coupled geomechanics and fluid flow fields. In this work, we present a fully implicit, fully coupled, and fully consistent finite element formulation for coupled geomechanics and fluid flow problems with finite deformation and nonlinear flow models. The Pelessone smooth cap plasticity model, an important material model to capture rock compaction behavior and a challenging material model for implicit numerical formulations, is incorporated in the proposed formulation. Furthermore, a stress-dependent permeability model is taken into account in the formulation. A co-rotational framework is adopted for finite deformation, and an implicit material integrator for cap plasticity models is consistently derived. Furthermore, the coupled field equations are consistently linearized including nonlinear flow models. The physical theories, nonlinear material and flow models, and numerical formulations are the focus of part I of this work. In part II, we verify the proposed numerical framework and demonstrate the performance of our numerical formulation using several numerical examples including a field reservoir with soft rocks undergoing serious compaction.     

Extended finite element framework for fault rupture dynamics including bulk plasticity

We present an explicit extended finite element framework for fault rupture dynamics accommodating bulk plasticity near the fault. The technique is more robust than the standard split‐node method because it can accommodate a fault propagating freely through the interior of finite elements. To fully exploit the explicit algorithmic framework, we perform mass lumping on the enriched finite elements that preserve the kinetic energy of the rigid body and enrichment modes. We show that with this technique, the extended FE solution reproduces the standard split‐node solution, but with the added advantage that it can also accommodate randomly propagating faults. We use different elastoplastic constitutive models appropriate for geomaterials, including the Mohr–Coulomb, Drucker–Prager, modified Cam‐Clay, and a conical plasticity model with a compression cap, to capture off‐fault bulk plasticity. More specifically, the cap model adds robustness to the framework because it can accommodate various modes of deformation, including compaction, dilatation, and shearing. Copyright © 2013 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.     

Fully coupled modeling of seabed subsidence and reservoir compaction of North Sea oil fields

This paper focuses on the aspects of fully coupled continuum modeling of multiphase poroelasticity applied to the three-dimensional numerical simulations of the Ekofisk oil reservoir in the North Sea (56°29′–34′N, 03°10′–14′E). A systematic presentation is chosen to present the methodology behind fully coupled, continuum modeling. First, a historical review of the subsidence phenomena above an oil and gas reservoir is given. This will serve as a background against which the relevance of the present approach to compaction and subsidence modeling will be demonstrated. Following this, the governing equations for a multiphase poroelasticity model are briefly presented. Particular attention is paid to the analysis of the pore-compressibility term usually used in an uncoupled approach for characterising the host-rock deformation. A comparative numerical analysis is carried out to contrast and highlight the difference between coupled and uncoupled reservoir simulators. Finally, a finite-element numerical model of the Ekofisk field is presented and a significant result is a contour map of seabed subsidence which is in general agreement with the shape of the subsidence contours based on past bathymetric surveys. Analysis of the simulation reveals that, due to the downward movement of the overburden, oil migration occurs from the crest of the anticline in which the field is situated, towards the flank. The pore-pressure depletion in the reservoir is significantly delayed due to the replenishment of the reservoir energy via the formational compaction. Horizontal movement in the reservoir, which is neglected in traditional modeling, can be significant and comparable in magnitude to the vertical subsidence. Electronic Publication     

Modelling mechanical behaviour of limestone under reservoir conditions

High porosity and low permeability limestone has presented pore collapse. As fluid is withdrawn from these reservoirs, the effective stresses acting on the rock increase. If the strength of the rock is overcome, pore collapse may occur, leading to irreversible compaction of porous media with permeability and porosity reduction. It impacts on fluid withdrawal. Most of reservoirs have been discovered in weak formations, which are susceptible to this phenomenon. This work presents a study on the mechanical behaviour of a porous limestone from a reservoir located in Campos Basin, offshore Brazil. An experimental program was undergone in order to define its elastic plastic behaviour. The tests reproduced the loading path conditions expected in a reservoir under production. Parameters of the cap model were fitted to these tests and numerical simulations were run. The numerical simulations presented a good agreement with the experimental tests. Copyright © 2006 John Wiley & Sons, Ltd.     

A constitutive model for the dynamic and high‐pressure behaviour of a propellant‐like material: Part II: Model development and applications

This paper presents the second part of a work that aims at developing a mechanical model for the behaviour of propellant‐like materials under high confining pressure and strain rate. The model is to be employed to determine the temperature rise due to mechanical dissipation during a dynamic penetration event. Using the kinematical and thermodynamic background derived in the first part, a viscoelastic–viscoplastic–compaction model is put forward and identified. Viscoelasticity and compaction refer to the elastomeric nature of the material, while viscoplasticity, including implicit damage through dilatancy, reflects its granular nature. Some numerical exercises are performed, in view of determining the major model forces and weaknesses, and of assessing the numerical algorithm robustness. A penetration event is preliminarily simulated, and a temperature rise field predicted. Necessary model improvements are finally discussed, together with several ways of research for a longer term approach. Copyright © 2001 John Wiley & Sons, Ltd.     

Implicit and explicit integration schemes in the anisotropic bounding surface plasticity model for cyclic behaviours of saturated clay

Two integration algorithms, namely the implicit return mapping and explicit sub-stepping schemes, are adopted in the anisotropic bounding surface plasticity model for cyclic behaviours of saturated clay and are implemented into finite element code. The model is a representative of a series of bounding surface models that have typical characteristics, including isotropic and kinematic hardening rules and a rotational bounding surface to capture complex but important cyclic behaviours of soils, such as cyclic shakedown and degradation. However, there is no explicit current yield surface in the model to which the conventional implicit algorithm returns the stress state back or the sub-stepping integration corrects the drift of the stress state. Hence, necessary modifications have been made for both of the integration schemes. First, the image stress point is mapped or corrected to the bounding surface instead of mapping back or correcting the stress state to the yield surface. Second, the unloading–loading criterion is checked to determine the image stress point rather than checking the yield criterion after giving the trial stress state in a conventional way. Comparative studies on the accuracy, stability and efficiency of the two integration schemes are conducted not only at the element level but also in solving boundary value problems of monotonic and cyclic bearing behaviours of rigid footings on saturated clay. For smaller strain increments, there is no significant difference in the accuracy between the two integration schemes, but the explicit integration shows a higher efficiency and accuracy. For relatively larger increments, the implicit return mapping algorithm presents good accuracy and more robustness, while the sub-stepping algorithm shows deteriorating accuracy and suffers the convergence problem. With the tolerance used in the present model, the bearing capacity of the rigid footing predicted by the return mapping algorithm is closer to the available analytical and numerical solutions, while the bearing capacity predicted by the sub-stepping algorithm shows a marginal increase.     

Implicit integration algorithm for Hoek-Brown elastic-plastic model

A realistic strength criterion often used to describe the yielding behaviour of a jointed rock mass at a continuum level is the well-known Hoek and Brown criterion. This paper is concerned with a 3-D stress generalization of the Hoek-Brown failure criterion by means of an elliptical functional which leads to a smooth deviatoric trace in the stress space. For its incorporation into a finite element analysis involving plasticity calculations, the formulation of an implicit stress integration algorithm is presented. The key computational methodology alludes to the notion of consistent tangent modulus and implicit return mapping schemes (radial and closest point return) for stress integration in a finite element analysis. Within the context of non-linear elastoplastic analysis, it is found that formulation of such consistent modulus and success into achieving numerical efficiency are closely intertwined. Indeed, the procedure results into accurate and rapid convergence of the displacement finite element scheme during the search for equilibrium. This means that considerable savings in computational time can be achieved for large scale problems. Numerical examples which focus on the Hoek-Brown plasticity model are presented in order to fully appreciate the robustness of the algorithm, and hence the viability of such method to practical problems.     

Numerical solution and convergence analysis of steam injection in heavy oil reservoirs

In this paper, the numerical methods for solving the problem of steam injection in the heavy oil reservoirs are presented. We consider a 3-dimensional model of 3-phase flow, oil, water, and steam, with the effect of 3-phase relative permeability. Interphase mass transfer of water and steam is considered; oil is assumed nonvolatile. We apply the simultaneous solution approach to solve the corresponding nonlinear discretized partial differential equation in the fully implicit form. The convergence of finite difference scheme is proved by the Rosinger theorem. The heuristic Jacobian-Free-Newton-Krylov (HJFNK) method is proposed for solving the system of algebraic equations. The result of this proposed numerical method is well compared with some experimental results. Our numerical results show that the first iteration of the full approximation scheme (FAS) provides a good initial guess for the Newton method. Therefore, we propose a new hybrid-FAS-HJFNK method while there is no steam in the reservoir. The numerical results show that the hybrid-FAS-HJFNK method converges faster than the HJFNK method.     

The influence of anisotropy on compaction bands: The case of coaxiality between stress and fabric anisotropy tensors

Compaction bands are localized failure patterns that appear in highly porous rock material under the effect of relatively high confining pressure. Being affected mainly by volumetric compression, these bands appear to be almost perpendicular to the most compressive principal stress of a stress state at the so-called “cap” of the yield surface (YS). In this study, we focus on the mechanism that leads to the onset of compaction bands by using a viscoplasticity model able to describe the post-localization response of these materials. The proposed constitutive framework is based on the overstress theory of Perzyna (1966) and the anisotropic clay plasticity model of Dafalias (1986), which provides not only the necessary “cap” of the YS, but introduces a rotational hardening (RH) mechanism, thus, accounting for the effect of fabric anisotropy. Following the analysis of Veveakis and Regenauer-Lieb (2015), we identify the compaction bands as “static” cnoidal wave formations in the medium that occur at a post-yield regime, and we study the effect of rotational and isotropic hardening on their onset. Moreover, we determine a theoretical range of confining pressures in triaxial compression tests for the compaction bands to develop. Under the assumption of coaxiality between stress and anisotropy tensors, the results show that the isotropic hardening promotes compaction localization, whereas the RH has a slightly negative effect on the onset of compaction localization.     

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.     

岩性圈闭盖层封闭动力学模型及其应用研究

岩性圈闭的泥页岩盖层阻止油气突破的过程是一个阻力克服动力的过程,本文在此基础上建立了岩性圈闭盖层封闭动力学模型及相应的数学模型。该模型认为岩性油气藏中的油气突破泥岩盖层的过程也是一个力学平衡现象,动力包括储集层地层压力、浮力,阻力主要包括盖层地层压力、毛管压力。通过对单井古孔隙度及古埋深的恢复,可以对比成藏期目的层油气突破泥岩盖层的动力与阻力之间的大小关系,定量研究泥岩盖层封闭性演化过程。应用该模型发现,白家海凸起彩43井侏罗系三工河组及西山窑组在侏罗纪末期缺乏有效泥岩盖层,白垩纪末期三工河组中下部存在有效盖层,而西山窑组盖层的封闭能力一般。     

Peridynamics simulations of geomaterial fragmentation by impulse loads

In this work, we have developed nonlinear peridynamics models of drained and saturated geomaterials, and applied them to simulations of dynamic fragmentation and ejecta formation due to impulse loads. First, we have re‐phrased and re‐interpreted the non‐local state‐based peridynamics formulation to connect the non‐local integral operator with the local differential operator. Second, we have implemented the Drucker–Prager plasticity model in state‐based peridynamics at finite strain. A peridynamics version of the Hughes and Winget algorithm is derived for the constitutive update. Third, we have developed a peridynamics U‐p formulation for saturated geomaterials. Fourth, numerical simulations have been carried out to verify the proposed peridynamics formulations in the simulation of geomaterial fragmentation induced by impulse loads. Comparisons of analytical and numerical results indicate that the peridynamics model has the ability to both match traditional continuum compression examples, as well as simulate complex geomaterial fragmentation processes resulting from impulse loads. Copyright © 2015 John Wiley & Sons, Ltd.     

Optimal experimental design for reservoir property estimates in geothermal exploration

During geothermal reservoir development, drilling deep boreholes turns out to be extremely expensive and risky. Thus, it is of great importance to work out the details of suitable borehole locations in advance. Here, given a set of existing boreholes, we demonstrate how a sophisticated numerical technique called optimal experimental design helps to find a location of an additional exploratory borehole that reduces risk and, ultimately, saves cost. More precisely, the approach minimizes the uncertainty when deducing the effective permeability of a buried reservoir layer from a temperature profile measured in this exploratory borehole. In this paper, we (1) outline the mathematical formulation in terms of an optimization problem, (2) describe the numerical implementation involving various software components, and (3) apply the method to a 3D numerical simulation model representing a real geothermal reservoir in northern Italy. Our results show that optimal experimental design is conceptually and computationally feasible for industrial-scale applications. For the particular reservoir and the estimation of permeability from temperature, the optimal location of the additional borehole coincides with regions of high flow rates and large deviations from the mean temperature of the reservoir layer in question. Finally, the presentation shows that, methodologically, the optimization method can be generalized from estimating permeability to finding any other reservoir properties.     

一种适用于饱和黏土循环动力分析边界面塑性模型的隐式积分算法

修正了传统隐式回映算法,建立了适用于饱和黏土循环动力分析的边界面塑性模型的完全隐式积分格式。该模型基于无弹性域概念和临界状态理论,采用各向同性、运动硬化准则、旋转的边界面,并引入表征土体结构损伤和重塑程度的损伤变量以反映循环载荷作用下饱和黏土的各向异性、刚度、强度软化及塑性变形累积等特征。针对等压固结 和偏压固结 的饱和高岭黏土的不排水三轴试验进行模拟,采用不同的应变增量步长进行计算,并与试验数据对比,结果表明,修正隐式回映算法应用于该类边界面模型的合理性、积分格式的精确性和稳定性;另外,结合有限元软件自动时间步长的增量迭代解法,对饱和黏土应力控制的不排水动三轴试验进行预测,结果表明,修正的适用于该边界面的塑性模型隐式回映算法可以得到比较合理的数值分析结果,能够反映饱和黏土的循环刚度的退化和强度的弱化等动力特性。     

3D geomechanical modeling and estimating the compaction and subsidence of Fahlian reservoir formation (X-field in SW of Iran)

A geomechanical model can reveal the mechanical behavior of rocks and be used to manage the reservoir programs in a better mode. Fluid pressure will be reduced during hydrocarbon production from a reservoir. This reduction of pressure will increase the effective stress due to overburden sediments and will cause porous media compaction and surface subsidence. In some oil fields, the compacting reservoir can support oil and gas production. However, the phenomena can also cause the loss of wells and reduced production and also cause irreparable damage to the surface structures and affect the surrounding environment. For a detailed study of the geomechanical behavior of a hydrocarbon field, a 3D numerical model to describe the reservoir geomechanical characteristics is essential. During this study, using available data and information, a coupled fluid flow-geomechanic model of Fahlian reservoir formation in X-field in SW of Iran was constructed to estimate the amount of land subsidence. According to the prepared model, in this field, the maximum amount of the vertical stress is 110 MPa and the maximum amount of the horizontal stress is 94 MPa. At last, this model is used for the prediction of reservoir compaction and subsidence of the surface. The maximum value of estimated ground subsidence in the study equals to 29 mm. It is considered that according to the obtained values of horizontal and vertical movement in the wall of different wells, those movements are not problematic for casing and well production and also the surrounding environment.     

Numerical modeling to predict the spread of landslide of schist’s area under climatic event of Ain El Hammam (Algeria)

The several reactivations of the landslide of Ain El Hammam (AEH) after each important weather event compel us to look closely at its triggering factors and predict its mechanisms and its evolution at the longer term. In this sense, the prediction of the slope behavior becomes necessary. This paper presents a numerical model of the AEH landslide using Plaxis® software. This model considers hydraulic effects such as precipitation and pore pressure even in the unsaturated parts. Soil and rock behaviors are described with proper elasto-plastic models named Hardening Soils and Jointed Rocks. The first model takes into account hardening on isotropic and deviatoric mechanism as well as a non-associated flow rule. The second model considers a non-isotropic elasticity with perfect plasticity along with given sliding directions. The hydraulic and mechanical models are coupled with an effective stress concept. To detect unstable areas in the landslide, we developed a Matlab® program to take into account the Hill’s bifurcation criterion, which is based on sign of the second-order work. It has been proved that this criterion allows detecting all failure modes that can appear in rate-independent materials and especially the ones that develop before the plasticity limit criterion. From such computations, we can predict the shape and position of slip surface responsible of the actual ground movement of the slope. To validate the numerical results, analysis of field measurement is included. We use high resolution of electrical tomography to delineate the geometry and position of failure surface and approve our results.     

Implicit numerical integration for a kinematic hardening soil plasticity model

Soil models based on kinematic hardening together with elements of bounding surface plasticity, provide a means of introducing some memory of recent history and stiffness variation in the predicted response of soils. Such models provide an improvement on simple elasto‐plastic models in describing soil behaviour under non‐monotonic loading. Routine use of such models requires robust numerical integration schemes. Explicit integration of highly non‐linear models requires extremely small steps in order to guarantee convergence. Here, a fully implicit scheme is presented for a simple kinematic hardening extension of the Cam clay soil model. The algorithm is based on the operator split methodology and the implicit Euler backward integration scheme is proposed to integrate the rate form of the constitutive relations. This algorithm maintains a quadratic rate of asymptotic convergence when used with a Newton–Raphson iterative procedure. Various strain‐driven axisymmetric triaxial paths are simulated in order to demonstrate the efficiency and good performance of the proposed algorithm. Copyright © 2001 John Wiley & Sons, Ltd.     

用流固耦合方法研究油藏压裂后应力应变和孔渗特性变化

油藏压裂后将引起地应力场发生变化，使岩石变形，导致孔隙度和渗透率变化，进而影响产量，为研究这一问题，作者建立了油藏压裂后流－固耦合渗流模型，考虑了以下因素：油藏岩石变形，地应力，孔隙度和渗透率变化，人工裂缝，流体渗流与岩石应变耦合，储藏渗流与裂缝渗流耦合，非达西效应等。较详细地给出了耦合方程及推导过程，控制方程包括的未知变量有压力，饱和度及位移，11个变量，和11个方程，用有限差分方法将流体渗流和岩石应变方程离散成主对角占优的七对角矩阵，可在修改已有三维二相渗流和三维固体力学程序的基础上，采用隐式迭代方法求解，示例分析表明，用此模型可以研究储层应力变变，孔隙度和渗透率随时间和空间变化规律，为开发方案制定，整体压裂设计，压后生产管理等方面提供定量分析技术。