FE Analysis of Rock with Hydraulic-Mechanical Coupling Based on Continuum Damage Evolution

A numerical finite element (FE) analysis technology is presented for efficient and reliable solutions of rock with hydraulicmechanical (HM) coupling, researching the seepage characteristics and simulating the damage evolution of rock. To be in accord with the actual situation, the rock is naturally viewed as heterogeneous material, in which Young’s modulus, permeability, and strength property obey the typical Weibull distribution function. The classic Biot constitutive relation for rock as porous medium is introduced to establish a set of equations coupling with elastic solid deformation and seepage flow. The rock is subsequently developed into a novel conceptual and practical model considering the damage evolution of Young’s modulus and permeability, in which comprehensive utilization of several other auxiliary technologies, for example, the Drucker-Prager strength criterion, the statistical strength theory, and the continuum damage evolution, yields the damage variable calculating technology. To this end, an effective and reliable numerical FE analysis strategy is established. Numerical examples are given to show that the proposed method can establish heterogeneous rock model and be suitable for different load conditions and furthermore to demonstrate the effectiveness and reliability in the seepage and damage characteristics analysis for rock.


Introduction
The finite element (FE) analysis technology presented in this paper takes rock as the research object, a typical porous medium with heterogeneous material property, dispersedly distributed damage region under applied load and in situ stress, and inorganic and organic pores to form the complex structure [1].
Heterogeneous distribution and degeneration of material property (i.e., Young's modulus and permeability) or strength property (i.e., compressive and tensile strength) make the assumption of homogeneousness and continuity is not consistent with the actual situation; for some extreme cases, the state-of-the-art strength theory of rock based on continuum mechanics also confronts severe problems challenging the researchers' best knowledge [2].Some technology, such as heterogeneous Young's modulus simulation to describe the heterogeneity of rock, is developed [3], and the authors of this paper also have adopted this technology to form the heterogeneous rock with heterogeneous Young's modulus and strength [4].The technology could be developed; the permeability is heterogeneous and has some relationship with the heterogeneity of Young's modulus to form a new technology to describe the heterogeneity of rock in this paper.
Unconventional rock, such as shale, has the characteristic of low porosity and permeability; particularly the key technology of hydraulic fracking could create artificial fracture network with centralized distribution to increase the porosity and permeability for large-scale extraction [5].The effective stress of porous medium will change with the fluid flow, pressure diffusion in pores, and solid deformation correspondingly; in other words, hydraulic-mechanical (HM) coupling makes the response of the rock reflect complex timedependent effect obviously.Also, for research, the singleporosity elastic model and dual-porosity elastic model [6] are developed successively; additionally, FE analysis technology 2 Mathematical Problems in Engineering for a new constitutive model considering the adsorption has been put forward recently [7].
For some conventional and simple problems, such as possessing regular solving domain with homogeneous material property, analytical method via the poroelastic theory has been proposed to solve a number of multifields coupling, wellbore stability, and fracture problems [8,9], however, which unfortunately does not consider heterogeneity, plastic deformation, and damage evolution and hence cannot serve as a practical and general method.Therefore, the authors of this paper were motivated to probe into this numerical area of rock fracture and damage evolution problems.The discrete fracture network (DFN) model [10] based on the fracture mechanics reflects both effective and efficient solving ability for single physical field and regular fracture system, belonging to this model; a typical extended finite element method (XFEM) has successfully developed.The authors of this paper have made a series of achievements on the XFEM theory and algorithm [11], that is, the problems of crack growth in pipes [12], multiscales computation [13], two-phase flows [14], error estimation [15], cracks in shells [16], and crack crossing [17]; additionally, the XFEM have been improved and applied to numerical simulation for hydraulic fracturing of shale rock in recent study [18].On the other hand, some researchers proposed the continuous model based on the damage mechanics to analyze the failure process, compared to the above DFN model, which more easily handles the multifields coupling and complex medium containing pores and cracks problems to be consistent with the material behavior [19].Inspired from that, some numerical methods and simulation technologies based on damage analysis are proposed [3,20].The damage variable that described the compressive and tensile strength of rock was developed by mesoscopic damage mechanics in [3]; a novel, simple FE analysis technology of seepage in rock has been presented by the authors of this paper with the further development of the continuous model and continuum damage variable [4].In this paper, the continuum damage variable was developed by introducing a correction factor to consider the part of the bearing capacity of damage rock effectively.
The well-known Biot constitutive theory [21,22] is introduced firstly as the basis of the proposed method; thereafter, solid and gas seepage equations coupled by the term of pore pressure are established to form the HM coupling.The FE analysis technology in this paper achieves the results for coupled HM model of rock based on continuum damage evolution by implementing the following three-step strategy.In the first step, for establishing physical HM coupling model, the heterogeneous rock is composed as Young's modulus and permeability on each FE node obey the Weibull distribution.Then FE solutions of HM model, such as effective stresses, permeability, and porosity, are obtained by a developed FE procedure in the second step.In the last step, the above effective stress solutions are used to calculate the damage variable on each FE node, and furthermore Young's modulus and permeability with damage on each FE node could be calculated, respectively.Young's modulus and permeability with damage under the current load would be used to form a new heterogeneous physical model; the procedure returns to the second step until the load steps are progressively completed.This yields a simple, efficient, reliable, and practical FE analysis technology that is able to make rock modeling and damage analysis in mining and petroleum engineering applications.
The numerical examples presented later have shown that utilizing petrophysical heterogeneity simulation is well effective and feasible; in addition to the reliability, the results also show that the method is almost consistent with damage evolution under different load conditions, considering HM coupling is necessary for effective stress analysis.

Assignment of Petrophysical Heterogeneity
Heterogeneity of rock is research focus with some simulated difficulties currently; taking simplifications of homogenizing treatment, obviously deviated results to the actual situation will be got.Various mineral particles and cementing materials constitute the common rock, in which forming process is extremely complex.Dispersedly distributed pores and microcosmic cracks in the rock exhibit inhomogeneity and heterogeneity of physical and mechanical properties.Petrophysical heterogeneity analysis widely adopts the Weibull distribution [23] to research the structural strength, fatigue, and other problems of rock, which have obtained some satisfactory numerical simulation effect and the good consistency with the experiments [3,24,25].
In order to characterize the heterogeneity of rock, in this study, the model is discrete as many elements by FEM, and the mechanical properties of each node on FE mesh are assumed to conform to a given Weibull distribution as defined in the following probability density function: where  is a variable representing the mechanical parameter (such as material modulus or strength),  is scale parameter representing the average range of , and  is shape parameter characterizing the homogeneity degree of .
From the properties of the Weibull distribution, a larger value of parameter  implies a more homogeneous material and vice versa.Therefore, the parameter  can be called the homogeneity index, for higher values of the homogeneity index, the material property, as variable  of more elements are concentrated closer to .An increase inhomogeneity index leads to more homogeneous numerical model.Utilizing of (1) in computer, simulation of a medium composed of many elements, one can produce numerically a heterogeneous rock model.The computationally produced heterogeneous model is analogous to a real specimen tested in laboratory, so in this investigation it is referred to as a numerical rock model.In this respect, the influence of heterogeneity of material properties can be examined based on the numerical simulations.
The material property (i.e., Young's modulus and permeability) or strength property could describe the heterogeneity of rock, which obey the above Weibull distribution in this paper.Utilizing of the Monte Carlo method, each node on FE mesh is assigned a random value   ( = 1, . . .,   ) within [0, 1] as the value of the probability distribution function of Weibull distribution; then the random value   could be obtained according to   ; here   is the total number of nodes.
According to the physical property of rock, the region with bigger value of Young's modulus could have weak permeability characteristic conversely in general cases.So the permeability is considered to have inversely proportional relationship with Young's modulus in this paper; in other words, if Young's modulus on one node is   computed by Weibull distribution, then the corresponding permeability on this node is   =  0 /(  / 0 ); here  0 and  0 are initial Young's modulus and permeability, respectively.

Governing Equations with Hydraulic-Mechanical Coupling
Rock is the typical porous medium in which seepage characteristic as an important factor affects the solid deformation and gas flow.In another aspect, rock has the obvious nonlinear characteristic as plastic deformation; however the nonlinear model will increase difficulty and cost of calculation with low numerical accuracy sometimes, so the linear elastic models combined with damage theory instead to describe the nonlinear behavior which is also adapted in this paper.The Biot constitutive relation is a well-developed and practical constitutive relation for porous and linear elastic medium, by which governing equations with HM coupling are derived as below.
The effective stress of rock could be introduced by considering pore pressure of gas with Biot constitutive theory; further utilizing of mechanical equilibrium and deformation compatibility conditions, the elastic deformation equation of rock could be obtained as where  is the displacement,  is gas pressure,  , is the component of the solid body load,  is Young's modulus, ] is Poisson's radio,  is effective stress parameter, and Ω is the solving domain.The boundary conditions (BCs) of ( 2) are given as =  0 () , where  0 and  0 are initial displacement and stress, respectively, in Ω and  0 () and   () are given displacement and stress, respectively, associated with time parameter  on Ω.
For simplicity, both subscripts  and  are taken as 1, 2, and 3 without special statement in this paper below.The hypothesis of gas in saturated state satisfying the gas mass conservation is taken in this paper; further utilizing of Darcy's law and the porosity have the cubic relationship with permeability according to experiments [26]; the control equation of gas seepage flow could be expressed as where  is porosity,  is permeability, f ℎ is hydraulic body load,  is kinematic viscosity coefficient,  V is the volumetric strain,   is gas quality source,   is /3(1 − 2), and  is  V + /  .The BCs of ( 4) are given as where  0 is initial gas pressure in Ω and () and () are given pressure and gas flow, respectively, with time parameter  on Ω.

Continuum Damage Evolution
In actual engineering, rock is typical brittle material in which damage will occur in some micro-units and expand continuously under load.Utilizing of continuum damage mechanics and statistical strength theory, the damage variable  can be defined by failed number   and total number  of micro-units as   /.In the conventional understanding, the micro-unit bearing capacity will lose once the damage happens, but in the practical cases, the bearing capacity decreases partly and also remains keeping capacity to bear the load.So a correction factor , which has the range of [0, 1], is introduced to reflect the above mechanical property; under the current external load, the failed number   of microunits is converted to   .The Drucker-Prager (DP) strength criterion is based on Mohr-Coulomb (MC) strength criterion, in which the inference of the intermediate principal stress and hydrostatic pressure is considered.The DP strength criterion ℎ() −   =  1 + √ 2 −   = 0 is used to judge the failure behavior; here  and   are real numbers,  1 is the first invariant of the stress tensor, and  2 is the second invariant of the stress deviator.So combined value of stress ℎ() can represent the current load status; according to the hypothesis of the micro-unit strength obeying the Weibull distribution, the damage variable could be derived as where parameter  is the rock compressive or tensile strength.When the micro-unit is in the dominant state of compression as 0 ≤ The permeability of rock will improve gradually due to the macroscopic damage generation and expansion; based on relationship between permeability and strain, an effective technology of permeability damage evolution is proposed [3,4] to reflect the damage characteristic as where k is permeability with damage,  0 and  0 are initial permeability and porosity, respectively,  is the permeability damage coefficient, and exp() is the continuum damage item to describe that the permeability is enhanced by damage.Through some numerical examples test including the results given below, ( 7) and ( 8) describing the damage of rock reveal excellent effect.

FE Analysis Strategy
The FE analysis in this paper achieves the results for coupled HM model of rock based on continuum damage evolution by implementing the following three-step strategy.
(1) Heterogeneous Model.The parameters  and  in Weibull distribution are given by user according to the heterogeneous property of rock; moreover the FE model of rock is established for generating the FE mesh with each node possessing the different value for Young's modulus  and permeability  to derive the heterogeneous rock, as described in Section 2.
(2) FE Solution.On the current load step, the effective stress   and gas pore pressure  could be calculated using ( 2) and ( 4) to obtain solutions of HM model by the standard FE analysis, described in Section 3.
(3) Damage Analysis.The above effective stress   is used to calculate the damage variable  on each node by ( 6), and furthermore Young's modulus Ẽ and permeability k with damage on each node could be calculated by ( 7) and ( 8), respectively, forming a new FE model, as described in Section 4. Then the procedure returns to the second step (i.e., FE solution) until the load steps are progressively completed.
The above three steps constitute a round of continuum damage analysis for rock with HM coupling, by which the gas seepage flow and rock deformation results could be obtained.
The numerical examples below show the proposed FE analysis strategy is effective and reliable.

Numerical Examples
The proposed analysis strategy has been coded into a MAT-LAB program and partly used the FE solver of COMSOL Multiphysics to obtain solutions of governing equations with HM coupling.This section presents three interrelated numerical examples showing the excellent performance of the procedure.Throughout, the program is run on a Dell OptiPlex 380 Intel5 Core6 2.93 GHz desktop computer.
The first example is chosen to discuss the appropriate parameter selection of Weibull distribution to describe the heterogeneity for rock, so that the rock model researched below with heterogeneous Young's modulus and permeability will be established reliably.The second example analyzes two damage evolution cases of heterogeneous rock with an inner circular hole under different load conditions, in which effective stress results of one typical model will be specially researched as the third example: considering HM coupling is compared with the single solid model without considering rock as pore structure.In the last two examples, the triangle element is used and the shape parameters  of Weibull distribution for heterogeneous material property and strength property are set to be the same throughout.
Example 1 (FE model of heterogeneous rock).Consider plane strain model in Figure 1; the geometric domain is square possessing an inner circular hole with  where  is side length of the square,  is diameter of the hole, and  is the thickness parameter.
The next discusses the influence of Weibull distribution parameters  and  on petrophysical heterogeneity.Take  = 1 and  = 2, 4, 6 as examples, the distribution results obtained by the proposed method as Figure 2. The distributed values in the entire domain are from 0.11 to 2.24, from 0.21 to 1.51, and from 0.38 to 1.34 in Figures 2(a), 2(b), and 2(c), respectively, in which it can be seen that the distributed values are getting closer together to the mean value  with homogeneity index parameter  taking bigger value.These results are well consistent with the Weibull distribution theory as introduction in Section 2.
Considering the above influence of Weibull distribution parameters,  = 5 is adapted to simulate the heterogeneous material parameter as a moderate distribution way to let values be neither too big nor too small, and  =  0 = 16.32 GPa is adapted as the initial, mean value of Young's modulus.As shown in Figure 3, the distributed result of Young's modulus   1, some of which are chosen from [7,27].Two different load conditions for damage evolution are analyzed comparatively: model I in Figure 5 has uniform extensile load  on the upper and lower edges and model II in Figure 6 has uniform pressure load  in inner circular hole and simply supported boundary outside.The computed damage evolution results of two different load conditions as load  from 0 MPa to 50 MPa are shown in Figures 7 and 8, in which region near the hole is selected because the damage variable  has relatively dramatic changes there.As shown in Figures 7 and 8, the damage regions expand discontinuously with  increasing as the heterogeneity of rock.In Figure 7, the damage region the hole expands gradually to the horizontal direction of the model by stress consideration, while, in Figure 8, the damage region around the hole expands outward consistently without significant damage propagation.These above results are consistent with damage theory and regular experiments to evaluate the effectiveness of rock damage analysis of the proposed method.
Example 3 (effective stress analysis by HM coupling influence).The effective stress characteristic of rock HM coupling will be analyzed by the proposed method, utilizing plane strain model II in Figure 6 widely used as wellbore model in mining and petroleum engineering.For comparative analysis, the stress result of single solid model was computed, which has the same geometry and physical parameters, considering petrophysical heterogeneity and damage analysis as HM coupling case in Figure 6.The radial and tangential stress results at cross section A-A in Figure 6 were analyzed as  = 50 MPa in Figures 9 and 10, respectively, where  is the radius of the wellbore and  is the distance from the wellbore center.As shown in Figure 9, it can be seen that both radial stress of HM model and single solid model are almost the same at near field and far field of wellbore, but the stress result at the middle field of HM model is bigger than solid model.The tangential stress is a decisive factor for wellbore stability, as shown in Figure 10, in which result of single solid model at near field is significantly smaller than HM coupling's; then the wellbore instability would tend to happen in single solid model case.Therefore, multifields coupling model and the effective stress analysis technology is significant and should be developed reasonably.

Conclusions
HM model of heterogeneous rock based on continuum damage evolution and FE analysis technology with suitable analysis strategy for effective and reliable computation of the seepage and damage of rock have been presented.Comprehensive utilization of continuum damage theory with a number of other auxiliary techniques (including the Weibull distribution, classic Biot constitutive relation, and DP strength criterion) has yielded a simple, efficient, and reliable FE procedure that describes the damage behavior for being consistent with the property of rock to be achieved.The techniques integration removes major difficulties as heterogeneity, multifields coupling, and damage evolution in the solution of the seepage flow, collapse, and fracture problems of rock, which enhances the overall efficiency and reliability of solution.Results for typical numerical examples, including ones known to be engineering problems, have shown that the present method is competitive and possesses the potential for further extension to more complex mechanic problems (e.g., anisotropy damage model and wellbore stability problems).A recent study of transversely isotropic hydraulic-mechanicaldamage coupling model successfully solved wellbore stability, as well as rock instability problems, so paving the way for corresponding anisotropic damage tensor and failure problems of brittle material in complex conditions to be solved

Figure 7 :
Figure 7: Damage variable  of model I with outer tensile load.
(6)  2 ≤  1 or  1 ≥ 0,  3 ≤ 0, and | 1 | > | 3 |,  takes the compressive strength ; on the other hand, when the micro-unit is in the dominant state of tension as  3 ≤  2 ≤  1 ≤ 0 or  1 ≥ 0,  3 ≤ 0, and | 3 | ≥ | 1 |,  takes the tensile strength .Here  1 ,  2 , and  3 are three principal stresses, whose sign convention used throughout this paper is that compressive stress and strain are positive.As shown in(6), the damage variable  is associated with both the rock strength and the stress state under the current load.It is noted that the damage variable  of micro-unit keeps the original value in the unloading case.

Table 1 :
Physical parameters of model.