Non-local continuum damage model for poro-viscoelastic porous media
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 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., , ), first and second order tensors and tensor-valued functions are written in bold letters (e.g., , ), fourth order tensors with blackboard bold letters (e.g., ). Additionally, the following notions of tensorial product are used, i.e., , . We also define identity tensor as , , , , where is the Kronecker delta. , and are the divergence, gradient and symmetric gradient of , respectively; , and denote the transpose, inverse and time derivative of , respectively; 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 where is the damaged porous media bulk module defined as , 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 as the damage grows to , 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.
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 .
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)
- et al.
Modeling of hydraulic fracturing in viscoelastic formations with the fractional Maxwell model
Comput Geotech
(2020) - et al.
Experimental study on creep of double-rock samples disturbed by dynamic impact
Int J Rock Mech Min Sci
(2021) - et al.
Time-dependent deformation of shale gas reservoir rocks and its long-term effect on the in situ state of stress
Int J Rock Mech Min Sci
(2014) - et al.
Experimental investigations of the creep–damage–rupture behaviour of rock salt
Int J Rock Mech Min Sci
(2014) - et al.
A temperature dependent creep damage model for polycrystalline ice
Mech Mater
(2012) - et al.
A prony-series type viscoelastic solid coupled with a continuum damage law for polar ice modeling
Mech Mater
(2016) - et al.
An equivalent stress-gradient regularization model for coupled damage-viscoelasticity
Comput Methods Appl Mech Engrg
(2017) - et al.
A thermodynamic framework for constitutive modeling of coupled moisture-mechanical induced damage in partially saturated viscous porous media
Mech Mater
(2016) - et al.
Determination of the main hydraulic and poro-elastic properties of a limestone from Bourgogne, France
Int J Rock Mech Min Sci
(2004) - et al.
Analytical and experimental investigation of pore pressure induced strain softening around boreholes
Int J Rock Mech Min Sci
(2019)