Three-Dimensional Numerical Investigation of Coupled Flow-Stress-Damage Failure Process in Heterogeneous Poroelastic Rocks

The failure mechanism of heterogeneous rocks (geological materials), especially under hydraulic conditions, is important in geological engineering. The coupled mechanism of flow-stress-damage should be determined for the stability of rock mass engineering under triaxial stress states. Based on poroelasticity and damage theory, a three-dimensional coupled model of the flow-stress-damage failure process is studied, focusing mainly on the coupled characteristics of permeability evolution and damage in nonhomogeneous rocks. The influences of numerous mesoscale mechanical and hydraulic properties, including homogeneity, residual strength coefficient, loading rates, and strength criteria, on the macro mechanical response are analyzed. Results reveal that the stress sensitive factor and damage coefficient are key variables for controlling the progress of permeability evolution, and these can reflect the hydraulic properties under pre-peak and post-peak separately. Moreover, several experiments are conducted to evaluate the method in terms of permeability evolution and failure process and to verify the proposed two-stage permeability evolution model. This model can be used to illustrate the failure mechanics under hydraulic conditions and match different rock types. The relation of permeability with strain can also help confirm appropriate rock mass hydraulic parameters, thereby enhancing our understanding of the coupled failure mechanism in rock mass engineering.


Introduction
The failure mechanism under hydraulic conditions in heterogeneous geotechnical materials is a significant issue and challenge in petroleum and mining engineering [1][2][3]. Failure initiation, propagation, coalescence, and out of control in-rock engineering are vital to ensuring safe and efficient production [4][5][6]. When new failures occur, the main migration pathways are formed and easily cause complicated dynamic mutation of the fluid flow, which should be considered seriously to prevent the induction of water inrush, gas outburst, and other dynamic disasters [7]. Experimental research results show that rock permeability is a function of confining pressures, pore pressures, stresses and Pore pressure only affects the mechanical responses of dilated and damaged rock salt [41], however, the rules of damage development and evolution and the relationship between permeability and damage on poroelasticity need further study. The poroelastic solutions for rock specimens under hydrostatic and triaxial conditions are necessary for precise predictions of the macroscopic mechanical behavior of the saturated samples.
The current study mainly focused on the coupling characteristics of flow-stress-damage in heterogeneous rocks based on poroelasticity. A 3D numerical model was established for study by using FEM with the software Comsol Multiphysics with Matlab, considering the mechanism of flow-stress-damage and different heterogeneity degrees, loading rates (displacement per unit time, LR), residual strength coefficients (the ratio of residual strength to initial strength, RSC), and strength criteria. Moreover, indoor rock mechanical tests were carried out to enhance the theory and method in terms of the failure process coupled with permeability evolution and to verify and compare them with the proposed model. This model could be applied to reveal the propagation of fractures and the coupled fluid flow and to illustrate the macroscopic mechanical characteristic of rock masses according to accumulative microcosmic failure. It is expected through the investigation to obtain profound understanding of the coupled flow-stress-damage mechanism and to provide a theoretical foundation for rock mass engineering application.

Mechanical Equilibrium Equation
The rock is assumed to behave as a linear elastic medium based on poroelasticity theory. The static equilibrium equation of the porous medium can be governed by the following: where, σ ij,j = ∂σ ij /∂x j is the component of total stress tensor and F i is the component of the body force in the ith direction. The constitutive relation for the deformed porous medium (positive for tension) can be expressed as where K and G are the bulk modulus and shear modulus of rocks, respectively. α is the Biot's coefficient which can be obtained by experiments, p represents the fluid pressure in the pore and δ ij represents the Kronecker delta, which is 1 for i = j and 0 for i = j.

Failure Criterion and Damage Evolution Equation
Based on damage mechanics theory, the elastic modulus of rock degrades gradually as damage occurs, which can be expressed as (as shown in Figure 1): where E and E 0 are the elastic modulus of the damaged and undamaged elements, respectively. D represents the damage variable. Here, the element and its damage are assumed to be isotropic, consequently, E, E 0 , and D are scalars. In this study, the maximum tensile stress criterion and the Mohr-Coulomb (MC) criterion are used to evaluate the tension damage and shear damage, respectively, which can be expressed as Equation (4) [42,43].
Energies 2018, 11, 1923 4 of 16 where σ 1 is the first principal stress, σ 3 is the third principal stresses, f t is the uniaxial tensile strengths, f c is the uniaxial compressive strengths, ϕ is the internal frictional angle, and F 1 and F 2 are tensile and shear damage threshold functions, respectively. Here, the tensile stress criterion is used preferentially, and is initially applied to evaluate whether elements are damaged under tension; only the undamaged elements under the tensile conditions are checked for shear damage by using the MC criterion [6]. Furthermore, the Drucker-Prager (DP) criterion does not consider the effects of the invariant J 3 on the cross-sectional shape of the yield surface. This criterion is a smooth function based on the invariants I 1 and J 2 together with two other constants α dp and K dp , as stated in Equation (5). This criterion is sometimes also named the extended Von Mises criterion because it is equivalent to the Von Mises criterion when setting α dp = 0.
Fα dp I 1 + J 2 − K dp = 0, where F is the damage threshold functions, α dp and K dp are the coefficients of DP strength criterion, and α s is the influence coefficient of the intermediate principal stress.
and is initially applied to evaluate whether elements are damaged under tension; only the undamaged elements under the tensile conditions are checked for shear damage by using the MC criterion [6]. Furthermore, the Drucker-Prager (DP) criterion does not consider the effects of the invariant J3 on the cross-sectional shape of the yield surface. This criterion is a smooth function based on the invariants I1 and J2 together with two other constants αdp and Kdp, as stated in Equation (5). This criterion is sometimes also named the extended Von Mises criterion because it is equivalent to the Von Mises criterion when setting αdp = 0. 1 2 0 − dp dp where F is the damage threshold functions, αdp and Kdp are the coefficients of DP strength criterion, and αs is the influence coefficient of the intermediate principal stress.   [44][45][46][47]. The DP and unified strength criteria (twin shear (TS) criterion) [47], both of which consider all principal stresses, are adopted instead of F2 with F in this study. The MC criterion is a popular criterion in soil/rock mechanics. This criterion states that failure occurs when the normal and the shear stresses acting on any element in the geotechnical material satisfy a yield function (not considering intermediate principal stress and other shear stress). The MC criterion provides an irregular hexagonal pyramid in the principal stress space, generating singularities in the derivatives of the yield function. The MC criterion causes numerical difficulties when treating plastic flow at the corners of the yield surface. The TS criterion considers all stress components causing the failure of any element. This criterion satisfies the property of geotechnical materials and provides a good solution to the tension problem, and the compression meridian limits are unequal in the geotechnical materials.
As illustrated in Figure 1, the damage variable can be expressed as   [44][45][46][47]. The DP and unified strength criteria (twin shear (TS) criterion) [47], both of which consider all principal stresses, are adopted instead of F 2 with F in this study. The MC criterion is a popular criterion in soil/rock mechanics. This criterion states that failure occurs when the normal and the shear stresses acting on any element in the geotechnical material satisfy a yield function (not considering intermediate principal stress and other shear stress). The MC criterion provides an irregular hexagonal pyramid in the principal stress space, generating singularities in the derivatives of the yield function. The MC criterion causes numerical difficulties when treating plastic flow at the corners of the yield surface. The TS criterion considers all stress components causing the failure of any element. This criterion satisfies the property of geotechnical materials and provides a good solution to the tension problem, and the compression meridian limits are unequal in the geotechnical materials.
The rock specimen in this study is assumed to comprise several mesoscopic elements to characterize rock heterogeneity, and the mechanical characteristic of these elements are assumed to follow a typical Weibull distribution, which is defined in Equation (9) [43]: where u and u 0 are the mechanical parameter of the elements and the scale parameter related to the material parameters average, respectively. m is the shape parameter reflecting the homogeneity degree and is defined as a homogeneity index. The value of m illustrates the statistics rule of the mechanical parameter.

Flow Equation
The fluid is assumed to conform to Darcy's law, and the mass conservation equation for fluid flow in porous rock material can be expressed as follows [39,48]: ∂p ∂t where k represents the rock permeability, ρ is the fluid density at the in-situ condition, which can be considered as incompressible, q and p represent the velocity vector and pore fluid pressure, µ and g are the fluid dynamic viscosity and gravity, respectively; t is the real time, Z is elevation, M represents the Biot modulus, α is the Biot-Willis coefficient and is the change rate in the volumetric strain of the porous matrix. The right-hand term of Equation (11) can be interpreted as the expansion rate of the pore space. As it increases, the volume fraction available for the fluid also increases, thereby allowing liquid to sink. Finally, the undrained Poisson ratio ν u should be noted. The drained shear modulus is identical to the undrained shear modulus; thus, the drained and undrained Poisson's ratios are [48]: where K u represents the undrained bulk modulus of rock. The Biot modulus M can be expressed as follows: According to the fundamentals of poroelasticity [48], porosity φ is related to the effective confining pressure P and the shear modulus of the solid phase G s . φ 0 is the porosity at the unstressed state. Mutation coefficient ξ is introduced in Equation (15) to express permeability change due to Energies 2018, 11,1923 6 of 16 damage [49,50], β represents the stress sensitive factor, and can be measured in-situ and in the laboratory on three different scales, more details are presented in Louis' work [51].
Rock permeability can be relevant to its porosity; the generalized cubic relationship of permeability and porosity can also be used for matrices [52][53][54][55][56], as shown in Equation (16).
Although the permeability is hereafter considered to be a stress-independent functional relationship, the above relationship is modified considering the flow-stress-damage process. Equation (17) indicates the exponential function for the absolute flow permeability considering porosity and damage [57]: where k 0 is the initial permeability, and k is the current permeability under non-zero stress conditions. In this study, the stress-permeability relationship is developed based on indoor tests and the combination of Equations (15) and (17) yields: The damage coefficient α D and the stress sensitive factor β are the key variables for controlling the progress of permeability evolution, which is a complicated and coupled progress; the damage coefficient α D shows the increases in permeability, whereas the stress sensitive factor β state the decrease in permeability. The results of our experimental study confirm the previous statement on the topic, and match well with the Equation (18), more details will be discussed later. In other words, the pre-peak and post-peak hydraulic properties are described by different parameters separately.

Rock Specimens and Testing Procedure
To investigate the influence of different test parameters (pore and confining pressure) on the permeability evolution during the complete stress-strain curve, we conducted the experiments using an electro-hydraulic servo-controlled rock mechanics testing system (MTS) 815.02. The maximum displacement and loading capacity are 25 mm and 2700 kN, respectively. The test system contains two subsystems, namely, axial load and pore water pressure load systems. The test process and detailed structures of the triaxial cell are shown in Figure 2. This servo-controlled system have two typical control (load and displacement loading) while the test data are recorded automatically. In the experiment, rock samples were tested at hydrostatic stresses of 6 and 8 MPa with a certain pore pressure gradient, and then the load was increased with the LR of approximately 0.005 mm/step to failure under conventional triaxial compressive stresses. Loads, deformations, and permeability of tested sandstone samples were recorded simultaneously. The experimental method adopted the transient method, and the implementation steps were as follows: (1) the cylindrical rock specimens were oven dried for 24 h or vacuumed and then saturated with water for 48 h before application. The sample was installed onto the base plate, and the piston was maintained in contact with the upper surface of sample, and then the triaxial cell was dropped, and oils were injected into the cell. (2) Axial pressure (P 1 ), confining pressure (P 2 ), and hole pressure (P 3 , P 4 ) were controlled by the computer to reach and maintain certain values. (3) The experiment was started with axial compression applied to the first experimental point, with P 2 , and P 3 remaining constant at the lower side of the pore pressure Energies 2018, 11,1923 7 of 16 P 4 , seepage results under pressure difference ∆P = P 3 -P 4 , and the ∆P-t curve was plotted. (4) Loading was increased to the next experimental point.
Meanwhile, mercury intrusion porosimetry (MIP) was conducted for determining the porosity of rock specimens, and scanning electron microscopy (SEM) was carried out for observing the characteristics and connectivity of pore structure. Sandstone from a coal mine in Tangshan City, Hebei Province, China, were selected for these tests. The sandstones were taupe, heterogeneous, and devoid of visible weak structures. The cylindrical samples were cored to the size of Φ50 × 100 mm approximately and oven dried for 24 h before application. Some other cylindrical rock samples were also prepared to determine the Young's modulus, Poisson's ratio, and uniaxial compression strength (UCS) of the rock materials according to the International Society for Rock Mechanics (ISRM) suggested method. surface of sample, and then the triaxial cell was dropped, and oils were injected into the cell. (2) Axial pressure (P1), confining pressure (P2), and hole pressure (P3, P4) were controlled by the computer to reach and maintain certain values. (3) The experiment was started with axial compression applied to the first experimental point, with P2, and P3 remaining constant at the lower side of the pore pressure P4, seepage results under pressure difference ΔP = P3-P4, and the ΔP-t curve was plotted. (4) Loading was increased to the next experimental point. Meanwhile, mercury intrusion porosimetry (MIP) was conducted for determining the porosity of rock specimens, and scanning electron microscopy (SEM) was carried out for observing the characteristics and connectivity of pore structure. Sandstone from a coal mine in Tangshan City, Hebei Province, China, were selected for these tests. The sandstones were taupe, heterogeneous, and devoid of visible weak structures. The cylindrical samples were cored to the size of Φ50 × 100 mm approximately and oven dried for 24 h before application. Some other cylindrical rock samples were also prepared to determine the Young's modulus, Poisson's ratio, and uniaxial compression strength (UCS) of the rock materials according to the International Society for Rock Mechanics (ISRM) suggested method.

Test results
Two sets of standard cylindrical samples were subjected to different confining and pore pressures. The parameters of these rocks obtained from the tests are presented in Table 1. Peak stress and strain increase couple with the increase in hydraulic gradient and confining pressure. The average UCS of the samples was 110.0 ± 10 MPa, and the permeability was approximately 0.2 × 10 −6 Darcy to 8 × 10 −6 Darcy. Figure 3 illustrates the SEM photograph of the sample F-3 and the curves of pore diameter versus the intrusion volume of F-1 and F-3. The fine-grained texture sandstone is mainly composed of quartz, calcite, iron oxide, and clay. The quartz crystal shape is principally granite-type single-crystal, the cementation and the contact method of particle are pore type and point-linear contact, respectively. A typical crystalline and blocky structure was shown in Figure 3a,b, the sandstone is a macroscopically heterogeneous fine grained rock material with an average bulk density of about 2530 kg/m 3 . Furthermore, the pore structure of the samples was tested using MIP, and the porosity was approximately 5% (Table 2), the vast majority of the pores were in the range of 0-1 µm (Figure 3c). Numerous studies demonstrate a correlation between pore geometry and strength of different type rocks [58][59][60]. Although specifically correlation for extrusive volcanic rocks, grain orientation may be an essential feature of sandstones and affect the strength anisotropy [59]. On this basis, we did not fully consider pore geometry and instead focused on grain anisotropy (heterogeneous) in this article.

Test Results
Two sets of standard cylindrical samples were subjected to different confining and pore pressures. The parameters of these rocks obtained from the tests are presented in Table 1. Peak stress and strain increase couple with the increase in hydraulic gradient and confining pressure. The average UCS of the samples was 110.0 ± 10 MPa, and the permeability was approximately 0.2 × 10 −6 Darcy to 8 × 10 −6 Darcy. Figure 3 illustrates the SEM photograph of the sample F-3 and the curves of pore diameter versus the intrusion volume of F-1 and F-3. The fine-grained texture sandstone is mainly composed of quartz, calcite, iron oxide, and clay. The quartz crystal shape is principally granite-type single-crystal, the cementation and the contact method of particle are pore type and point-linear contact, respectively. A typical crystalline and blocky structure was shown in Figure 3a,b, the sandstone is a macroscopically heterogeneous fine grained rock material with an average bulk density of about 2530 kg/m 3 . Furthermore, the pore structure of the samples was tested using MIP, and the porosity was approximately 5% (Table 2), the vast majority of the pores were in the range of 0-1 µm (Figure 3c). Numerous studies demonstrate a correlation between pore geometry and strength of different type rocks [58][59][60]. Although specifically correlation for extrusive volcanic rocks, grain orientation may be an essential feature of sandstones and affect the strength anisotropy [59]. On this basis, we did not fully consider pore geometry and instead focused on grain anisotropy (heterogeneous) in this article.  Logarithm of Pore Diameter (μm)   Figure 4 reveals the complete stress-strain-permeability curve of the typical sample during permeability testing. The two types of curves are stress-strain and permeability-strain. The permeability-strain curve initially decreases due to compression, gradually increases before the point of peak strength, and then increases sharply thereafter. Reference [8] showed three typical permeability-strain curves, but pointed out the peak value of permeability occur in the strain softening region with a frequency of 11/16. Our test results are consistent with these typical permeability-strain curves. Data from the laboratory observations can be used to calculate the mesosize parameters of stress-seepage coupled progress, especially in Equations (15) and (17), except for the methods mentioned in [51]. The material parameters of the numerical simulation in this article were selected accordingly ( Table 3). The parameters are equivalent mesosize parameters, which are different from normal laboratory tests [57].    Figure 4 reveals the complete stress-strain-permeability curve of the typical sample during permeability testing. The two types of curves are stress-strain and permeability-strain. The permeabilitystrain curve initially decreases due to compression, gradually increases before the point of peak strength, and then increases sharply thereafter. Reference [8] showed three typical permeability-strain curves, but pointed out the peak value of permeability occur in the strain softening region with a frequency of 11/16. Our test results are consistent with these typical permeability-strain curves. Data from the laboratory observations can be used to calculate the mesosize parameters of stress-seepage coupled progress, especially in Equations (15) and (17), except for the methods mentioned in [51]. The material parameters of the numerical simulation in this article were selected accordingly ( Table 3). The parameters are equivalent mesosize parameters, which are different from normal laboratory tests [57].

Numerical Simulation
The main experimental observations were further investigated by applying numerical simulations using flow-stress-damage failure process based on poroelasticity. The differential equations of the flow-stress-damage model are solved numerically by using the FEM with the software Comsol Multiphysics with Matlab. The adaptive LR was used to match the laboratory test when certain element failures occurred: the modulus and strength parameters reduced following Equation (7), the numerical model began a loop calculation, and all other conditions remained unchanged until no new failure occurred. Figure 5 shows the modeling setup, discrete mesh, and element quality histogram of the 3D numerical model in this article. The multifrontal massivelyparallel sparse-direct solver was selected, and the time-stepping method involved backward differentiation formulas (BDF) with the order of accuracy varying from one (that is, backward Euler) to five [61]. BDF methods have been used for a long time and are known for their stability. Analytical solutions were implemented to reveal the indoor tests to a certain extent. Coupled 3D model was applied for experiment simulation to study stresses-seepage and the failure in triaxial stress. The input mechanics parameters and geometry structure for this model were rely on properties of test specimens and an experimental setup. The results are discussed, and additional details are presented and compared with the laboratory observations. According to the cylindrical rock samples, the numerical model with the same size of Φ50 × 100 mm was used. For the stress boundary, the bottom boundary was fixed, whereas the top boundary was loaded at 0.005 mm/s with a confining pressure of σ2 = σ3 = 6 MPa. True triaxial conditions (σ2 ≠ σ3) were not considered in this study because of space limitations. Such conditions will be studied in future publications. For the

Numerical Simulation
The main experimental observations were further investigated by applying numerical simulations using flow-stress-damage failure process based on poroelasticity. The differential equations of the flow-stress-damage model are solved numerically by using the FEM with the software Comsol Multiphysics with Matlab. The adaptive LR was used to match the laboratory test when certain element failures occurred: the modulus and strength parameters reduced following Equation (7), the numerical model began a loop calculation, and all other conditions remained unchanged until no new failure occurred. Figure 5 shows the modeling setup, discrete mesh, and element quality histogram of the 3D numerical model in this article. The multifrontal massively-parallel sparse-direct solver was selected, and the time-stepping method involved backward differentiation formulas (BDF) with the order of accuracy varying from one (that is, backward Euler) to five [61]. BDF methods have been used for a long time and are known for their stability.

Numerical Simulation
The main experimental observations were further investigated by applying numerical simulations using flow-stress-damage failure process based on poroelasticity. The differential equations of the flow-stress-damage model are solved numerically by using the FEM with the software Comsol Multiphysics with Matlab. The adaptive LR was used to match the laboratory test when certain element failures occurred: the modulus and strength parameters reduced following Equation (7), the numerical model began a loop calculation, and all other conditions remained unchanged until no new failure occurred. Figure 5 shows the modeling setup, discrete mesh, and element quality histogram of the 3D numerical model in this article. The multifrontal massivelyparallel sparse-direct solver was selected, and the time-stepping method involved backward differentiation formulas (BDF) with the order of accuracy varying from one (that is, backward Euler) to five [61]. BDF methods have been used for a long time and are known for their stability. Analytical solutions were implemented to reveal the indoor tests to a certain extent. Coupled 3D model was applied for experiment simulation to study stresses-seepage and the failure in triaxial stress. The input mechanics parameters and geometry structure for this model were rely on properties of test specimens and an experimental setup. The results are discussed, and additional details are presented and compared with the laboratory observations. According to the cylindrical rock samples, the numerical model with the same size of Φ50 × 100 mm was used. For the stress boundary, the bottom boundary was fixed, whereas the top boundary was loaded at 0.005 mm/s with a confining pressure of σ2 = σ3 = 6 MPa. True triaxial conditions (σ2 ≠ σ3) were not considered in this study because of space limitations. Such conditions will be studied in future publications. For the Analytical solutions were implemented to reveal the indoor tests to a certain extent. Coupled 3D model was applied for experiment simulation to study stresses-seepage and the failure in triaxial stress. The input mechanics parameters and geometry structure for this model were rely on properties of test specimens and an experimental setup. The results are discussed, and additional details are presented and compared with the laboratory observations. According to the cylindrical rock samples, the numerical model with the same size of Φ50 × 100 mm was used. For the stress boundary, the bottom boundary was fixed, whereas the top boundary was loaded at 0.005 mm/s with a confining pressure of σ 2 = σ 3 = 6 MPa. True triaxial conditions (σ 2 = σ 3 ) were not considered in this study Energies 2018, 11,1923 10 of 16 because of space limitations. Such conditions will be studied in future publications. For the flow boundary, the bottom boundary showed a high pressure at p down = 3.1 MPa, whereas the top boundary presented a low pressure at p up = 1.3 MPa, and flow did not occur in the other boundary. Figure 6a,b show the numerically simulated progressive failure process and flow field distribution under confining and hydraulic pressures, respectively, with heterogeneity degrees of m = 6. First, failure initiation is at a low strength element because of the distribution of spatially statistical material parameters. Tensile and shear damages (0 < D ≤ 1 for shear damage and −1 ≤ D < 0 for tensile damage) occur during this process [57]. Figure 6a Figure 6a,b show the numerically simulated progressive failure process and flow field distribution under confining and hydraulic pressures, respectively, with heterogeneity degrees of m = 6. First, failure initiation is at a low strength element because of the distribution of spatially statistical material parameters. Tensile and shear damages (0 < D ≤ 1 for shear damage and −1 ≤ D < 0 for tensile damage ) occur during this process [57]. Figure 6a  At Step_143_13, shear damage propagates and results in this shear damage zone extends (at Step_143_16) into the nearby area, and complete failure occurs at Step_143_19. Meanwhile, shear damage and the tensile damage appear alternately. Shear damage mainly appears at the beginning, whereas tension damage mainly appears and acts on the rock's typical properties to induce complete failure. Increasing shear damage appears with small m (Figure 6c), and additional tension damage occurs with large m even with initial failure emerging on the boundaries (the area marked with A0 in Figure 6d). The discreteness exceeds that with small m, and the strength exhibits the inverse trend accordingly. This outcome means that fewer elements of the sample exist near the mathematical expectation of the strength with small m, thereby resulting in lower strength. However, no exact relations between m and strength exist. A certain study supported the determination of the relations under some conditions [62]. Complete stress-strain curves with/without water were obtained via our 3D numerical mode (Figure 7a). The results show brittle/ductile characteristics with the At Step_143_13, shear damage propagates and results in this shear damage zone extends (at Step_143_16) into the nearby area, and complete failure occurs at Step_143_19. Meanwhile, shear damage and the tensile damage appear alternately. Shear damage mainly appears at the beginning, whereas tension damage mainly appears and acts on the rock's typical properties to induce complete failure. Increasing shear damage appears with small m (Figure 6c), and additional tension damage occurs with large m even with initial failure emerging on the boundaries (the area marked with A0 in Figure 6d). The discreteness exceeds that with small m, and the strength exhibits the inverse trend accordingly. This outcome means that fewer elements of the sample exist near the mathematical expectation of the strength with small m, thereby resulting in lower strength. However, no exact relations between m and strength exist. A certain study supported the determination of the relations Energies 2018, 11, 1923 11 of 16 under some conditions [62]. Complete stress-strain curves with/without water were obtained via our 3D numerical mode (Figure 7a). The results show brittle/ductile characteristics with the decrease/increase in heterogeneity degree. The stress and strain of dry samples are, notably, slightly higher than those of the saturated and wet samples. The results explain well the weakening effect of the rock strength by water. Interestingly, stress-strain curves enter the unstable fracture stage in advance with a small m. Rock brittleness gradually decreases and the nonlinear mechanical characteristics are gradually enhanced. Flow field distribution is more evident in the new cracks. When an element of rock sample undergoes compression or failure, its permeability changes as the relationship in Equation 15 stated above. Hydraulic properties can be changed with the stress state. Meanwhile, new group of fractures form a macroscopical seepage pathway (area marked with A, B, and C in Figure 6b), and the formed seepage pathway throughout the sandstone sample at complete failure. Moreover, the same situation occurs under confining and hydraulic pressures with different heterogeneity degrees. The trend also indicates that the rock sample will immediately reach the last failure, which occurs in a short time. This phenomenon may aid in warning and forecasting in mining and tuning engineering; moreover, it can help in the analysis of the premonition information, especially in hydraulic fracturing. The pressure for crack initiation in hydraulic fracturing is approximately 10 MPa [1], which means increased tension damage will appear in advance along with a varying failure pattern. Before failure, permeability and porosity decrease with the increase in stress. However, the trend may be reversed, and permeability maybe increase by more than one to three orders of magnitude. Several experimental tests have illustrated that a strong correlation between permeability and volumetric change exists.

Effect of Different Parameters and Strength Criteria
Rock failure is a gradual and non-continuous process that shows a localization phenomenon. Rock/rock mass possesses residual strength after failure. RSC illustrates the value of residual strength compared with initial strength. The stress-strain curve may be different when applied to different types of rock. Figure 7b illustrates the rules of different RSC values. The stage of the curve after peak Flow field distribution is more evident in the new cracks. When an element of rock sample undergoes compression or failure, its permeability changes as the relationship in Equation (15) stated above. Hydraulic properties can be changed with the stress state. Meanwhile, new group of fractures form a macroscopical seepage pathway (area marked with A, B, and C in Figure 6b), and the formed seepage pathway throughout the sandstone sample at complete failure. Moreover, the same situation occurs under confining and hydraulic pressures with different heterogeneity degrees. The trend also indicates that the rock sample will immediately reach the last failure, which occurs in a short time. This phenomenon may aid in warning and forecasting in mining and tuning engineering; moreover, it can help in the analysis of the premonition information, especially in hydraulic fracturing. The pressure for crack initiation in hydraulic fracturing is approximately 10 MPa [1], which means increased tension damage will appear in advance along with a varying failure pattern. Before failure, permeability and porosity decrease with the increase in stress. However, the trend may be reversed, and permeability maybe increase by more than one to three orders of magnitude. Several experimental tests have illustrated that a strong correlation between permeability and volumetric change exists.

Effect of Different Parameters and Strength Criteria
Rock failure is a gradual and non-continuous process that shows a localization phenomenon. Rock/rock mass possesses residual strength after failure. RSC illustrates the value of residual strength compared with initial strength. The stress-strain curve may be different when applied to different types of rock. Figure 7b illustrates the rules of different RSC values. The stage of the curve after peak strength changed with different RSC values. RSC minimally affected the elastic modulus during the initiation of loading on account of the accumulated damage process. Then, the trend changed near the peak strength point. However, apparent stress drop did not occur with higher RSC values. Inversely, brittleness was revealed under small values of RSC. Furthermore, the failure pattern changed with different RSC values and can be used to explain the marked shear smooth crack under certain conditions.
The ISRM suggests a loading rate of 0.002-0.005 mm/step in rock laboratory tests. The results may be changed under different LR. A large LR induced impact failure, whereas the low LR showed rheological behavior. Both conditions are needed to specialize rock mechanics. Thus, this factor does not merit any additional discussion in this work. Several numerical tests were conducted under different constant and adaptive LRs for comparison with the laboratory test and to verify the perfect accuracy of the numerical simulations. The adaptive step is more applicable in controlling the computation time effectively than is the constant LR. Moreover, the adaptive step changed from 0.00002 mm/step to 0.02 mm/step automatically with the development of damage. Figure 7c shows the comparison between constant and adaptive LR, which yielded virtually identical types of stress-strain curves. However, slight differences exist on the stress or strain level. Meanwhile, the results of the adaptive LR and LR = 0.005 mm/step almost overlap. In general, the constant LR can match the accuracy requirement, but should be studied further under certain specific conditions. Such conditions include an exceedingly small constant LR for inducing new damage and massive amounts of computation on 3D real damage modeling, which entails long computation time and large resources.
The widely used strength criteria are MC, DP, TS, and Hoek-Brown criteria. Figure 7d illustrates the complete stress-strain curve by using different strength criteria. The results from the TS criterion are the most conservative, and the stress-strain level under the DP criterion is enhanced and the highest. As stated above, the production and dissipation of excess pore water pressure causes damage development. When excess pore pressure appears and is sufficient to induce new damage, the failure criterion is implemented and used to estimate the failure of the element. Strictly speaking, the Hoek-Brown criterion is an empirical function that is commonly used while dealing with rock masses of variable quality because the material parameters can be estimated based on simple field observations together with knowledge of the UCS of the intact rock. This article did not cover the topic given the uncertainty regarding experiential parameters.  Figure 8c,d show the complete permeability-strain curve by using different combined variables. By keeping other variables constant and changing the damage coefficient α D , the stress sensitive factor β, and RSC, the influence of the three items on the permeability-strain can be determined. Figure 8c illustrates that β = 0.3 MPa −1 , α D = 3, and RSC = 0.03 can match well with sample F-3 in the laboratory test. α D controls the level of the ultimate permeability, and the stress sensitive factor β controls the entire trend and influences the contribution of stress to permeability. Figure 8d reveals the results under different α D and RSC values; RSC and α D play opposite roles in influencing the level of the last permeability, except for the experimental error and rock sample properties. Thus, the three items affect the permeability-strain relation. Although the permeability notably increased by more than three to four orders of magnitude in actual mining, we obtained results of approximately one to two orders of magnitude, which clearly underestimates the values obtained from the laboratory tests and needs further study. The combined variables should be tested and selected reasonably in actual engineering. Notably, no significant difference in comparison with the laboratory results was observed. Numerical results on the permeability-strain curve show the same change trend without regard to downsized flow channels, and other channels are completely blocked. The results perfectly explained the influence on permeability-strain under different combined variables. In addition, the stress sensitive factor may be changed in accordance with failure initiation, propagation, coalescence, and out of control. Thus, further specific and in-depth research should be conducted. (c) (d)

Conclusions
A 3D flow-stress-damage method based on poroelasticity was applied to investigate the failure mechanism of rock samples under hydraulic conditions; the method was then compared with the laboratory observations and was analyzed to determine the influence of different material parameters and strength criteria. The study found that the stress sensitive factor β and the damage coefficient αD are key variables of permeability evolution-characterization under pre-peak and post-peak, separately. In addition, RSC influences post-peak permeability change. Based on the test data, a twostage permeability evolution model was proposed and verified, and the model can be used to illustrate the failure mechanics under hydraulic conditions and to match different rock types for considering heterogeneity.
Different heterogeneity degrees affect the mechanics beyond the elastic limit strain. The strength decreased with the decrease in heterogeneous degree. The complete stress-strain curve shows ductility at small m values and brittleness at large m values. Moreover, a large RSC means evident plasticity. Inversely, the brittleness is revealed under smaller values of RSC. In the present numerical investigations, all parameters mentioned above exert important influence on the stress-strainpermeability relation. Meanwhile, whether the selected strength criterion is conservative should depend on the different rock types in actual engineering, given that hydraulic properties are beneficial in determining the appropriate parameters for a certain rock.
Author Contributions: S.C. and C.W. conceived this model; S.C. performed the simulations and wrote the paper; C.W. and H.L. analyzed the numerical simulation and test data; T.Y. and W.Z. contributed analysis tools; P.G.R. revised the paper and gave final approval of the version to be submitted.
Funding: This study was supported by the National Natural Science Foundation of China (grant numbers

Conclusions
A 3D flow-stress-damage method based on poroelasticity was applied to investigate the failure mechanism of rock samples under hydraulic conditions; the method was then compared with the laboratory observations and was analyzed to determine the influence of different material parameters and strength criteria. The study found that the stress sensitive factor β and the damage coefficient α D are key variables of permeability evolution-characterization under pre-peak and post-peak, separately. In addition, RSC influences post-peak permeability change. Based on the test data, a two-stage permeability evolution model was proposed and verified, and the model can be used to illustrate the failure mechanics under hydraulic conditions and to match different rock types for considering heterogeneity.
Different heterogeneity degrees affect the mechanics beyond the elastic limit strain. The strength decreased with the decrease in heterogeneous degree. The complete stress-strain curve shows ductility at small m values and brittleness at large m values. Moreover, a large RSC means evident plasticity. Inversely, the brittleness is revealed under smaller values of RSC. In the present numerical investigations, all parameters mentioned above exert important influence on the stress-strain-permeability relation. Meanwhile, whether the selected strength criterion is conservative should depend on the different rock types in actual engineering, given that hydraulic properties are beneficial in determining the appropriate parameters for a certain rock.