A fluid‐solid coupling method for the simulation of gas transport in porous coal and rock media

Understanding coal and rock permeability, and the corresponding influence on stress, is important in the field of energy development. In applied engineering, there is a tendency to employ three‐dimensional methods that are simpler, less time‐ and cost‐expensive, less computationally expensive, and larger scale. Thus, a numerical simulation fluid‐solid coupling method is proposed in this paper. The proposed numerical simulation method utilizes the stress‐permeability test results for elastic and plastic coal samples. The permeability models of the elastic and plastic coal samples under loading and unloading were obtained by fitting the experimental results and embedded them into FLAC3D software by using FISH language. The results of the uniaxial and triaxial flow simulations are consistent with the experimental results, thereby confirming the accuracy and feasibility of the numerical method proposed in this paper. This allows the permeability of the numerical reservoir model to be continuously updated according to the current stress level in the production process.


| INTRODUCTION
Understanding coal and rock permeability, and the corresponding influence on stress, is important for processes such as oil and gas resource development, coal mine safety production, and nuclear waste storage. [1][2][3] The main research methods used to evaluate the permeability of overlying strata during coal-seam group mining are field measurement, laboratory simulation, and numerical simulation. Field measurement is generally regarded as the most direct and reliable way to describe the fracture evolution laws and permeability of overlying strata, and has been widely employed. The slug method is the most common method for testing the permeability and fracture development characteristics of overlying strata, especially for saturated coal seams 4,5 . In addition, pumping and flowmeter tests are used in field permeability testing. 6 Field measurements have many advantages and are the best way to obtain accurate information about an actual field. However, they are not widely performed because they are difficult, relatively expensive, and have a poor safety record and low feasibility. With the rapid increase of computing capacity, numerical simulation continues to garner more attention in the field of engineering application. Its cost is relatively low, it is flexible and direct, and it can quantitatively analyze the relative influence of each factor; furthermore, it yields results that are highly consistent with the corresponding measured results. Understanding the respective influences of stresses on the overall permeability and flow pathways is critical to the development of a numerical model of the coupled hydromechanical process of fractured rocks in various fields of rock engineering. [7][8][9][10] Considering this, the aforementioned methods can be broadly categorized as either a continuum model-based large-scale field engineering simulation, or a discrete element model-based laboratory-scale inversion simulation. The discrete element model explicitly simulates the deformation and flow in rock masses, and can easily depict flow anisotropy and heterogeneity owing to the random distribution of fractures and their hydraulic apertures. 11,12 This type of model has been well developed in recent years. The basic assumptions are that the rock matrix is impermeable and linearly elastic, and that the fluid only flows through fractures. 13,14 At present, DEM-based stress-flow simulations can generally be performed in two steps. The first step entails the use of a random fracture network that is generated as based on the statistical data from site investigations to obtain the representative elementary volume. In the second step, cubic laws are applied to use the obtained model to study the effects of stress on flow and transport. 15,16 Regarding application of the finite element method to practical engineering problems, the respective relationships between permeability and bending-, fractured-, and caved-zone stresses have been summarized in previous studies. [17][18][19][20] FISH language was scripted in FLAC and subsequently used to simulate the permeability variability of surrounding rocks during mining of a longwall face and to analyze the effects of gas extraction. Si et al. 21 used FLAC3D to simulate stress distribution in the extremely thick coal seams resulting from fully mechanized coal mining; they then used MATLAB to introduce the stress results to ECLIPSE300 to study gas migration. Although the DEM-based fractured-rock coupled hydromechanical simulation can yield results that well match the laboratory test results, it is rarely used in large engineering applications. 22 For an engineering problem, a more realistic three-dimensional numerical simulation is clearly better than a plane simulation; however, at present, most discrete element fluidsolid coupling simulations are two-dimensional. The finite element method applied to practical engineering problems has yet to be validated on a laboratory scale 24,25 ; this directly affects the reliability of the numerical simulation. Thus, in this study, a simple fluid-solid coupling simulation technique was developed as based on the results of laboratory tests. The coupling model was then verified by using the simulated uniaxial and triaxial compression results.

NUMERICAL SIMULATION
The steady-state permeability of coal is known to be highly stress-dependent. Laboratory experiments and field measurements performed at various burial depths have shown that, in general, there is an exponential decrease with increasing net confining/overburden stress. Depending on whether the model equations take into account changes in effective reservoir stress or coal-bed cleat porosity as the coal-bed reservoir is being depleted, the permeability models may be broadly classified as either stress or strain based. 26,27 Models, such as the Palmer and Mansoori (P and M) model, Shi and Durucan (S and D) model, Cui and Bustin (C and B) model, and improved P and M model, have been used to predict permeability change during the process of coal-seam gas extraction. [28][29][30] All models have been developed to describe permeability change under the uniaxial strain conditions that are commonly observed in coal-bed reservoirs. It was also assumed that the total vertical stress remains unchanged during coal-bed methane (CBM) recovery. However, how the permeability varies throughout the coal mining process significantly differs from how it varies during CBM extraction. The applied stress continuously changes throughout the process of coal-seam mining; this variability occurs in coincidence with the plastic yielding of coal. [31][32][33] The new permeability model for coal under elastic, plastic, and post-failure state is given by Chen et al. 34,35 . In this model, the authors stated When the coal turns into plastic state, the modified logistic function needs to be incorporated to account for the contribution of fracture growth and generation to the permeability surge. Once the deviatoric effective stress surpasses a certain level and then reaches a plateau after failure. This equation was derived from matchstick model ( Figure 1A) proposed by Seidle et al. 36 : where k f is the fracture permeability of the coal, md; σ 1 is the effective stress, MPa; c f is the compressibility of the fracture, and k f0 and σ 10 are the initial states of the permeability (md) and effective stress (MPa), respectively. However, the model proposed by Chen et al. 34,35 ignores a case in which the stress will decrease dramatically when plastic failure occurs. It is no significant decrease was observed of the deviatoric effective stress and mean effective stress in stress state. The method to model the permeability change from elastic state to plastic state by Si et al. 21 directly increases the permeability value to a specific level. The permeability models of the elastic and plastic state were the same as Equation .
For the plastic state, the initial permeability k f0 was enlarged several times (5, 10 or 20) to the elastic state. Besides, the permeability compressibility factor, c f , was decreased from 0.7 to 0.25, 0.45, and 0.65. However, it requires constant attempts to match field measured data to determine final k f0 and c f . In this paper, in order to determine the above parameters more accurately, laboratory experiments in two states have been carried out. Note that c f is not a constant value, but varies according to the effective stress. More accurate simulated results can be obtained by using the mean fracture compressibility, c f , as suggested by Mckee et al. 37 and Chen et al. 38 : where c f0 is the initial state of fracture compressibility and α f is the rate of fracture compressibility change as a function of effective stress. Recently, Chen et al. (2015) and Chen et al. 34,35 verified that the permeability model, as expressed in Equations and , is still applicable to rocks with irregular fracture network ( Figure 1B) and porous rock ( Figure 1C). Therefore, the unified exponential form model can be applied to any type of fractured or porous rock, and is not limited to idealized rocks. Thus, the above-described model can be used to fit the experimental data of the plastic coal sample in this study. Assume that the initial effective stress state for the permeability test is zero and substitute Equation into Equation ; the simplified formula is as follows: As can be ascertained from Equation , k f0 , c f0 , and α f can be obtained via laboratory tests. Therefore, by incorporating Equation into the numerical simulation, the permeability can be updated in real time according to the calculated stress.

| Coal sample selection and preparation
The coal samples were extracted from the Huainan coalfield in China. Raw coal samples with dimensions of ϕ50 mm × 50 mm were used to study the stress-permeability relationship. During coal-seam mining or triaxial compression testing, the coal will yield to the high stress. When the coal yields, a large number of fractures are generated, and permeability sharply increases.
(1)  Therefore, it is necessary to carry out stress-seepage testing on the coal before and after yielding. In this study, shear tests were performed on standard coal samples; the coal samples with a shear-type failure, specifically, a failure with a shear fracture plane that lies perpendicular to the coal sample edges, were selected for further experimental study of the coal sample in plastic state, as shown in Figure 2. The bias (related to accuracy) of the permeability test equipment is less than 5% within the measurement range of 0.001 mD-1000 mD.

| Loading and unloading stress paths
During uniaxial or triaxial compression, the stress of coal initially increases until the ultimate strength is reached; it then proceeds to gradually decrease. The stress-permeability relationship observed during loading differs from that observed during unloading. 39 Thus, it is necessary to obtain the permeability model of coal under loading and unloading conditions. In this study, the experimental stress paths of the loading-unloading cycles were constructed under the initial stress conditions (Figure 3). The maximal axial and confining stresses were increased to 20 MPa at intervals of 2 MPa, and the upstream gas pressure (ie the absolute pressure) was maintained at 0.2 MPa. The main boundary conditions are depicted in Figure 2. The in situ stress measurement results show that the horizontal stress was equal to the vertical stress. Thus, the experimental axial stress was equal to the confining stress.

| EXPERIMENTAL PROCEDURES
The coal sample was covered with a rubber sleeve and two O-rings and horse-shoe clamps on each top cap and base pedestals were placed over the rubber sleeve to avoid leakage. Before the permeability test, the gas (CH 4 ) pressure was applied to the predetermined value, and the sample was adsorbed to the adsorption equilibrium under the initial stress condition. Then, permeability tests were performed on coal sample at different stress path as presented in Figure 3. Once the sample was exposed to the target stress, the gas was injected from the bottom of the sample and collected from the top of the sample. The test was allowed to continue until the flow became stabilized. After flowing through the sample, the gas flow rates were recorded at every 5 s interval by the flow rate transducer. The test time of the each stress state was approximately about 30 min.

| Results analysis
The experimental results corresponding to the stress paths shown in Figure 3 are shown in Figure 4. Large permeability loss was observed in the elastic and plastic coal samples during the loading and unloading processes. When the stress was unloaded to the initial stress level, the permeability of the coal sample significantly decreased, especially for plastic coal samples. The permeability of the elastic coal sample during loading and unloading was 59.65% of the initial permeability. The loss of the plastic coal samples reached 84.33%. The main reasons were that mutual embedding of fracture surfaces and pulverized coal increased in the single fracture during the loading, providing a significant reduction of the permeability in the loading-unloading cycle, as shown  Figure 5. 39 At the beginning of loading, the two fracture surfaces are linearly compacted with elastic deformations ( Figure 5A). When the stress exceeds the point contact strength of the fracture, the fracture contacts were crushed and reconsolidated, the pulverized coal begins to fill the fracture space ( Figure 5B). Meanwhile, the self-adjustment by fracture face contact crush and friction slipping makes the mutual embedding of fracture surfaces and stabilized in a piled position ( Figure 5B). As more fracture contact is broken, the pulverized coal gradually fills full the fracture; the fracture surface is gradually stabilized and no longer sliding ( Figure 5C). In this stage, the void between two fracture surfaces substantially reduces which lead to a drastic drop in the permeability at the beginning of the loading. The experimental results were fitted by using Equation ; the fitted results are shown in Figure 4 and Table 1. Figure 4 and Table 1 show that Equation resulted in very good fitting, with the coefficient of determination R 2 exceeding 0.99. This confirms that the fitting formulas applied to the coal samples can be used for further calculations. Using the average values presented in Table 1, the permeability equations for the elastic and plastic coal samples can be expressed as follows: where subscripts E and S represent the elastic and plastic coal samples, respectively, and subscripts L and U denote loading and unloading, respectively.

CALCULATION METHOD
Based on the permeability model, given as Equations and , FISH language was used to write a simple algorithm to incorporate into the flow model. FISH language is a built-in FLAC3D programming language. The fluid-solid coupling mining simulation was carried out by using this method.

| Model establishment and parameter selection
The coal of the numerical model built by using FLAC3D had the same shape and size as the laboratory-tested cylindrical coal samples. The model had a height of 100 mm and diameter of 50 mm. The mesh contained 138 240 elements and 140 605 grid points, as shown in Figure 6. A constant velocity (or loading rate) of 2.5 × 10 −7 m/step was applied to the top surfaces in the axial direction. In the uniaxial compression test, the side surface of the numerical model was free; conversely, in the triaxial compression tests, uniformly distributed confining pressure was applied normally to the side surface. In addition to the boundary conditions introduced to the mechanical parameter calculations, boundary conditions were applied for the seepage calculations. A constant gas pressure difference of 0.5 MPa was applied to the top and bottom end faces in the axial direction. The top and bottom end faces were permeable, and the side surface of the coal sample was impermeable. In the model, the strain-softening failure criterion was used to evaluate the coal sample performance. The mechanical properties of the elastic coal sample are listed in Table 2. The cohesion and friction angles of the coal sample degraded as the plastic shear strain increased; these factors were assigned residual values (shown in Table 2) when the plastic shear strain reached 0.01. Using the mechanical properties in Table 2, the uniaxial compression test stress-strain curves were derived as shown in Figure 7. It can be seen that the curves derived as based on the simulated results well match those derived as based on the laboratory test results, and that the failure mode is consistent. This confirms that the mechanical properties given in Table 2 can be used to design large-scale simulations and triaxial compression and flow tests for engineering applications. It should be noted that the unit of permeability in the numerical simulation software FLAC3D is m 2 /(Pa·s), and 1 md equals 9.036 × 10 −11 m 2 /(Pa·s).

| Employing the fluid-solid coupling method by using FLAC3D
The permeability-stress models of the elastic and plastic coal samples were obtained by performing laboratory tests. Thus, in the numerical simulation, the permeability can be updated according to the current stress. Yes the block after identifying the mechanical characteristics. In addition, the permeability-stress relationship was also found to be different. Thus, the stress state of the model zone should be determined after the plastic state has been observed. The permeability was then updated according to Equations and . In this study, the flow calculations were carried out after the mechanical parameters were calculated, and the final permeability of the entire model was calculated by using Darcy's law. The specific program that was scripted in FISH language and implemented in the FLAC3D software is shown in Figure 8.

| RESULTS AND DISCUSSION
Uniaxial and triaxial compression flow tests with a confining stress of 2 MPa were carried out in the numerical simulation; the simulated results are shown in Figure 9.
As can be seen in Figure 9, the uniaxial and triaxial compression flow simulation yielded stress-strain curves and permeability-strain curves that are consistent with those obtained via laboratory testing. This confirms that the proposed model can well simulate uniaxial and triaxial compression flow experiments.
In the uniaxial compression elastic phase (point a to point c) in Figure 9, the permeability of the model decreased slowly from 13.65 md to 0.13 md with axial stress increasing from 0 MPa to 12.01 MPa. When the model entered the yield phase (point c to point d), the permeability began to increase from 0.13 md to 5.5 md. When the axial stress exceeded the peak value (15.79 MPa), the permeability increased significantly. When the stress decreases to 1.468 MPa, the permeability increases to 88.91 md, which is 6.5 times of the original permeability. In the triaxial compression test with a 2-MPa confining stress, the compressive strength and residual strength of the model increased obviously. The initial permeability and residual stage permeability are also significantly reduced. This is mainly due to the increase of initial confining stress. With the increase of axial strain, the change of permeability curve presented an "S" type change. This is basically consistent with the experimental results of Chen et al. 34,35 and the experimental results of the laboratory, indicating the accuracy and feasibility of the numerical simulation in this paper. To study the evolution law as it pertains to the relationship between stress-induced damage and permeability during the compression seepage process, the permeability, stress, and plastic state at five different points throughout the process were selected for further study ( Figure 9A); the results are shown in Figure 10.
During loading, the yield block (new fracture) of the coal sample gradually developed as the stress increased. In the early stage of uniaxial compression, the stress at the midpoint between the two coal sample ends significantly increased, thereby causing the plastic phase to dominate in this area. With the development of a plastic zone at the left and right sides, the stress gradually shifted to the middle of the coal sample, and the stress in the plastic zone began to unload. The permeability of the coal sample was found to be highly correlated with the plastic zone (fracture development) and stress, and the permeability of the elastic block gradually decreased as the stress increased. In addition, the permeability of the block increased as the block yield increased, but the increase in permeability was not significant because the stress within the plastic block was still relatively high, as shown in Figure 10B. Furthermore, because the plastic zone in Figure 10B was small and does not go through the model, the permeability of the entire model gradually changed. With the gradual expansion of the yield zone toward the two sides, the stress significantly decreased and shifted to the center of the coal sample. This led to a sharp increase in permeability at both sides of the coal sample. However, the permeability of the center area continued to decrease, and the plastic zone was not thoroughly connected each other. This caused the permeability of the coal sample to begin to gradually increase, but maintain a small value. When the stress on the coal sample reached its peak value, the plastic zone connected each other, and a tensile fracture occurred at the junction of the elastic and plastic zones (similar to the X-type region). It shows obvious brittleness characteristics and is basically the same as the laboratory test results of Li et al. 40,41 The sharp decreases in stress and more block yield (new fracture generation) caused the overall permeability of the coal samples to sharply increase. In the residual stress phase, the continuous reduction of stress and development of new fractures caused the permeability of the coal sample to continuously increase. 22,42 The permeability of the full model was obtained by performing flow calculations. The gas pressure at the top and bottom of the coal sample was 0.5 MPa and 1 MPa, respectively. For the flow calculations, the pore pressure at every point in the model was monitored. The flow calculations were considered to be balanced when the pore pressure was constant. The results of flow calculations for points a and d of the uniaxial compression curve are shown in Figure 11. The fullmodel permeability at points a and d is, respectively, shown in Figure 9. The permeability at point d is noticeably higher than that at point a. Furthermore, the pore pressure distribution at these two points is also different. The pore pressure at point a was observed to be parallel each other, except at the upper and lower ends. In contrast, there was very irregular pore pressure distribution at point d. By comparing the permeability distribution at each point, the pore pressure distribution was found to be relatively dense at low permeability. This means that higher permeability corresponded to thinner pore pressure distribution. This is because lower permeability requires a larger pore pressure difference to achieve the same flow distance.
Thus, the fluid-solid coupling method proposed in this paper can accurately represent the laboratory uniaxial and triaxial compression flow experiments to quantitatively analyze the relationship between the stress, plastic state, and permeability of coal samples. The proposed method can provide a basis for developing fluid-solid coupling simulations for coal mining and gas drainage.

| CONCLUSIONS
Regarding engineering applications, there is a tendency to use three-dimensional methods that are simpler, less cost-and time-expensive, larger scale, and require less computation as compared to DEM. However, the accuracy of these numerical methods may be lower than that of DEM. Therefore, a simple fluid-solid coupling model was developed to simulate uniaxial and triaxial compression. As compared to existing numerical methods, the proposed method is faster and requires fewer parameters. Moreover, by using the proposed method, the permeability can be evaluated during coal-seam mining, and fluid-solid coupling engineering problems can be calculated. Figure 12 shows the establishment and application of fluid-solid coupling simulation method. The main conclusions of this study are as follows: (1) Stress-permeability tests were carried out on elastic and plastic coal samples. The permeability models of elastic and plastic coal samples under loading and unloading were obtained by fitting the experimental results.
(2) FLAC3D was used to develop a simple fluid-solid coupling method. The specific program was scripted by using FISH language and implemented in FLAC3D. With the proposed method, the permeability of the elastic model, and that of the plastic model, can be updated according to the stress as the mechanical parameters are calculated.
(3) Uniaxial and triaxial compression flow simulations were carried out by using the proposed fluid-solid coupling method. The simulated results were consistent with the experimental results, demonstrating the accuracy and feasibility of the numerical models proposed in this paper.
(4) The evolution law as it pertains to the relationship between stress-induced damage and permeability during the compression flow process was illustrated via the results of numerical simulation. The permeability of the coal sample was highly correlated with the plastic zone and stress. The stress within the plastic zone was transferred to the elastic zone after the block yielded; subsequently, the permeability of the plastic zone significantly improved. When the plastic F I G U R E 1 1 Gas pressure results for the uniaxial compression