Investigation on Coupled Fluid-Flow and Stress in Dual Model Rock Mass with Time-Dependent Effect and Its Simulation

The coupled fluid-flow and stress model with time-dependent effect is developed. In this model, rock mass is simulated as a homogeneous and isotropic dual-porosity, dual-permeability continuum. Darcy’s law is used to describe the fluid flow in a porous medium, and cubic law is used to describe the flow in fractures. The finite element method is introduced to solve the system, and the effects of fractures and fracture spacing are considered numerically. The numerical results indicate that fractures have a significant impact on the rocks’ displacement and pore pressure. It also shows an increase in plastic strain with decreasing fracture spacing, as does the creep strain. The coupled process is highly sensitive to the fracture spacing. A series of numerical simulations are conducted to better understand the complex coupled processes, which leads to an improved understanding of different aspects of naturally fractured reservoirs and may impact on experimental design to explore these attributes in a real reservoir situation.


Introduction
The analysis of those problems related to the fractures in rock mass is intrinsically difficult due to their complex characteristics.All of these solid rock engineering problems have two distinct porosities, one in the matrix and the other in the fractures.Because of the different storage and conductivity characteristics of the matrix and fractures, these reservoirs can be represented by dual model reservoirs.Dual-porosity, single permeability model and dual-porosity, dual-permeability model can be used to describe the fluid solid coupled effect for not only the weakly compressible fluid but also the slightly compressible fluid with adsorption mechanism, such as water flowing though fractures in deforming aquifers [1][2][3], naturally fractured reservoirs [4,5], and other mechanical influences.Furthermore, a number of papers have extensively addressed and investigated in the CO 2 sequestration in deep coal seams and coalbed methane recovery with dual porosity or dual permeability model [6][7][8][9][10][11][12][13].Wu et al. used finite element model to quantify the net change in permeability, the gas flow, and the resultant deformation in a prototypical coal seam with a dual poroelastic model [6].Wei and Zhang [7,8] used finite volume method and fully implicit finite-difference method fully implicit finite-difference method to simulate the coalbed methane (CBM) recovery with a triple-porosity/dual-permeability model, respectively.Masoudian and co-researchers [9][10][11]13] did a lot of work to study the effect of CO 2 injection on reservoir and geomechanical performance of the coal seam in a more comprehensive way.
Ever since Barenblatt et al. firstly established the dual-porosity model by introducing the theory of mixtures [14], there have been numerous studies reporting the theory and numerical analysis on modeling fractured formations with the dual-porosity, dual-permeability approach [5,[15][16][17][18].Some researchers followed the conventional fluid flow modeling coupled with geomechanics through stress dependent rock properties and during the development process [19][20][21].Considering pore deformation, Valliappan and Khalili-Naghadeh [20] derived a set of coupled differential equations governing the behaviour of fissured porous media.Chen and Teufel [19] focused on theoretical developments with emphasis on the identification of the critical coupling and physical interpretations of the parameters involved.Taking advantage of the joint-mechanics theory, Bagheri and Settari [21] developed general, rigorous coupling between the fluid-flow equation and deformation of fractured media.The need for coupled modeling of hydro-mechanical model has been addressed and investigated in various studies whether the analytical and numerical modelling [22][23][24][25].For example, Zhang et al. [22] used the finite element method to investigate the effects of fracture spaces on the displacement, stress and pressure in dual-porosity media.However, they neglected the matrix-fracture interaction transfer function.Rutqvist [23] proposed a linked multicontinuum and crack tensor approach for modelling of coupled geomechanics, fluid flow and transport in fractured rock.Hu et al. [25] used a numerical manifold method for analyzing fully coupled hydro-mechanical processes in porous masses with discrete fractures.Some even considered the chemo-hydro-mechanical model for the coupled multiphysics of carbon dioxide sequestration [10,11].
Although all of these works provide better extension to the existence theory and models, the long-term creep effect on the matrix with coupled fields is seldom considered.In practice, the time-dependent effect has a great impact on the characteristics of rock material and the stability of rock engineering, such as deep underground rock engineering, high slope engineering, dam base engineering and nuclear waste disposition engineering [26][27][28][29][30].This influence has been investigated through experimental and theoretical [27,28,31,32] and numerical modeling [33].Some researchers considered the influence of pore pressure on the creep mechanical behavior [34,35].However, these studies have improved our understanding on the time-dependent influence on rock mass, we still should know its influence on dual-porosity, dual-permeability system with important coupling between mechanical and fluid-flow transport behaviors.
Therefore, the focus of this paper is to investigate the time-dependent response to coupled hydro-mechanical with dual-porosity, dual-permeability medium.In order to achieve this target, dual-porosity, dual-permeability model is applied to illustrate the system with coupling between mechanical and fluid-flow transport both in fractures and in pores.The physical and mathematical models of dual-porosity, dual-permeability medium are used to derive the governing equations of rock mass deformation and coupled fluid flow.The time-dependent influence on the deformation is also considered.From the numerical simulation, four different cases are used to describe the difference with consideration of the effects of fractures and the spacing of the fractures.The numerical results are performed in a comprehensive way including displacement contour, displacement vector, pore pressure contour, equivalent plastic strain contour, and equivalent creep strain contour.

Governing Equations
The research on dual-porosity, dual-permeability model is transformed from the homogeneous isotropic medium to an orthotropic medium, an even more complicated medium.The formulations have two properties: the first is the porous rock storage capacity is high and transmissivity is low; fracture storage capacity is low and fracture transmissivity is relatively high.The second is the possibility that exists on flow through both the blocks and the fractures with a transfer function describing the fluid exchange between the two continua.
In the following, based on the references [20,22], a set of field equations for rock time-dependent deformation and fluid flow are defined.These equations are coupled with the interaction effect of pore pressure and rock deformation in rock matrix and fractures.In addition, these equations are based on the following assumptions: (1) Rock mass is a dual-porosity, dual-permeability homogeneous isotropic continuum.
(3) Rock material has the time-dependent mechanical behavior.(4) Darcy's law is used to describe fluid flow in porous medium, and cubic law is used to describe the flow in fractures.(5) The density of fluid in rock mass is constant.(6) Isothermal condition is used.

Rock Mass Deformation
There are two distinct types of porosity and permeability coexisting in rock volume.The traditional effective stress form is used with single fluid pressure in the system.Here, by extending Terzaghi's effective stress saturated rock, and considering the effects of pore pressure in matrix and fractures separately, the total stress in a general form is expressed as [15], where σ ij is total stress tensor, a 1 , a 2 and p 1 , p 2 are the Biot coefficient of effective stress and fluid pressure of blocks and fractures, respectively, and δ ij is the unit tensor.Equilibrium equation is where b i is the body force vector.The stress-strain relationship of homogeneous isotropic continuum is where λ and G are lame constants and shear modulus, respectively.
The strain-displacement relationship is defined as The total strain can be assumed to be as follows: where ε e ij is elastic strain, ε p ij is viscous strain, and ε v ij is viscous (creep) strain.The above expression should be determined by the actually rheology equation [36][37][38].Semenov et al. [36] discussed the classical and generalized rheological models for the materials in civil and mechanical engineering.Masad et al. [37] used Perzyna's viscoplasticity along with the Druker-Prager yield surface for modeling hot mix asphalt.Thus, the above Equation ( 5) is only one of the assumptions.In Section 3, the Druker-Prager model is used to exhibit the pressure-dependent yield and the time hardening form of the power law creep model for viscous behavior.
Substituting the Equations ( 3) and (4) into Equation (2), we yield the stress field governing equation of stress-seepage model for the dual-porosity medium Using the principle of λ = 2Gv/(1 − 2v) (here, v is the Poisson's ratio), Equation ( 6) can be rewritten in the form

Fluid Flow Governing Equation
Assuming that the fluid flow through porous media can be described by the Darcy's law, if the gravity effect is considered, then the Darcy volumetric flow rate can be defined as the following equation: where v is the fluid's relative velocity in porous media, ρ f and µ is the fluid's density and viscosity, respectively, k 1 is the permeability of the porous media, p 1 is the fluid pressure on porous media, g = g∇z denotes the gravitational acceleration, and z is the vertical coordinate.Here, the permeability variation in the matrix is not considered.For the sake of simplicity, we assumed that the fracture wall is smooth and the flow is in the laminar range.Here, we use the cubic law to describe the tangential flow through a fracture, which is expressed by [39] where v 2 is the fluid's average velocity in fracture, d is the average fracture width, and p 2 is the fluid pressure on fracture media.The permeability in the fracture is defined as [39].
The relative fluid velocity can be defined as where v ωi , φ ω , v ωi f , v is (ω = 1 standing for the porous medium and ω = 2 for the fracture) represent relative fluid velocity, porosity, absolute fluid velocity and velocity of the blocks, respectively.The solid mass continuous equation is defined as where φ = φ 1 + φ 2 is the total porosity of rock mass, and ρ s is the solid density.Due to v is ∂/∂x i << ∂/∂t, we can get Considering the fluid exchange between the porous media and fractures, the fluid mass continuous equation can be written as where is the interporosity flow coefficient, which is related to the permeability and geometric characteristic of the matrix) is defined as the matrix-fracture interaction function.The interaction function has many different equations and has been investigated by many authors [40][41][42].The state equation of fluid and solid blocks can be expressed as where c f , c s are the compressibilities of fluid and solid blocks, and p = a 1 p 1 + a 2 p 2 is the total pressure: where ω = 1, 2, a ω is the Biot coefficient shown in Equation ( 1) of a ω 1 , a ω 2 , which are the constants related to the coefficient of compressibility and porosity [20].The above equation is the governing equation for fluid flow in dual-porosity medium.

Boundary Conditions
Generally, there are three kinds of boundary conditions that are the displacement boundary condition, stress boundary condition and seepage boundary conditions.
The displacement boundary initial condition is The stress boundary condition can be written as where l, m, n are the three exterior normal direction cosine to the boundary, σ i , τ ij is the component of stress tensor, and Fx , Fy , Fz are the pressure forces.This boundary condition can be rewritten as where L is the coefficient matrix of the above set of equations, and F is the vector of pressure force.
There are three kinds of seepage boundary conditions.If we already know the fluid pressure at each point at any time, then we have The second seepage condition is used to describe the fluid flow velocity where v n is the component of v normal to the boundary for the impervious boundary v = 0. Cauchy boundary condition is both the potential and its normal derivative prescribed on the boundary.Thus, p verifies the following equation: where α, β are Cauchy data.

Model Description and Input Parameters
Using Galerkin's weighted residual method in each element and finally assembling into the whole models, we can get the global matrix of the system and the increment schemes.In order to study fractures and the time-dependent effects on the system, a two-dimensional model [22], which has a length of 12 m and a height of 10 m, is constructed using ABAQUS (Dassault Systèmes, Vélizy-Villacoublay, France) as shown in Figure 1.This model is used to describe the rock mass deformation in a saturated dual-porosity, dual-permeability model under the external load.A concentrated force P is applied as shown on the middle part of the top surface.The bottom of the model is rollered and the displacement constraints on the x-direction to the left and right sides.At the top, an outlet boundary condition is specified with the fluid pressure set to zero, and the pressures on the left and right boundaries are γ w h.Here, γ w is the unit weight of water and h is the hydraulic height.A no-flow boundary condition is applied to the bottom of the system.Two fracture groups in horizontal and vertical direction are contained in this model with the same fracture aperture.The same four-node bilinear displacement and pore pressure element are used in matrix and fractures ignoring the cracking effect on fractures in case 1.In cases 2-4, four-node bilinear displacement and pore pressure element in matrix and cohesive element type in fractures are used to simulate the stress and strain effect in a dual medium system.Different fracture spacing b sp with 2 m, 1 m and 0.5 m and the fracture aperture b = 0.01 m are set in cases 2-4.Tangential flow is the fluid flow in the length direction of the fracture and normal flow refers to the fluid leak-off into the surrounding medium.Fluid leakage coefficient defines a pressure flow leak-off relationship like the coefficient a d ρ f /µ of the matrix-fracture interaction function.Since the tangential flow is absent from the porous medium, the tangential flow will enhance the crack at the fractures if the tangential permeability is defined.The tangential flow and the normal flow are only existed in cases 2-4.The values of the properties are listed in Table 1 and the simulated conditions are shown in Table 2.

Results and Discussion
The displacement contour, displacement vector diagraph and pore pressure contour of the four cases simulation with the rock consolidation under the top pressure are shown in Figures 2-4.From Figures 2-4, the difference in the displacement contour, displacement vector and pore pressure under different conditions can be seen.The continuity being weakened from case 1 to cases 2-4 can be demonstrated because of the existence of fractures.The cracking effect of tangential fluid flow plays a dominant role in enhancing the discontinuity with the decrease of fracture spacing.In addition, notice that the displacement increases with the decrease of fracture spacing.This implies that the fractures decrease the stiffness of the system and increase the flexibility of the rock mass.Figure 3 shows the displacement of the points around the fractures is greater than that of other points, which means that the crack is stress concentration and more dangerous in the whole system, which has been proved by many experiments.However, case 1 in Figure 3 has no evidence of cracking in addition to the obvious downward displacement of the middle section.In fact, with the uploading and the effect of fluid, the displacement becomes larger when the fracture spacing is taken into account.
Figure 4 shows the different situations of pore pressure diagram of single medium case 1 and dual-porosity cases 2-4.Because the permeability of fracture is larger than that of the pore by several orders of magnitude, most of the water flows out of the cracks.Due to the middle part of the applied load, the fractures in this part are more easily cracked.As the cracks increase, the volume of water discharged from the cracks will be more.The vertical displacement with respect to time curves of the center point on the top surface is shown in Figure 5. Specifically, the displacement of the center point at the end of consolidation with time 1 × 10 9 s from case 1 to case 4 is 9.55 × 10 −3 m, 9.60 × 10 −3 m, 2.14 × 10 −2 m, 1.30 × 10 −2 m, respectively.The figure shows that the displacement considering the fracture effect in case 2 to case 4 is larger than that in case 1.The figure also shows the displacement increases as the fracture spacing decreases, indicating the significance of considering the fracture effect in the whole system.Similar trends are observed from Figure 6, where the displacement over the depth of the center line is displayed in different cases.The Drucker-Prager model with associated flow rule is utilized to describe the plasticity behavior of the matrix.The Drucker-Prager yield function is written as [43] F s = q − ptanβ − d = 0, where β is the material's angle of friction d is its cohesion, p = −1/3σ ii is the equivalent pressure stress and q is the Mises equivalent stress.Here, we set the friction angle 35 • and dilation angle 0 • in the above four cases .The plasticity strain contours are shown in Figure 7. Figure 8 compares the creep strain contours by using the time hardening form of the power law creep model with the creep material parameters A = 2.4 × 10 −8 , n = 0.15, m = −0.8,t = 1 × 10 6 .The time hardening form of the power law model is εcr = A( σcr ) n t m , where εcr is the equivalent creep strain rate, which is defined as εcr = | εcr 11 |; σcr is the equivalent creep stress; t is the total or the creep time; and A, n, and m are user-defined creep material parameters.It is observed from Figure 7 that the presence of cracks weakens the load-bearing capacity and enlarge the plastic area.As is shown in this figure, the decrease in the length of fracture spacing has a significant impact on the increase of plastic area.Therefore, the ignorance of the crack effect will lead to a wrong evaluation during the calculation of elastic-plastic materials.Similar trends are observed from Figure 8 for creep strain on the four cases, which increases as the fracture spacing decreases, with very significant variations by the creep strain contours.These observations indicate that neglecting the crack effect can lead to significant errors.We still should care about the creep deformation even when considering the crack effect.In order to illustrate the fracture effect deeply, we give the creep strain of the monitoring point at node 442 with x = 8.01 m, y = 9.00 m in case 1 and case 3 as shown in Figure 9. Figures 8 and 9 indicate that the creep strain of considering the crack effect in case 3 is greater than that of ignoring the crack effect in case 1.The above numerical simulations show that whether fracture spacing or time-dependent effect factor should be considered to assess the safety of practical engineering.

Conclusions
Numerical investigations have been conducted to enhance the understanding of the coupled fluid-flow and geomechanics processes with a time-dependent effect.Four cases are carried out to study the effects of fractures and fracture spacing on the single medium model and the double medium model.
It is shown that fractures and fracture spacing have a significant impact on the rock mass system.Compared to a double medium model, a single medium model has a more continuous change, which can be seen from the displacement contours and pore pressure contours.In addition, the displacement of double medium model is greater than that of a single medium model.The physical mechanism is that fractures reduce the stiffness of rock and increase the flexibility of rock mass.The concentrated force and tangential flow enhance the crack of fractures.Therefore, the permeability coefficient of fractures is several orders of magnitude larger than that of pores.As a result, the majority of the fluid flow through the fractures, and the flow of fluid in the pores is relatively slow.
Fracture spacing is also a critical factor in the rock mass modeling because the plastic area and creep strain increase when the fracture spacing decreases.The maximum creep deformation of the double medium is greater than that of the single medium model, and the maximum creep deformation increases with decreasing fracture spacing.The conclusion that can be drawn is that the factors of fracture distribution and time-dependent effect should not be ignored in a realistic rock mass modeling.
This work can help to provide a further insight in the time-dependent effect on coupled fluid-flow and stress in dual medium.Furthermore, these results have significant implications for the characterization of natural fracture reservoirs, and may impact experimental design to understand such attributes in a real reservoir situation, especially in rocks with distinct rheological features.
the cracking effect in fracture, b ap (m) 1 Case 4 Fracture spacing considering the cracking effect in fracture, b ap (m) 0.5

Figure 5 .Figure 6 .
Figure 5. Displacement time curves at the center point of different cases.

Figure 9 .
Figure 9. Displacement time curves at node 442 of case 1 and case 3.

Table 1 .
Property parameters of simulation model.

Table 2 .
Simulation models under different conditions.