Lijing Zhng, Hu Zhng, Yngung Yun, Shunde Yin

a School of Mathematics and Physics, University of Science and Technology Beijing, Beijing, China

b Department of Petroleum Engineering, University of Wyoming, Laramie, USA

c BitCan Geosciences and Engineering Inc., Calgary, Alberta, T2A 2L5, Canada

d Department of Civil and Environmental Engineering, University of Waterloo, Waterloo, Canada

Keywords: Poroelastoplasticity Reservoir simulation Finite element methods Tangent stiffness method

ABSTRACT As oil and gas extraction activities move into deeper rock formations, many experimental studies and field investigations indicate rock exhibits a plastic behavior rather than a pure linear elastic behavior, so poroelastoplasticity must be taken into account in the reservoir simulation.Because reservoir rock is a porous material consisting of a compressible solid matrix and number of compressible fluids occupying the pore space, fully coupled modeling is required for reservoir simulation considering solid-fluid interaction, complex stress conditions and nonlinear behaviors.But the computational process could be cumbersome when constant tangent stiffness method is used to address the poroelastoplastic behavior.In this paper, a fully coupled poroelastoplasticity reservoir model based on Drucker-Prager yield criterion is implemented the tangent stiffness method, and the computational efficiency is compared with the constant stiffness method.The accuracy of these two methods is demonstrated in one-dimensional consolidation.In a case study, these two methods are used to analyze the stresses and pore pressure of a reservoir and computing results and running efficiency are compared.Also, the linear elastic and nonlinear solutions are compared in one-dimensional consolidation and reservoir modeling.It shows that the difference between results by constant stiffness method and tangent stiffness method is very small, while the tangent stiffness method shows significantly fewer iteration numbers and shorter running time than the constant stiffness method.

1.Introduction

The theory of poroelasticity, originally developed by Terzaghi in one-dimensional consolidation, is the first theory unfolding coupled fluid-soil interactions [1,2].Biot extended this theory to 3-dimensional case [3].Geertsma originally applied poroelasticity to address coupled geomechanical processes during petroleum production operations [4].Later, using the mixture theory, Bowen derived a similar theory for the incompressible and compressible porous media models, which was proven to be consistent with Biot's consolidation theory that has been well accepted [5,6].Subsequently, several studies on poroelasticity have been presented in the literature [7-15].

As oil and gas extraction activities move into deeper rock formations, many experimental studies and field investigations indicate rock exhibits a plastic behavior rather than a pure linear elastic behavior.This requires poroelastoplasticity to be considered in the reservoir simulation.Biot's original theory that assumes linear behavior for the solid matrix has been generalized to complex models dealing with nonlinear problems [16,17].And because of the complexity of the nonlinear constitutive relationship involved, full elastoplastic solutions were obtained by means of numerical methods [18-21].Borja investigated the behavior of linear multistep methods to deal with nonlinear consolidation [22].By applying disturbed state concept (DSC) method, Desai evaluated the solid contact stress and fluid pressure affected by the change of void ratio during deformation for porous saturated materials [23,24].For the steam assisted gravity drainage (SAGD) process, Li and Chalaturnyk utilized a coupled method to conduct reservoir geomechanical simulations [25,26].Ouria et al.analyzed the consolidation of plastic clays under cyclic loading [27].Wong et al.developed a methodology to construct the tensile and shear failure envelopes for compacted unsaturated clayey soils [28].And they also proposed a methodology to decouple the creep-deformation component from the total deformation [29].Some papers considered the change of compressibility and permeability [30-36].Lewis [37], Choo et al.[38], Wan and Eghbalian [39] used double porosity model to analyze the fractured porous media.

Fig.1.Constant stiffness method.

Fig.2.Tangent stiffness method.

Finite element methods are widely used in fully coupled reservoir modeling, but the computational process could be cumbersome when addressing the poroelastoplastic behavior.Computationally efficient solution is the key to analyzing elastoplastic and poroelastoplastic problems.In practical finite element analysis, there are two main types of solution approaches for plasticity simulation [40,41].The first one is the constant stiffness method, in which nonlinearity is considered by iteratively modifying the loads vector, and the global stiffness (usually elastic) matrix is formed once only.Iteration is terminated as long as stresses satisfy certain yield criterion within prescribed tolerances.The second approach is the tangent stiffness method in which the reduction in stiffness of the material is considered.In computer methods for nonlinear poromechanics, it is expected that the tangent stiffness method is more efficient as the number of iterations can be reduced, however, the speculation has not been corroborated, especially when applied in reservoir engineering.

In this paper, fully coupled finite element poroelastoplastic reservoir modeling is implemented by the tangent stiffness method, and its performance is compared with the model implemented by the constant stiffness method.In the first part of the following sections, the accuracy of constant stiffness method and tangent stiffness method is analyzed in a one-dimensional consolidation problem.And in the second part, these two different models are used to analyze the stresses and pore pressure of reservoir, and results and efficiency are compared.

Fig.3.Finite element mesh for 1-D consolidation.

Here we hypothesize that the results between constant stiffness method and tangent stiffness method are similar, and the iteration numbers and running time for the tangent stiffness method is superior to that for the constant stiffness method.

2.Model structure and methodology

2.1.Drucker-Prager yield criterion

In order to get a smooth yield surface approximate to Mohr- Coulomb surface, Drucker and Prager modified Mises yield criterion and put forward the following yield criterion [42].

whereI1is the first stress invariant,J2is the second deviatoric stress invariant, and the material constantskwhich are related to the angle of internal frictionand cohesion of the materialc, are shown in Equations (2) and (3).

Two of the most common approximations used are obtained by making the yield surfaces of the Drucker-Prager and Mohr-Coulomb criteria coincident either at the outer or inner edges of the Mohr- Coulomb surface.Coincidence at the outer edges is obtained when

and coincidence at the inner edges is given by the choice

wherecis the cohesion of the material, andis the internal friction angle.

Table 1 Geometric and mechanical parameters for 1-D consolidation analysis.

Table 2 Comparison of pore pressures between analytical and numerical solution.

2.2.Finite element implementation

The general theory of poroelasticity is based on Biot's theory.With a compressible fluid flowing through a saturated porous medium, the governing equations for the problem of fluid flow in deforming rock can be described as:

Fig.4.Comparison of pore pressures between analytical and numerical solution.

Table 3 Comparison of displacements (unit: mm).

whereGandare Lamé constants.kis the permeability of the porous medium,μis the viscosity of the fluid,uandpdenote the displacement of the porous medium and the pore pressure respectively, the subscripttdenotes time derivative,is the porosity of the porous medium,K,KfandKmare the bulk modulus of the skeleton, fluid and matrix, respectively.Furthermore,iT= [1, 1, 1, 0, 0, 0], andDis the elastic stiffness matrix.

The Galerkin finite element method is used herein to approximate above governing equations.The final form of the FEM solution to the poroelastic equations is as follows [15].

whereM,H,SandCare the elastic stiffness, the flow stiffness, the flow capacity and coupling matrixes, respectively.are the vectors of unknown variablesu,pand corresponding time derivatives.is the vector for the nodal loads and flow source.The explicit expressions of above matrixes are as follows:

To integrate the above equations with respect to time (method), the equation becomes:

Fig.5.Displacement for (a) Z = 0.5H; and (b) Z = H.

Table 4 Comparison of pore pressures (unit: MPa).

If the stresses are in the plastic state at a special timetabased on the yield criterion, equation (11) become equation (12) at the special timeta.

where,is the set of element nodal forces round the node [10],Qis the plastic potential function.

In the process of iteration, the convergence condition can be described as:

wherenormis the norm of vectoris unbalanced force, V is convergence precision.

For different iterative methods, the global stiffness matrices are different, which are shown in Fig.1 and Fig.2.WhereK0is stiffness matrix for constant stiffness method,is stiffness matrix for tangent stiffness method, and K=For the constant stiffness method,(equations (7)-(10)).For the tangent stiffness method, the explicit expressions ofMa,SaandCaare as follows.

Fig.6.Pore pressure for (a) Z = 0.1H; and (b) Z = 0.7H.

Fig.7.Mesh of finite element model for reservoir.

Fig.8.Schematic diagram of reservoir.

whereDepis elastic-plastic matrix.

3.Numerical experiments

3.1.One-dimensional consolidation analysis

Fig.3 shows the finite element mesh used for 1-D consolidation analysis; the number of elements is 10, and each element has a side length of 1 m.In the implementation of current FEM modeling, 20-node isoparametric brick elements are employed for solid field, and 8-node isoparametric brick elements are employed for liquid field.The geometric and mechanical parameters and their values are provided in Table 1.

Table 5 Geometric and mechanical parameters for a reservoir analysis.

Table 6 Horizontal displacements of nodes A and B.

Table 7 Horizontal displacements of nodes C and D.

For linear elastic material 1-D consolidation, the analytical solution of pore pressure is given below [43].

Fig.9.Horizontal displacement of (a) node A; (b) node B; (c) node C; and (d) node D.

Table 8 Surface subsidence for injection well and production well.

Table 2 and Fig.4 show a comparison for pore pressure between the linear elastic analytical and numerical solutions.

For 1-D nonlinear consolidation, Drucker-Prager yield criterion is adopted.The geometric and mechanical parameters and their values are given in Table 1.The comparisons between linear elastic and nonlinear behaviors are shown as follows, and the nonlinear behaviors are based on the models implemented by the constant stiffness and tangent stiffness methods.

Fig.10.Surface subsidence for (a) injection well; and (b) production well.

Table 9 Pore pressures of nodes A and B.

Table 10 Pore pressures of nodes C and D.

Table 3 and Fig.5 show a comparison for vertical displacement between the linear elastic and two different nonlinear solutions at a typical height.It can be seen that the nonlinear methods lead to greater displacements than elastic method, and for constant stiffness and tangent stiffness method, the results are almost the same, which indicates the accuracy of both the constant stiffness and tangent stiffness methods.

Fig.11.Pore pressures of nodes A, B, C and D.

Table 4 and Fig.6 show a comparison for pore pressures between the linear elastic and two different nonlinear solutions at a typical height.It can be seen that the plastic action slows the dissipation of pore water pressure; near final consolidation the elastic and plastic results are approximately equal.At any instant, the predicted pore pressures from the plastic analysis are higher than the predictions from the elastic analysis, giving rise to lower effective stresses.Although the effective stresses in the plastic analysis are lower, the deformations are greater than the elastic predictions (Table 3 and Fig.5).And the plastic flow causes larger deformations at smaller effective stresses, indicating that the net effect of the nonlinearity causes the medium to be softer and hence larger settlements.For constant stiffness and tangent stiffness method, the results are almost the same, which again indicates the accuracy of both constant stiff-ness and tangent stiffness methods.

From the comparison of the results of the poroelastic consolidation model and two poroelastoplastic consolidation models, the conclusion can be made as follows (Tables 3 and 4 and Figs.5 and 6).

There was once a fine gentleman whose entire worldly possessions consisted of a boot-jack and a hair-brush; but he had the most beautiful shirt-collar in the world, and it is about this that we are going to hear a story

(1) The nonlinear models lead to greater final displacements than elastic method.

(2) The plastic deformation slows the dissipation of pore water pressure, but finally the pore pressures of the elastic and plastic methods are approximately equal.

(3) In the middle moment, the predicted pore pressures from the plastic analysis are higher than the predictions from the elastic analysis, giving rise to lower effective stresses.

(4) Both constant stiffness and tangent stiffness method leads to same displacements and pore pressures which indicates the accuracy of both methods.

Table 11 Horizontal and vertical stresses of element (a).

Table 12 Horizontal and vertical stresses of element (b).

3.2.Poroelastoplastic analysis for a reservoir

3D mesh of a reservoir and the surroundings for this model is shown in Figs.7 and 8.In this model, the overburden, sideburden and underburden are obtained by extending 440 m, 700 m and 700 m from top, side and bottom of the reservoir, respectively.The reservoir dimension is 600 m × 600 m × 440 m.The model is fixed in the horizontal direction at each side.The bottom part of boundary is also fixed in the vertical direction.And the top part of the boundary is free.The surrounding rock of reservoir is considered to be elastic deformation, whereas reservoir rock is assumed to be elastoplastic deformation.In the implementation of current FEM modeling, 20- node isoparametric brick elements are employed for solid field, and 8- node isoparametric brick elements are employed for liquid field.The geometric and mechanical parameters and their values are provided in Table 5.

Tables 6 and 7 and Fig.9 show comparison of horizontal displacements of nodes A, B, C and D for linear elastic and nonlinear behavior.The results for the nonlinear plastic behavior are based on the reservoir models implemented by the initial (constant) stiffness and tangent stiffness methods, respectively.It can be seen that the nonlinear modeling leads to greater displacement than elastic modeling.It can also be seen that the difference between the displacements (at nodes A, B, C and D) obtained by initial stiffness and tangent stiffness methods is less than 5%.

Table 8 and Fig.10 show comparison of surface subsidence for injection well and production well for linear elastic and two nonlinear modeling.It can be seen that nonlinear modeling leads to greater surface subsidence than elastic modeling.For two nonlinear methods, which are initial stiffness and tangent stiffness methods, the difference of surface subsidence for injection well and production well is less than 5%.

Tables 9 and 10 and Fig.11 show comparison of pore pressure of nodes A, B, C and D for linear elastic and two nonlinear modeling including initial stiffness and tangent stiffness approaches.It can be seen that the pore pressures of nodes A, B, C and D among linear elastic and two nonlinear methods are almost the same, and plasticity of a few local elements cannot change the dissipation velocity of pore water.

Tables 11-13 and Fig.12 show the comparison of effective stresses in elements (a), (b) and (c) for poroelastic and poroelastoplastic behavior, and the latter is based on the reservoir models implemented by both initial stiffness and tangent stiffness methods.It can be seen that effective stresses of the poroelastoplastic reservoir model are greater than those of the poroelastic reservoir model.The results also show that the difference between effective stresses of elements (a), (b) and (c) obtained by the initial stiffness and tangent stiffness methods is less than 2%.

Comparison of the number of iterations and running time incurred by the initial and tangent stiffness methods for each timestep in the transient reservoir modeling is shown in Table 14.It can be seen that the difference of average iteration time per iteration is not much between initial stiffness method and tangent stiffness method, but the total number of iterations and the total running time incurred by the initial stiffness method are much larger than those of the tangent stiffness method.

Fig.12.(a) Horizontal stress of elements (a); (b) vertical stress of elements (a); (c) Horizontal stress of elements (b); (d) vertical stress of elements (b); (e) Horizontal stress of elements (c); and (f) vertical stress of elements (c).

In order to further compare the computational efficiency of the initial and tangent stiffness methods, the number of iterations and running time incurred by each method are shown in Table 15 and Figs.14 and 15, for different scenarios in terms of number of elements and convergence precision in poroelastoplastic reservoir modeling.It can be seen that as the number of elements and convergence precision increase, the number of iterations and running time are significantly increased for the initial stiffness method while are moderately increased for the tangent stiffness method.So, for the poroelastoplastic reservoir modeling that requires higher convergence precision and larger number of elements, the tangent stiffness method performs much better than the initial stiffness method.

4.Conclusions

In this paper, a fully coupled reservoir model is used to deal with an ideal poroelastoplastic problem based on Drucker-Prager yield criterion by both constant and tangent stiffness methods.

In solving a 1-D consolidation problem for verification purpose, the results by the constant and tangent stiffness methods are almost the same, which indicates the accuracy of both the constant stiffness and tangent stiffness methods.

Fig.13.Schematic diagram of reservoir for plastic region.

Table 13 Horizontal and vertical stresses of element (c).

Table 14 Iterations and running time for initial and tangent stiffness methods.

Table 15 Iterations and running time for different element number and V value.

Fig.14.Comparison of iteration number for initial and tangent stiffness methods in poroelastoplastic reservoir modeling.

In a case study of reservoir simulation, the difference between the results by the constant and tangent stiffness methods is very small.For initial stiffness method, each iteration time is less than tangent stiffness method, but the number of iterations and running time are much larger than tangent stiffness method.As the number of elements and the convergence precision increase, the number of iterations and running time are significantly increased.Therefore, in the case of a higher convergence precision and a larger number of elements required for a model, the tangent stiffness method is significantly more efficient than the constant or initial stiffness method.

Fig.15.Running time: (a) number of elements = 294; (b) number of elements = 384; (c) number of elements = 486; and (d) number of elements = 600.

References

[1] K.V.Terzaghi, Die.Berecbnung der Durchlassigkeitsziffer des Tones aus dem Verlauf der hydrodynamischen Spannungserscheinungen.Akademie der Wissenschaften in Wien, Sitzungsberichte Mathematisch-naturwissenschaftliche Klasse 132 (1923) 125-138 Part 11a.

[2] K.Terzaghi, Erdbaumechanik auf bodenphysikalischer Grundlage, Deuticke.Wien, 1925.

[3] M.A.Biot, General theory of three-dimensional consolidation, J.Appl.Phys.12 (1941) 155-164.

[4] J.Geertsma, A remark on the analogy between thermoelasticity and the elasticity of saturated porous media, J.Mech.Phys.Solid.6 (1957) 13-16.

[5] R.M.Bowen, Incompressible porous media models by use of the theory of mixtures, Int.J.Eng.Sci.18 (9) (1980) 1129-1148.

[6] R.M.Bowen, Compressible porous media models by use of the theory of mixtures, Int.J.Eng.Sci.20 (6) (1982) 697-735.

[7] J.R.Rice, M.P.Cleary, Some basic stress diffusion solutions for fluid-saturated elastic porous media with compressible constituents, Rev.Geophys.Space Phys.14 (1976) 227-241.

[8] E.Detournay, A.H.D.Cheng, Poroelastic response of a borehole in a non-hydrostatic stress field, Int.J.Rock Mech.Min.Sci.Geomech.Abstr.25 (3) (1988) 171-182.

[9] M.Gutierrez, H.Hansteen, Fully coupled analysis of reservoir compaction and subsidence, Proc.Eur.Pet.Conf.2 (1994) 339-347.

[10] R.W.Lewis, B.A.Schrefler, The Finite Element Method in the Static and Dynamic Deformation and Consolidation of Porous Media, John Wiley, New York, 1998.

[11] L.Y.Chin, R.Raghavan, L.K.Thomas, Fully coupled geomechanics and fluid-flow analysis of wells with stress dependent permeability, SPE J.5 (1) (2000) 32-45.

[12] L.Y.Chin, L.K.Thomas, J.E.Sylte, R.G.Pierson, Iterative coupled analysis of geomechanics and fluid flow for rock compaction in reservoir simulation, Oil Gas Sci.Technol.57 (5) (2002) 485-497.

[13] G.Gambolati, M.Ferronato, P.Teatini, R.Deidda, G.Lecca, Finite element analysis of land subsidence above depleted reservoirs with pore pressure gradient and total stress formulations, Int.J.Numer.Anal.Methods GeoMech.25 (2001) 307-327.

[14] S.Yin, L.Rothenburg, M.B.Dusseault, 3D coupled displacement discontinuity and finite element analysis of reservoir behavior during production in semi-infinite domain, J.Transp.Porous Media 65 (2006) 425-441.

[15] S.Yin, Geomechanics-Reservoir Modeling by Displacement Discontinuity-Finite Element Method, [D].Waterloo, University of Waterloo, 2008.

[16] J.C.Small, J.R.Booker, E.H.Davis, Elastoplastic consolidation of soil, Int.J.Solid Struct.12 (1976) 431-448.

[17] O.Coussy, A general theory of thermo-poroelastoplasticity for saturated porous materials, Transport Porous Media 4 (1989) 281-293.

[18] O.C.Zienkiewicz, T.Shiomi, Dynamic behaviour of saturated porous media; the generalized Biot formulation and its numerical solution, Int.J.Numer.Anal.Methods GeoMech.8 (1984) 71-96.

[19] D.M.Potts, L.Zdravkovic, Finite Element Analysis in Geotechnical Engineering: Theory, Thomas Telford, London, 1999.

[20] Y.Yang, H.S.Yu, K.K.Muraleetharan, Solution existence conditions for elastoplastic constitutive models of granular materials, Int.J.Plast.21 (12) (2005) 2406-2425.

[21] Y.Yang, K.K.Muraleetharan, H.S.Yu, Generalized trapezoidal numerical integration of an advanced soil model, Int.J.Numer.Anal.Methods GeoMech.32 (1) (2008) 43-64.

[22] R.I.Borja, One-step and linear multistep methods for nonlinear consolidation, Comput.Methods Appl.Mech.Eng.85 (1991) 239-272.

[23] C.S.Desai, Mechanics of Materials and Interfaces: the Disturbed State Concept, CRC Press, Boca Raton, Fla, 2001.

[24] C.S.Desai, Z.Wang, Disturbed state model for porous saturated materials, Int.J.GeoMech.3 (2) (2003) 260-265.

[25] P.Li, R.J.Chalaturnyk, M.Polikar, Issues with reservoir geomechanical simulations of the SAGD process, Paper 2002-130 Presented at the Petroleum Society's Canadian International Petroleum Conference, Canada, Calgary, Alberta, 2002June 11-13.

[26] P.Li, R.J.Chalaturnyk, T.B.Tan, Coupled reservoir geomechanical simulations for the SAGD process, Paper 2003-083 Presented at the Petroleum Society's Canadian International Petroleum Conference, Canada, Calgary, Alberta, 2003June 10-12.

[27] A.Ouria, C.S.Desai, V.Toufigh, Disturbed state concept-based solution for consolidation of plastic clays under cyclic loading, Int.J.GeoMech.15 (1) (2015) 04014039.

[28] C.K.Wong, R.G.Wan, R.C.K.Wong, Tensile and shear failure behaviour of compacted clay-hybrid failure mode, Int.J.Geotech.Eng.14 (3) (2020) 231-241.

[29] C.K.Wong, R.G.Wan, R.C.K.Wong, Methodology for estimating creep deformation from consolidation deformation in 1D compression, Int.J.GeoMech.(6) (2018) 04018042.

[30] P.J.Fox, J.D.Berles, CS2: a piecewise-linear model for large strain consolidation, Int.J.Numer.Anal.Methods GeoMech.21 (7) (1997) 453-475.

[31] P.J.Fox, H.Pu, Enhanced CS2 model for large strain consolidation, Int.J.GeoMech.(5) (2012) 574-583.

[32] R.C.K.Wong, A model for strain-induced permeability anisotropy in deformable granular media, Can.Geotech.J.40 (2003) 95-106.

[33] J.F.Shao, G.Duveau, F.Bourgeois, W.Z.Chen, Elastoplastic damage modeling in unsaturated rocks and applications, Int.J.GeoMech.6 (2) (2006) 119-130.

[34] K.M.Hall, P.J.Fox, Large strain consolidation model for land subsidence, Int.J.GeoMech.18 (11) (2018) 06018028.

[35] J.C.Chai, Y.Zhou, Method for considering the effect of nonuniform consolidation, Int.J.GeoMech.18 (2) (2018) 04017151.

[36] B.Li, Y.G.Fang, Z.F.Ou, Asymptotic solution for the one-dimensional nonlinear consolidation equation including the pore evolution effect, Int.J.GeoMech.18 (10) (2018) 04018125.

[37] R.W.Lewis, H.R.Ghafouri, A novel finite element double porosity model for multiphase flow through deformable fractured porous media, Int.J.Numer.Anal.Methods GeoMech.21 (11) (1997) 789-816.

[38] J.Choo, J.A.White, R.I.Borja, Hydromechanical modeling of unsaturated flow in double porosity media, Int.J.GeoMech.(2016) D4016002.

[39] R.G.Wan, M.Eghbalian, Hydro-mechanical description of fractured porous media based on microporomechanics, Int.J.Solid Struct.96 (2016) 274-287.

[40] H.J.Siriwardane, C.S.Desai, Two numerical schemes for nonlinear consolidation, Int.J.Numer.Methods Eng.17 (1981) 405-426.

[41] N.Manoharan, S.P.Dasgupta, Consolidation analysis of elasto-plastic soil, Comput.Struct.54 (6) (1995) 1005-1021.

[42] D.C.Drucker, W.Prager, Soil mechanics and plastic analysis for limit design, Q.Appl.Math.10 (2) (1952) 157-165.

[43] K.Terzaghi, B.R.Peck, G.Mesri, Soli Mechanics in Engineering Practice, third ed., Wiley Inc., New York, 1999.