Non-local continuum damage model for poro-viscoelastic porous media

https://doi.org/10.1016/j.ijrmms.2022.105212Get rights and content

Highlights

  • A non-local damage model for poro-viscoelastic saturated porous media is proposed.

  • Von-Mises strain measure is used to drive the non-local integral type damage growth.

  • Damage-response includes nonlinear Biot’s coefficients and anisotropic permeability.

  • A mixed FEM solution is developed, including a monolithic iterative solution.

  • Analytically derived Jacobian is used to drive the non-linear solution of equations.

  • Benchmark problems include poro-viscoelastic column and hydraulic fracture problems.

  • Objective mesh-independent and oscillation free numerical results are presented.

  • The model has the potential to represent time-dependent damage in several rock types.

Abstract

We present a novel poro-damage-viscoelastic model for predicting the failure response of fluid-saturated porous geomaterials. The Generalized Maxwell model is introduced for the representation of the viscoelastic behavior of the solid skeleton, which is achieved by a standard Prony-series type expansion. Damage regularization is obtained by an non-local integral-type formulation and damage behavior is described by Mazars model with the modified von Mises-type equivalent strain measure. The poromechanics parameters (Biot’s coefficient, Biot’s modulus) are functions of damage, and the fluid flow obeys Darcy’s seepage law in the entire domain, while the permeability is assumed to be anisotropic and strain dependent. The coupled system is discretized in time using a backward Euler scheme. The non-linear system is linearized using a Newton Raphson scheme and solved monolithically every time step. A consistent Jacobian matrix and residual vector are derived analytically. Several numerical examples are studied in order to investigate the performance of the proposed approach, including (i) a column undergoing hysteresis from cyclic loading, stress relaxation, creep and variable strain rate loading tests and (ii) fluid-driven fracturing in a 2D poro-viscoelastic domain. The numerical time-dependent results exhibit mesh insensitivity for all field variables, and confirm the feasibility and applicability of the proposed non-local damage model for simulating hydraulic fracture.

Introduction

Poromechanics deals with a multi-physics coupled problems involving fluid–saturated porous media. In recent years, poromechanics has become a prominent research topic especially in geomechanics.1,2 The earliest mathematical foundations of poromechanics is attributed to Terzaghi,3 who studied one-dimensional soil consolidation problems. Extension of the work to three-dimensional porous media was later developed by Biot,4 which combined the mechanical equilibrium equation and fluid flow continuity equation, while assuming linear elastic constitutive law for the solid deformation and Darcy’s law for the fluid transport. The first finite element method (FEM) for solving the coupled solid deformation and fluid transport of the saturated porous mixture, was formulated by Zienkiewicz.[5], [6], [7] A thorough review of the mixed finite elements for poroelasticity is presented in Ref. 8. While most applications of poromechanics neglect the time-dependence of the solid skeleton, i.e. assume elastic behavior of the solid, in some applications evidence show that the solid skeleton can creep and deform over time.9,10 Hence, for those applications, it is necessary to introduce viscoelastic behavior of the solid in order to capture more accurately the response of the porous media over time.

To this end, viscoelasticity has been investigated in a variety of topics, for instance, in rock mechanics,[11], [12], [13] geotechnical engineering,14,15 biological engineering,16,17 polymers,18 frozen soils,19 ice[20], [21], [22] and others. In reality, under long-term loading conditions of the porous media, the mechanical behavior of the solid skeleton will be accompanied by a slow creep deformation, and a damage-zone may be initiated leading to a localized failure response.13,23 Furthermore, the interstitial fluid transport is also influenced by creep, stress-relaxation and damage propagation behavior of the solid skeleton, which results in a more complex multi-physics problem.23 Compared with poro-elastic24,25 or poro-elastoplastic26,27 models, poro-damage-viscoelasticity has received far less attention in the literature.28 Accurate modeling of such systems with time-dependent failure is therefore important at the material and structure levels.29,30

A thorough review of integral and differential forms of viscoelastic theories are presented by Vyalov.31 Early forms of viscoelastic constitutive models applied the spring (Hooke model) and dashpot (Newton model) elements32 to construct the Maxwell model33 in series and the Kelvin–Voigt model15 in parallel. Furthermore, research efforts have also explored the use of a combination of these simple units in series and/or parallel to better characterize the viscoelastic behavior of porous media. Examples include the Generalized Maxwell model,[34], [35], [36] Generalized Kelvin–Voigt model37 or the Burgers model.38,39 In this work, the Generalized Maxwell model is utilized to model the viscoelastic behavior, in which the shear components of the deviatoric stress tensor are time dependent while the volumetric part remains time-independent, i.e. purely elastic.22,40 In particular we employ a Prony series-type formulation which provides a convenient a semi-analytical time integration scheme for the stresses, as suggested by Taylor.41 Simo and Hughes.42

Modeling damage behavior in porous media due to mechanical loads, e.g. consolidation or impact problems, or fluid loads, e.g. hydraulic fracture or liquefaction is an important research topic on its own right. In recent decades, many theoretical and computational methods have been presented to simulate damage processes in porous media, which could in general be modeled by discrete or continuum approaches. In the context of discrete models, the discontinuous displacement field to represent the fracture domain is introduced. One of the common options to model cracks is via the extended finite element method (XFEM)[43], [44], [45] or generalized finite element method (GFEM),46,47 which add a set of discontinuous shape functions to the standard parts of FEM. An alternative approach is to employ cohesive zone methods (CZM), which also succeeds in capturing the fracture surfaces and has been widely applied to simulate fracturing processes in porous media.48,49 Nevertheless, one difficulty with discrete crack modeling methods is the need to track complex fracturing patterns, e.g. curved cracks, crack branching, coalescence and crossing, especially in the case of 3D crack surfaces.50,51

In contrast to the discrete approaches, Continuum Damage Mechanics (CDM) is another approach to represent fracture using a narrow continuous band. The earliest efforts in continuum models of fracture are attributed to Rabotnov52 and Kachanov.53 Continuum damage models consider cracks as continuous entities where the damage propagation in the material is identified by a damage variable, as depicted in Fig. 1. Several well known options to represent the damage field, which include one scalar damage,54 two independent scalars[55], [56], [57] or a damage tensor58,59 have been established to model failure in solids by inserting those into the material constitutive law to characterize the degree of stiffness degradation.60 Research efforts have explored the use of classical local damage models for different class of porous media, e.g., in geomaterials61,62 and biomaterials.63

In mathematical terms, the classical local damage leads to a loss of the elliptic character of the set of partial differential equations due to strain softening and lead to a discontinuous rate of deformation and lack of uniqueness.64,65 Further, in this case the numerical solution suffers from dependence on the discretization after reaching a certain damage level.66 In order to overcome the above-mentioned drawbacks introduced by the classical local damage, research efforts have been put forth to develop regularization strategies that restore ellipticity and provide physically meaningful results. Among these methods, integral-64,67 and gradient-type[68], [69], [70], [71] non-local damage models have been extensively developed in recent years. Another closely related technique is the so-called phase field method,[72], [73], [74] which are based on energy minimization of a continuum crack. A comparison between gradient non-local models and phase field models can be found in Refs. 75, 76.

In this research, we adopt the non-local integral-type formulation proposed by Pijaudier-Cabot and Bažant 64 in which the non-local damage variable at a given Gauss point is evaluated as a weighted average of local damage variables within a certain representative volume defined by a characteristic length scale. Non-local regularization methods introduce a characteristic length to effectively alleviate numerical instabilities, spurious post-bifurcation and mesh dependency solutions.65 The non-local integral technique has been widely used to regularize various variables, for instance, damage,77,78 damage rate,79 equivalent strain,80,81 plastic strain,82,83 the rate of plastic multiplier,82 energy dissipation rate,84 strain invariants85 and others. Additionally, an evident advantage of the non-local integral approach is the significant reduction of the computational cost especially for large scale models than that in gradient methods, which require an additional set of degrees of freedom.86

In the context of fluid transport, Darcy’s law87 has been widely employed to describe the fluid flow in porous media. However, a higher speed flow may occur within the fracture domain, which indicate that a constant permeability assumption may be a poor estimate. To that end, an anisotropic strain-dependent non-linear permeability variation of Darcy’s law is adopted.[88], [89], [90] In addition, both Biot’s coefficient and Biot’s modulus in our model are not constant but rather vary with the damage value, which has been shown to be effective in Refs. [91], [92], [93], [94].

In this paper, we aim at developing a novel non-local approach for modeling damage propagation in poro-viscoelastic and fluid-saturated porous media. The paper is organized as follows. Section 2 describes the linear poro-viscoelastic constitutive law and the strain-based non-local integral-type damage approach. In Section 4, the finite element formulation is introduced, and the coupled system of equations is discretized in time using a backward Euler scheme, and solved by an incremental iterative Newton Raphson method for which a consistent Jacobian matrix is derived analytically. Viscoelastic parameters of a typical rock material are determined from experimental tests in Section 5. In Section 6, the performance of the method is examined numerically on a porous media considering hysteresis, creep, relaxation, and strain rates sensitivity. Following, a 2D specimen subjected to fluid-driven fracture are studied to validate the predictive capacity of the proposed model. Finally, some concluding remarks are summarized and future developments are noted.

Throughout this paper, scalars and scalar-valued functions are written in light letters, (e.g., p, D), first and second order tensors and tensor-valued functions are written in bold letters (e.g., u, σ), fourth order tensors with blackboard bold letters (e.g., ̃). Additionally, the following notions of tensorial product are used, i.e., AB=AijBkl, A:B=AijBij. We also define identity tensor as I=δij, I=12(δikδjl+δilδjk), Ivol=13δijδkl, Idev=IIvol, where δij is the Kronecker delta. (), () and s() are the divergence, gradient and symmetric gradient of (), respectively; ()T, ()1 and ()̇ denote the transpose, inverse and time derivative of (), respectively; tr() is the trace operator of (). As for the sign convention, unless otherwise specified, we consider tensile stresses as positive.

Section snippets

Balance laws

Physics of fluid saturated porous media involves the mechanical behavior of a solid skeleton and the fluid percolating throughout the pores. Such problems lead to solid–fluid coupled systems although one does not intend to resolve exactly the interface between the two. Under the actions of mechanical loads or fluid flux excitations, the porous media will deform to a point where a damage-zone may initiate and propagate leading to fluid transport, which likely include fluid transport into the

Poro-viscoelasticity damage relationships

Following previous research efforts,91,93,101 it is worth noting that in the current work Biot’s coefficient and modulus are assumed to vary with the damage growth, such that α(Dˆ)=1K(Dˆ)KsM(Dˆ)=KuK(Dˆ)α(Dˆ)2 where K(Dˆ) is the damaged porous media bulk module defined as K(Dˆ)=(1Dˆ)K, Ks denotes the solid grain bulk modulus. Note that Biot’s coefficient α in Eq. (14) will increase to 1 and Biot’s modulus in Eq. (15) will tend to M=Ku as the damage grows to Dˆ=1, which means that the material

Summary of the boundary value problem

The boundary value problem for poro-damage-viscoelasticity of fluid-saturated porous media, based on the non-local integral-type approach, consist of two partial differential equations: (1) the momentum balance equation (40), and (2) the fluid mass balance equation (41), assuming that inertia terms are negligible. (1Dˆ(t))σ̃(t)α(Dˆ(t))Ip+b=0, in ΩtpM(Dˆ(t))+α(Dˆ(t))tr(ɛ(t))k(κ)p=0, in Ω

Eq. (40) is derived by substituting the definition of the total Cauchy stress σ from Eq. (12) into

Viscoelastic parameter calibration under tri-axial constant stress compression

Identification of parameter values for the Generalized Maxwell model of geomaterials require physical experiments, however, those are rarely reported in the literature. Thus, in order to simplify the parameters choice in the model, we consider an analytical solution for a creep strain response of only the solid-phase (i.e. excluding the fluid and assuming homogeneous solid without damage), in which case we take the Generalized Maxwell model with only one Maxwell element (n=1).

In a conventional

Numerical examples

This section demonstrates the performance of the proposed poro-viscoelastic non-local damage model with strain-dependent anisotropic permeability to capture fracture propagation in fluid-saturated porous media. First, we investigate the viability of the proposed model by means of modeling a polar ice column subjected to cyclic loading, relaxation, creep and various strain rates. Next, the model is used to simulate the hydraulic fracturing process driven by fluid volume injection in porous

Conclusions

We propose a novel non-local integral-type method to simulate fracture failure in poro-viscoelastic media. A Generalized Maxwell model is established and incorporated into constitutive law for solid skeleton behavior to model the time-dependent response. The material deterioration and failure is characterized by a Mazar’s damage law with a modified von Mises-type equivalent strain measure. The fluid flow in the entire domain is described by Darcy’s law with a non-linear strain-dependent

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors are grateful to the funding support provided by the Fundamental Research Funds for the Central Universities (B200203059), Postgraduate Research and Practice Innovation Program of Jiangsu Province (Grant No. KYCX20_0472). The first author is supported by China Scholarship Council (CSC) (201906710104) for his two-year visiting scholar appointment at Columbia University. The authors thank the anonymous reviewers for their constructive comments to improve the quality of the paper.

References (154)

  • SongH. et al.

    Numerical modelling of hydraulic fracture propagation in poro-viscoelastic formation

    J Pet Sci Eng

    (2021)
  • NguyenS.-T. et al.

    A homogenization approach for effective viscoelastic properties of porous media

    Mech Mater

    (2016)
  • NguyenS.-T.

    Effect of pore shape on the effective behavior of viscoelastic porous media

    Int J Solids Struct

    (2017)
  • BaoT. et al.

    Experimental workflow to estimate model parameters for evaluating long term viscoelastic response of CO2 storage caprocks

    Int J Rock Mech Min Sci

    (2021)
  • BaxevanisT. et al.

    Bifurcation and creep effects in a viscoelastic non-local damageable continuum

    Eur J Mech A Solids

    (2008)
  • AroraK. et al.

    Viscous-elastic-plastic response of tunnels in squeezing ground conditions: Analytical modeling and experimental validation

    Int J Rock Mech Min Sci

    (2021)
  • ShenR. et al.

    Fracture of viscoelastic solids modeled with a modified phase field method

    Comput Methods Appl Mech Engrg

    (2019)
  • Parchei-EsfahaniM. et al.

    Dynamic hydraulic stimulation and fracturing from a wellbore using pressure pulsing

    Eng Fract Mech

    (2020)
  • KomijaniM. et al.

    Enriched mixed finite element models for dynamic analysis of continuous and fractured porous media

    Comput Methods Appl Mech Engrg

    (2019)
  • CarrierB. et al.

    Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model

    Eng Fract Mech

    (2012)
  • GuptaV. et al.

    Stable GFEM (SGFEM): Improved conditioning and accuracy of GFEM/XFEM for three-dimensional fracture mechanics

    Comput Methods Appl Mech Engrg

    (2015)
  • YiL.-P. et al.

    A fully coupled fluid flow and rock damage model for hydraulic fracture of porous media

    J Pet Sci Eng

    (2019)
  • Le BellégoC. et al.

    Calibration of nonlocal damage model from size effect tests

    Eur J Mech A Solids

    (2003)
  • CerveraM. et al.

    Cracking of quasi-brittle structures under monotonic and cyclic loadings: A d+/d- damage model with stiffness recovery in shear

    Int J Solids Struct

    (2018)
  • LubardaV. et al.

    Damage tensors and the crack density distribution

    Int J Solids Struct

    (1993)
  • SaharaD.P. et al.

    Analysis of borehole breakout development using continuum damage mechanics

    Int J Rock Mech Min Sci

    (2017)
  • LeiQ. et al.

    Modelling fluid injection-induced fracture activation, damage growth, seismicity occurrence and connectivity change in naturally fractured rocks

    Int J Rock Mech Min Sci

    (2021)
  • JeffersJ.R. et al.

    Damage accumulation, fatigue and creep behaviour of vacuum mixed bone cement

    Biomaterials

    (2005)
  • PeerlingsR. et al.

    A critical comparison of nonlocal and gradient-enhanced softening continua

    Int J Solids Struct

    (2001)
  • PeerlingsR. et al.

    A thermodynamically motivated implicit gradient damage framework and its application to brick masonry cracking

    Comput Methods Appl Mech Engrg

    (2004)
  • MobasherM.E. et al.

    Dual length scale non-local model to represent damage and transport in porous media

    Comput Methods Appl Mech Engrg

    (2021)
  • MieheC. et al.

    A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits

    Comput Methods Appl Mech Engrg

    (2010)
  • YuZ. et al.

    Numerical study of thermo-hydro-mechanical responses of in situ heating test with phase-field model

    Int J Rock Mech Min Sci

    (2021)
  • de BorstR. et al.

    Gradient damage vs phase-field approaches for fracture: Similarities and differences

    Comput Methods Appl Mech Engrg

    (2016)
  • JamesK.A. et al.

    Failure mitigation in optimal topology design using a coupled nonlinear continuum damage model

    Comput Methods Appl Mech Engrg

    (2014)
  • De VreeJ. et al.

    Comparison of nonlocal approaches in continuum damage mechanics

    Comput Struct

    (1995)
  • TotiJ. et al.

    Nonlocal damage propagation in the dynamics of masonry elements

    Comput Struct

    (2015)
  • Pijaudier-CabotG. et al.

    Finite element analysis of bifurcation in nonlocal strain softening solids

    Comput Methods Appl Mech Engrg

    (1991)
  • ComiC.

    A non-local model with tension and compression damage mechanisms

    Eur J Mech A Solids

    (2001)
  • SelvaduraiA. et al.

    Mandel–Cryer effects in fluid inclusions in damage-susceptible poroelastic geologic media

    Comput Geotech

    (2004)
  • SelvaduraiA.

    Stationary damage modelling of poroelastic contact

    Int J Solids Struct

    (2004)
  • MobasherM.E. et al.

    Thermodynamic framework for non-local transport-damage modeling of fluid driven fracture in porous media

    Int J Rock Mech Min Sci

    (2018)
  • ShaoJ.

    Poroelastic behaviour of brittle rock materials with anisotropic damage

    Mech Mater

    (1998)
  • ChengA.-D.

    Material coefficients of anisotropic poroelasticity

    Int J Rock Mech Min Sci

    (1997)
  • MobasherM.E. et al.

    Non-local formulation for transport and damage in porous media

    Comput Methods Appl Mech Engrg

    (2017)
  • YiL.-P. et al.

    A consistent phase field model for hydraulic fracture propagation in poroelastic media

    Comput Methods Appl Mech Engrg

    (2020)
  • PolizzottoC.

    Unified thermodynamic framework for nonlocal/gradient continuum theories

    Eur J Mech A Solids

    (2003)
  • BorinoG. et al.

    A symmetric nonlocal damage theory

    Int J Solids Struct

    (2003)
  • De BoerR.

    Trends in Continuum Mechanics of Porous Media, Vol. 18

    (2006)
  • CoussyO.

    Mechanics and Physics of Porous Solids

    (2011)
  • Cited by (0)

    View full text