Numerical simulation of fluid‐solid coupling of fractured rock mass considering changes in fracture stiffness

The fluid‐solid coupling simulation of fractured coal and rock is helpful to study the influence of the fracture structure on the stress sensitivity of the permeability. Therefore, a key factor in improving the reliability of such a numerical simulation is selecting parameters that match measured laboratory results. This paper firstly illustrates the use of a common inversion method to select the fluid‐solid coupling parameters through the constant joint stiffness method (CJSM). The main principle, application range, and the shortcomings of this method (CJSM) are analyzed: It can be used to simulate the measured results in a lower stress range. However, the stress sensitivity of permeability is zero when the stress reaches the critical stress (the stress where the hydraulic aperture reaches the residual hydraulic aperture). In order to solve this problem, a numerical simulation method for fluid‐solid coupling using the varying joint stiffness method (VJSM) is proposed based on the characteristic that the compressibility of a fracture decreases with increasing stress. In this method, the compressibility of the hydraulic aperture can be changed by increasing the joint stiffness during the loading process. Thus, the hydraulic aperture can continuously change and there is no need to set the residual hydraulic aperture. Finally, fluid‐solid coupling simulations under different stress states were performed using the VJSM in this study. The numerical simulation results were basically consistent with the experimental results, which indicated that a fluid‐solid coupling simulation using the VJSM can match experimental results very well and can be further used to study the influence of the persistent fracture angle on the permeability.

are attracting more attention in the field of the fractured rock mass permeability. [8][9][10][11][12] The main purpose of using a numerical methodology for fractured rocks is to create a realistic model for fractures with a complex geometry, as well as the interaction between joints and rock blocks. The hybrid numerical discrete fracture network-distinct element method (DFN-DEM) and bonded particle (or grain) model-distinct element method (BPM-DEM) provide effective approaches for calculating the mechanical properties of fractured rocks. [13][14][15][16][17][18] A nonlinear constitutive model of fractures was employed by Yao et al, 19 in which the nonlinear normal stress deformation relationship, tangential slipping, and dilation effects were considered. It well simulated the fracture stiffness changes with an increase in the normal stress, but set the joint stiffness by subsection. Six patterns of fracture sets were generated and studied by Ren et al 20 to show the effects of the fracture distribution on the hydro-geometric anisotropy and permeability anisotropy. The numerical simulation of solute transport in fractured rocks carried out by Zhao et al 21 showed that stress not only significantly changes the solute residence time through the fracture networks, but also changes the solute travel paths. The fracture roughness characteristics also affect the seepage characteristics in the fracture rock found in a numerical simulation. 21,22 However, because of the limitations on numerical simulation, such as simplification and assumptions for the numerical model, many influencing factors cannot be considered in a numerical simulation. Thus, under many conditions, numerical modeling in conjunction with numerous laboratory experiments is efficiently used to assess the impact of mining on a regional scale. The fluid-solid coupling simulation of fractured coal and rock is helpful to study the influence of the fracture structure on the stress sensitivity of the permeability. Therefore, the key factor to improve the reliability of numerical simulation is selecting parameters for the numerical simulation that match the measured results. At present, most of the commonly used methods for selecting fluid-solid coupling parameters adopt CJSM. The matching effect of the experimental results is relatively acceptable when the stress changes little. However, for some cases when the stress variation is relatively large, the matching effect of this method is substantial decline. This paper firstly illustrates an inversion method to select the fluid-solid coupling parameters using CJSM. The main principle, application range and the shortcomings of this method (CJSM) are analyzed. Based on this, a numerical simulation method for fluid-solid coupling using the varying joint stiffness method (VJSM) is proposed based on the characteristic that the compressibility of a fracture decreases with increasing stress. Fluid-solid coupling simulations under different stress states were performed using the VJSM in this study.

OF CONSTANT JOINT STIFFNESS METHOD
Universal Distinct Element Code (UDEC), a 2D numerical program based on the distinct element method (DEM) for discontinuum modeling, is used to simulate the response of a jointed rock mass subjected to static loading and fluid flow (Itasca Consulting Group Inc, 2011). A stress-fracture-flow coupling simulation can be carried out using Equations (1)-(3), and the entire flow process is shown in Figure 1.
where a is the contact hydraulic aperture, m; μ is the fluid viscosity, Pa s; ΔP is the pressure difference between the two adjacent domains, Pa; l is the length assigned to the contact between the domains, m; a 0 is the contact hydraulic aperture at zero normal stress, m; σ n is the normal stress at the contact, Pa; k n is the contact normal interface stiffness, Pa/m; P is the new pore pressure of the domain, Pa; P 0 is the original pore pressure, Pa; K w is the bulk modulus of the fluid, Pa; Q is the sum of the flow rates of the regional joint, m 2 /s; and ΔV = V − V 0 and V m = (V + V 0 )/2, where V and V 0 are the new and original areas (m 2 ), respectively, and Δt is the time step, s.
The numerical simulation for the CJSM first determines the influence of the parameters on the flow rate. The formulas where T i is the permeability coefficient; l s is the length of the seepage section, m; and l a is the block length, m. According to Darcy's law, Equation (4) can be transformed into a formula for permeability k j (m 2 /MPa s).
According to Equation (5), the joint stiffness and initial joint aperture of elastic coal samples and a persistent fracture coal sample ( Figure 2) can be calculated by combining the experimental curves of the stress permeability. It must be recognized that in many cases, reservoir rocks are anisotropic in permeability due to different pore or layering structures in different directions. 20,23,24 However, if the heterogeneity of reservoir rocks is considered, the general conditions of the axial and radial flow processes cannot be assumed to be the same. The same stress (axial and confining) in different directions has different effects on the contact hydraulic aperture. In general, using an isotropic model is intended to remove the influence of anisotropy in the study on effect of seepage direction on permeability stress tests. The parameters of the experimental coal samples (l = 0.05 m, l s = 0.05 m, l a = 0.002 m) and experimental conditions (ΔP = 0.2 MPa) are shown in Figures 3 and 4. Meanwhile, there is a residual joint hydraulic aperture a res to ensure that the joint aperture is always greater than zero. The selection of the stiffness and initial aperture for the persistent fracture is based on Equations (1) and (2), as well as the permeability stress experimental results of persistent fracture coal samples. Therefore, the results of the above method are a good match to laboratory test results under the condition of a relatively low effective stress (within 5-10 MPa). However, the stress sensitivity of the model permeability is zero at the high-stress stage because of the residual hydraulic aperture. This results in a certain difference between the stress sensitivities of the numerical simulation and experimental test in the high-stress state. Of course, the experimental results show that the stress sensitivity of the permeability is very small under the condition of high stress. Therefore, the CJSM is still useful in most cases, and the parameters used for a simulation do not need to change, which allows for a faster calculation. However, in the case of a coal seam with a large buried depth, the in situ stress is high. If the CJSM is adopted, the hydraulic aperture of the coal seam joint will remain almost at the residual hydraulic aperture. The main reason for this phenomenon is that the joint stiffness is a constant value in the simulation process. In other words, the compressibility of the joints is invariable. When the axial pressure or confining pressure exceeds a certain stress under the condition of deviatoric stress, the hydraulic aperture does not change and is equal to the residual hydraulic aperture, which directly affects the precision of the calculation, especially for fractures with angles, as shown in Figure 5. Under the condition of different axial and confining stresses, the joint inclination angle directly affects the stress on the joints.
The upper part of Figure 5 shows the aperture changes of cutters with the CJSM. As seen in the figure, a horizontal stress of 10 MPa will cause the joint to immediately enter the residual hydraulic aperture state. After that, the vertical stress has no influence on the joint hydraulic aperture. In other words, for the hydraulic aperture, a vertical stress of 5 MPa is equal to a vertical stress of 0 MPa under the condition of 10 MPa horizontal stresses. This is obviously different from the actual situation. The stiffness of the joint is further increased by an increase in stress under the varying joint stiffness condition. Thus, it is far from the residual hydraulic aperture in the case of a horizontal stress of 10 MPa. At the same time, an increase in the axial stress can further decrease the hydraulic aperture, as shown in Figure 2. Therefore, the results of a simulation with the VJSM are closer to those found in an experimental test, but in the calculation process, the joint stiffness should be updated according to the stress state, which decreases the simulation speed.

OF VARYING JOINT STIFFNESS METHOD
The experimental results show that the compressibility of a crack changes with a change in the effective stress. Better data matching can be obtained using the mean fracture compressibility c f , as suggested by Mckee et al 25 and Chen et al 26 : where c f0 is the initial state of fracture compressibility, and α f is the changing rate of fracture compressibility with the effective stress. The relationship between the compressibility coefficient and the joint stiffness is as follows: The initial state stress σ 10 is assumed to be zero, and Equation (7) can be simplified as follows: Equation (8) is brought into Equation (5) and Equation (1), and the formula of permeability can be further obtained as follows: where k j is the permeability of the joint, and k f is the permeability of the persistent fracture. In order to match the laboratory research, the parameters are l = 0.05 m, l s = 0.05 m, l a = 0.002 m, and ΔP = 0.2 MPa. The above parameters are introduced into Equations (9) and (10) to obtain a formula for the permeability (md), stress, and joint normal stiffness.
According to Equation (11), a 0 and k n can be obtained by fitting the measured data. The fitting results for elastic coal samples and persistent fracture coal samples are shown in Figures 6 and 7, and Table 1.
It can be seen from the fitting curves and correlation coefficient R 2 (more than 0.99) that the fitting effect with the VJSM is better than that with the CJSM. In addition, the fitting range for the stress is wider: the fitting range of the CJSM is generally less than 10 MPa, while the fitting range of the VJSM is much larger than 20 MPa. This indicates that the VJSM results are closer to the experimental results and have a higher reliability. The curves of the joint stiffness and hydraulic aperture with changes in the effective stress using the average value of the two kinds of coal samples are shown in Figure 8.
As can be seen from Figure 8, with an increase in the effective stress, the joint stiffness gradually increases, and the hydraulic aperture gradually decreases. The stress sensitivity of the joint stiffness slowly increases with the stress, while the stress sensitivity of the hydraulic aperture sharply decreases with an increase in the stress, but is almost unchanged under high-stress conditions. In fact, it is not necessary to set the residual hydraulic aperture according to Figure 8. In order to meet the needs of a numerical simulation, the hydraulic aperture of the 100 MPa effective stress is set as the residual hydraulic aperture, and the residual hydraulic apertures of the elastic and persistent fracture coal samples are 1.17 × 10 −6 m and 1.22 × 10 −5 m, respectively. As seen in Figure 8, the hydraulic aperture at 100 MPa is basically equal  to that at 20 MPa. In this study, a numerical simulation of the stress permeability was performed using the average fitting data and embedding Equation (8) in the UDEC with the FISH language. The simulation results are shown in Figure 9.
In order to reduce the effects of contingency, each condition (from model setup to model solution) was simulated at least 5 times for averaging in this paper. The laboratory measured results in Figure 9 show the average permeability of the coal samples. As can be seen from Figure 9, the numerical simulation results are in good agreement with the measured results. This indicates that the parameters selected by the VJSM can match the experimental results very well, and these parameters can be further used to study the gas seepage characteristic during pressure relief mining.
In the previously mentioned numerical simulation, the axial pressure was equal to the confining pressure. Thus, the normal stress acting on the joint was not related to the    Figure 11. The joint stiffness changes with the axial stress or confining stress. Because of the different inclination angles of the joints, the influences of the axial and confining stresses are different. When the confining stress increases, the stiffness of the joints in the vertical direction of the model rapidly increases, while the stiffness of the joints in the horizontal direction increases more gradually. The changes in the joint stiffness with an increase in the axial pressure are just the opposite. Thus, the hydraulic aperture at a certain angle does not reach the residual hydraulic aperture, increasing the reliability of the simulation.
To further verify the correctness of the simulation, numerical simulations and laboratory experiments were carried out under the conditions of confining stresses of 3 MPa, 5 MPa, 8 MPa, 10 MPa, and 15 MPa and axial stresses of 1-15 MPa. A comparison between the measured results and numerical results is shown in Figure 12. As can be seen from Figure 12, the numerical simulation results are basically consistent with the experimental results under deviatoric stress conditions, especially in the case of a high axial stress or high confining stress. Thus, discrete element numerical simulation using the VJSM can be used to simulate the relationship between the stress and permeability in a greater stress range and provides a research method for studying the influence of the fracture structure on the stress sensitivity of the permeability.

| DISCUSSIONS AND CONCLUSIONS
The fluid-solid coupling simulation of fractured coal and rock is helpful to study the influence of the fracture structure on the stress sensitivity of the permeability. Therefore, the key factor to improve the reliability of numerical simulation is selecting numerical simulation parameters that match the measured results. This paper firstly illustrated the use of an inversion method to select the fluid-solid coupling parameters by the CJSM. This method can be used to simulate the measured results in a lower stress range. However, because of the constant joint stiffness of the joint or fracture, it is necessary to set the residual hydraulic aperture to ensure that the hydraulic aperture is always positive. Permeability/md Axial stress/MPa This results in zero stress sensitivity for the permeability when the stress reaches the critical stress. Under the deviatoric stress condition, once the axial stress or confining stress is greater than the critical stress, the other stress has no influence on the permeability. In order to solve the above problems, a numerical simulation method for fluid-solid coupling by the VJSM was proposed based on the characteristic that the compressibility of a fracture decreases with increasing stress. In this method, the compressibility of the hydraulic aperture can be changed by increasing the joint stiffness during the loading process. Thus, the hydraulic apertures can change continuously, and there is no need to set the residual hydraulic apertures. It can not only greatly improve the simulation accuracy, but also make the application range wider than CJSM. At the same time, the changing formulas of joint and facture stiffness with stress are given according to the results of laboratory permeability tests (Equation11). It can be applied to numerical simulation of fluid-solid coupling in fracture development during hydraulic fracturing, pressure relief mining, and other engineering practices. The main process of fluid-solid coupling simulation with VJSM is shown in Figure 13.
Embedded the Equation (11) into the numerical simulation software with FISH language, fluid-solid coupling simulations under different stress states were performed using the VJSM in this study. The numerical simulation results were basically consistent with the experimental results, which indicated that the results of a fluid-solid coupling simulation using the VJSM could match the experimental results very well under the different stress states. In addition, the proposed method is also applicable to three-dimensional numerical simulation; as long as the hydraulic apertures of the lines are converted to the contact faces hydraulic apertures.