Abstract
Trying to reveal the mechanism of gas seepage in coal is of significance to both safe mining and methane exploitation. A series of FEM numerical models were built up and studied so as to explore the mesoscale mechanism of seepage in coal fractures. The proposed mesoscale FEM model is a cube with micron fractures along three orthogonal directions. The distribution of velocity and pressure under fluid-solid coupling was obtained, and furthermore, the seepage flow flux and an equivalent permeability of the whole model were calculated. The influences of fracture width, outlet velocity, and in situ stress level on seepage were investigated. The numerical results show that nonlinear Darcy seepage occurs during low velocity zone. The permeability is increased linearly with the increasing of facture width and outlet velocity. A certain change of lateral coefficient of in situ stress also affects seepage. The permeability is increased sharply once deviating the isotropic spherical stress state, but it is no longer changed obviously after the lateral coefficient has been increased or decreased more than 20%. The mesoscale seepage mechanism in coal fractures has been preliminarily revealed by considering fluid-solid coupling effect, and the key factors influencing fluid seepage in coal fractures were demonstrated. The proposed methods and results will be helpful to the further study of seepage behaviour in coal with more complex structures.
1. Introduction
In the process of coal mining, there is a kind of unconventional natural gas stored in coal seams, commonly known as coalbed methane (CBM). As an important energy in coal, environment pollution produced by using CBM is relatively much less than directly using coal resources, so its status is more and more important [1]. On the one hand, as an energy resource, CBM has a high utilization value; on the other hand, the outburst of gas during coal mining has become a serious geological disaster. Therefore, it is of great engineering significance to study the migration and seepage law of CBM in coal body and formulate reasonable extraction measures.
The migration and seepage laws of CBM are closely related to the deformation of coal matrix. Chen et al. [2] and Peng et al. [3] pointed out that stress and strain have a great influence on coal permeability. Geng et al. [4] and Cheng et al. [5] systematically studied the relationship between the stress sensitivity coefficient and permeability changes of briquette with different particle sizes. In coal reservoirs, the migration channel in the reservoir is regarded as a system of pores and fractures, and the permeability of this dual-pore and fracture structures has become one of the indicators of the success of CBM extraction engineering. The coal-bearing basin has undergone various periods and various degrees of compression and extension and other geological factors, which have reshaped the pore and fracture morphology of the coal, which have an important influence on the permeability of the coal seam. In addition, during the loading process of coal and rock, the microfracture structure inside coal and rock will change significantly [6–10]. From the perspective of acoustic emission, it shows that the cracks in the coal and rock masses gradually expand with the increase of stress [11, 12]. Meng et al. [13] conducted a three-dimensional mesoscopic simulation study on the crack propagation inside the soil-rock mixture. In the process of CBM extraction, the change of internal pressure will cause the stress of coal to change, which in turn affects the characteristics of coal pores and cracks. Due to changes in the structural characteristics of pores and fissures, the deformation of coal will have an impact on fluid seepage [14]. Therefore, in the process of CBM extraction, fluid-solid coupling will be quite prominent. This means that the influence of interaction between fluid and solid matrix must be added in the process of various simulation studies. After recognizing this problem, most scholars have fully considered the influence of fluid-solid coupling on the seepage of CBM and other fluids and have achieved certain research results [15–18]. Majewska et al. [19] found that coal experienced a contraction trend after the initial expansion of adsorption, which was attributed to the compression effect of injection pressure. Yin et al. [20] studied the influence of different adsorbent gases on coal deformation and permeability and pointed out that the stronger the adsorption under the same gas pressure, the greater the deformation of the coal and the lower the permeability. Based on the double pore structure characteristics of coal and considering the influence of moisture on the adsorption characteristics of coal and rock, the seepage model under the condition of solid-liquid-gas coexistence was established [21]. Mitra et al. and Gu et al. [22–24] studied the non-Darcy phenomenon of gas migration in coal seams based on simplified coal seam physical models. These models mainly consider the compressibility of cleats and the impact of coal matrix shrinkage on the cleat opening. Shi et al. [25] conducted a large number of experimental studies and theoretical model derivation on the direction of coal matrix desorption-induced deformation and permeability changes under effective stress and tested coal under different uniaxial strain, displacement control, and pore pressure. Based on the mechanical properties of the coal sample and the deformation characteristics of the porous medium, a mathematical model of CBM flow solidification under various conditions is obtained.
According to the research results of many scholars mentioned above, both the experimental results and the numerical simulation results show that the seepage law of CBM in coal is subject to the interaction of fluid and solid, which is due to fluid-solid coupling. The study of fluid-solid coupling mechanism between fracture and solid matrix is helpful to better understand the seepage law of CBM in coal. However, most of those researches are resulted from phenomenological coupling between the pore and the solid matrix at the macro scale. Therefore, some mesoscale FEM models will be built to represent the characteristics of fluid-solid coupling, and the seepage law of CBM along microfractures will be investigated by changing different geometry features and boundary conditions.
2. Theoretical Equations and Model Construction
2.1. Equation for Fluid-Solid Coupling Calculation
The fluid-solid coupling effect in coal and rock mass is a kind of interaction between seepage field of fluid and stress field of solid, which belongs the interdisciplinary of solid mechanics and fluid mechanics. Therefore, the influence of fluid and solid interaction can be derived by coupling flow law in fractures and deformation law of coal matrices.
2.1.1. Governing Equations for Fluid Calculation
Fluid flow should follow the basic conservation principles, including the law of mass conservation, momentum conservation, and energy conservation law. If the fluid includes other different components of the mixture, the system also follows the component conservation law. The general compressible Newtonian conservation law is described by the following governing equations: in which is time, is volume force vector, is fluid density, is fluid velocity vector, and is shear force tensor, which can be expressed by the following equation: in which represents pressure, represents dynamic viscosity, is unit tensor, and represents velocity stress tensor.
2.1.2. Solid Governing Equation
The equation of conservation of solids can be derived from Newton’s second law:
The deformation of the solid matrix is assumed to be elastic, following Hooke’s law:
In Equations (4) and (5), is the density of the solid; is the Cauchy stress tensor; is the local acceleration vector in solid domain; are the principal stresses; are the principal strains, respectively; and and are the Lame constants.
2.1.3. Fluid-Solid Coupling Equation
On the fluid-solid coupling interface, the fluid pressure and solid stress () and the displacements () must be equal or conserved; that is, the following equations should be satisfied: in which the parameters of the lower corner band represent the characteristics of solid parameters, while the parameters of the lower corner band represent the characteristics of fluid parameters.
2.2. Fluid-Solid Coupling Analysis Method
At present, solving control equations of fluid-solid coupling problems are basically divided into two categories: direct coupled solving and separated solving. Direct coupled solving means integrate fluid constitutive equations, and solid constitutive equations combine into one equation matrix for solution, so the fluid constitutive equation and the solid constitutive equation are solved together in a solver. Although direct coupled solving method is perfect and advanced in theory, its calculation process is time-consuming and needs too large memory. Therefore, the direct coupled solving method is only suitable for some very simple situation. Separated solving does not require a fluid-solid coupling control equation. It calculates fluid control equations and solid control equations in one solver alternately or two different solvers parallel; hence, the solution can be obtained after several iterations. The sequence solving the flow field and the solid field need to be specified in advance. Although the synchronous solution is difficult to converge because the energy on the fluid-solid coupling interface cannot be completely conserved due to time lag, the separated solving can fully utilize the existing procedures of current computational fluid dynamics and computational solid mechanics, so as to retain the module of available computing program and reduce the required memory greatly. In this study, a separated solving method was adopted to handle the fluid-solid coupling problem.
The pressure of fluid field in fractures will act on solid field and result matrix deformation; meanwhile, the deformed matrix will transform the geometry of fractures and cause a change of fluid flow which means the pressure and velocity of the fluid will be altered. Therefore, bidirectional fluid-solid coupling calculation needs to be considered. At the same time, it is necessary to impose large deformation effect of solid field.
2.3. Model Constructing and Parameter Setting
A cube with fractures along three orthogonal directions was constructed as shown in Figure 1. The side length of the model is tens of microns, and the width of the fracture is several microns. The fracture domain is designated fluid field indicated by the yellow part, and the residual matrix domain is designated solid field indicated by the blue part.
The flow field mesh needs to be divided more densely. In comparison, the solid field mesh can be appropriately relaxed. In addition, considering the calculation rate and convergence effect, dense tetrahedral meshes are used for the solid field near the fluid-solid coupling surface, and hexahedral meshes are used for the solid field away from the fluid-solid coupling surface. The elements and boundary conditions of one model are shown in Figure 2.
(a)
(b)
(c)
(d)
The fluid-solid interfaces are set as nonslip surfaces. The seepage channels are designed as three directions along the fractures, and hence, the inlets and outlets are assigned on different surfaces as shown in Figure 2(c). In order to simulate various seepage situations of these models, the kinds of pressure differences are assigned and studied. Three groups of inlet and outlet boundaries are defined as follows: (1)The seepage along -axis-negative direction is from top (named as inletz) to bottom (named as outletz)(2)The seepage along -axis-positive direction is from front (named as inlety) to back (named as outlety)(3)The seepage along -axis-positive direction is from left (named as inletx) to right (named as outletx)
The solid matrices are denoted as coal, in which elastic modulus and Poisson’s ratio are 1.4 GPa and 0.3, respectively. The in situ stress on each surface of a cube needs to be divided into matrix stress and fluid pressure, which can be calculated according to the area ratio of solid matrix and fluid field: where and are the in situ stress and the total area of a cube side surface, respectively; and are the effective stress and the area of surface solid matrix; and and are, respectively, the fluid pressure and the area of surface fracture fluid. Commonly, the in situ stress adopted in this simulation is selected as 5 MPa.
3. Results and Discussion
3.1. Influence of Geometric Width of Fractures
In order to study the relationship between fracture structure and seepage at mesoscopic scale, it is necessary to investigate various models with different fracture width. Five groups of models were created under the same boundary conditions. Each cube has a size of , but the fracture width in these five groups of models was 2 μm, 3 μm, 4 μm, 5 μm, and 6 μm, respectively.
The pressure on each inlet is set as 5 MPa. The inletz⟶outletz channel is set as the main flow channel. The velocity on outletz is set as 0.1 m/s, and the velocity on other two outlets is 0.9 . The differential pressure of each flow channel between the inlet and the outlet can be calculated by subtracting the inlet pressure from the outlet pressure obtained from the numerical simulation results. The differential pressures under various fracture width are shown in Figure 3.
The pressure difference gradually decreases with the increasing of fracture width. The variation of differential pressure with fracture width increasing is basically coincident along three different directions, i.e., inletx-ouletx, inlety-oulety, and inletz-outletz. An inflection point at the fracture width of about 4 μm can be noticed from Figure 3. When the fracture width is greater than 4 μm, the differential pressures along three directions are almost same.
According to the traditional Darcy’s linear seepage law, the pressure difference should be a certain value under the specified velocity. In order to explain the aforementioned variation of pressure difference, both the deformation of coal matrix and the flow of fracture fluid need to be learned.
Figure 4 shows the total displacements of coal matrix near fractures with different width. The displacements under large fracture width are smaller than those under small fracture width. It means the impact of matrix deformation becomes weaker with the increasing of fracture width.
(a)
(b)
The pressure on each outlet decreases evenly from one side to another side as shown in Figure 5(a). Therefore, the pressure on the fluid-solid interface (named as FSI) changes gradually along the diagonal direction as shown in Figure 5(b).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
It can be seen from Figure 5(a) that the two horizontal flow channels are similar except that their flow directions are different, so a section parallel to the YOZ plane located at μm can be selected to represent fluid flow, which is the middle plane of the inlety⟶outlety flow channel (named as MPF). The velocity components along the X, Y, and Z directions on MPF are shown in Figures 5(c)–5(e), respectively, and Figure 5(f) shows the contour of total velocity.
The velocity contours on MPF indicate that flow in fractures is irregular. The velocity mostly appears negative, so the fluid along Z direction is mainly from the top inlet to the bottom outlet. On the whole, the vertical flow along the Z direction is relatively slower near the left horizontal inlet than in the right part. It is noticed that the velocity of fluid near the horizontal inlet or the horizontal outlet is smaller owing to the restrictions of horizontal boundary conditions. The velocity mostly appears positive, and thus, the fluid along the X direction is mainly from the left inlet to the right outlet. On the whole, the horizontal flow along the X direction is relative slower near the top inlet than in the bottom part. The velocity mostly appears almost zero except in the middle intersection of two horizontal flow channels. It infers that some turbulent flow may occur in this middle zone as shown in Figures 5(e) and 5(f). The total maximum velocity appears in the corner between the inlet and the outlet, and the minimum velocity is in the corner between two inlets. The velocity decreases and increases gradually around the maximum and the minimum, respectively, and then an almost fixed value was retained which is distributed on the whole right-bottom part near the two outlets.
The distribution of fluid velocity on bottom outlet is shown in Figure 5(g). Since the outlet velocity has been specified as boundary conditions in advance, its distribution obeys laminar flow as expected. The velocity changes from the maximum in the middle to zero on both sides, and the average value is just the specified value.
Figure 5(h) shows the distribution of fluid velocity on top inlet. The changes along the X direction and Y direction are similar. The maximum velocity appears near the horizontal outlet, around which the velocity decreases continuously to zero. The velocity near horizontal inlets also retains zero owing that no pressure alteration occurs over there. The velocity distributions in fractures with different width are similar. However, the impacts of nonslip solid-fluid interfaces and the turbulent flow in middle zones become weaker with the increasing of fracture width, so the flow capacity is enhanced in large fractures. The actual velocity on each inlet or outlet can be derived from the numerical simulation results. Furthermore, the flow flux of inlet or outlet can be calculated by surface integral of these velocity values. Figure 6 shows the flow flux of each inlet and outlet under different fracture widths. The flow flux increases with the increasing of fracture width. The increasing multiple of flow flux is not equal to the increasing multiple of the inlet or outlet area. Consequently, the impact of fracture width on the fluid flow cannot be neglected.
According to the Darcy’s seepage law, an equivalent permeability of fracture seepage can be defined as in which is fluid viscosity, is flow flux, is model length from the inlet to the outlet, and represents differential pressure between the inlet and outlet.
The equivalent permeability under different fracture width can be calculated according to Equation (8), where the flow flux is from each outlet. The results are shown in Figure 7.
It can be seen that the equivalent permeability increases almost linearly with the increasing of fracture width, which indicates fluid seepage becomes easier for larger fractures. As mentioned previously, the horizontal flows along inlety⟶outlety and along inletx⟶outletx are similar, so their equivalent permeability is also the same. Because the velocity on outletz is specified as a larger value than that on outletx and outlety, the equivalent permeability of inletz⟶outletz flow channel is greater than those of the other two flow channels. Therefore, such equivalent permeability will be influenced by both geometric width of fracture and flow pressure and velocity. It also means that the traditional linear Darcy’s seepage law will become nonlinear due to the fact that the coefficient, i.e., permeability, is no longer a constant.
3.2. Seepage Law of the Model under Different Outlet Velocities
A series of models with a fracture width of 4 μm were taken to investigate the influences of outlet velocities. The inlet pressures of the X, Y, and Z directions are still set as 5 MPa. The inletz⟶outletz channel is set as the main flow channel, and the corresponding outlet velocity is 0.01 m/s and 0.025 m/s, 0.05 m/s, 0.075 m/s, 0.1 m/s, 0.5 m/s, 1.00 m/s, 5.00 m/s, or 10.00 m/s. The outlet velocity of the two horizontal flow channels inlety⟶outlety and inletx⟶outletx is 0.9 .
The distribution of pressure and velocity under different outlet velocities is roughly the same to those contours shown in Figure 8 only except that the values increase with the increasing of outlet velocity.
(a)
(b)
Since the pressure on inlet has been specified as 5 MPa and the pressure on outlet can be obtained from the numerical simulation results, the differential pressure along each flow channel can be calculated by subtraction. Figure 6 shows the differential pressure between each inlet and outlet under different outlet velocities.
It can be seen from Figure 8 that the differential pressure on each flow channel is basically proportional to the outlet velocity. Especially when the velocity is relatively large, the relationship between the increments of differential pressure and the increments of velocity is linear according a constant proportionality factor. It agrees with linear Darcy’s law during these ranges. However, the relationship curve deviates from the straight line gradually when the velocity is relatively small. It indicates a phenomenon of nonlinear Darcy flow during low velocity zone.
The equivalent permeability under different velocity can be calculated according to Equation (8). The results are shown in Figure 9. The equivalent permeability increases linearly with the increasing of velocity. It is also noticed that the equivalent permeability of inletz⟶outletz flow channel is greater than those of the other two flow channels because the velocity on outletz is specified as a larger value than that on outletx and outlety.
3.3. Seepage Law under Different in-Situ Stress
In situ stress is the natural stress that exists in the stratum that greatly influences the mechanical behaviours of coal rock. According to current researches, it is generally believed that the formation of in situ stress is caused by the compression of continental plate boundaries, thermal convection in the mantle, and gravity and often disturbed by kinds of engineering activity. A large amount of measured data show that the ratio of the maximum horizontal principal stress to the vertical principal stress is about 0.5~5.0. In this study, a lateral coefficient is defined as the ratio of the principal stress in the horizontal plane to the principal stress in the vertical plane. The vertical stress is specified as 5 MPa, and the horizontal stress is specified as 3 MPa, 4 MPa, 5 MPa, 6 MPa, or 7 MPa; i.e., the value of is 0.6, 0.8, 1.0, 1.2, and 1.4, respectively.
The previous models have fractures with equal width in which a certain difference pressure must be applied between the inlet and the outlet so that fluid can flow. In order to investigate the influences of in situ stresses, the pressure on each inlet needs to be equal to that on the opposite outlet so as to retain the specified in situ stress. Considering fluid in fractures with equal width cannot flow under constant pressure, the fractures with various width were introduced as shown in Figure 10. A series of models with fracture which width changes from 6.0 μm to 2.0 μm were constructed. Both the inlet pressure and the outlet pressure of each seepage channel are specified to equal to the in situ stress applied on the corresponding direction.
The distribution of velocity under different stress and pressure is roughly similar to that mentioned previously. Figure 11 shows the calculated velocity on each inlet and outlet under different lateral coefficient. Furthermore, the corresponding equivalent permeability can be calculated according Equation (8) and shown in Figure 12.
When the in situ horizontal stress and the in situ vertical stress are not equal, the flow velocity of fracture fluids increases and the equivalent permeability of model becomes higher. A jump is noticed near the lateral coefficient of 1.0, which means the seepage fluid under the isotropic spherical stress is the weakest; i.e., the deviator stress can enhance the seepage in fractures obviously. But once the lateral coefficient has been increased or decreased more than 20%, its influences become very slight so that little variations can be noticed for higher or lower lateral coefficients. The seepage at the lateral coefficient of 0.7 and 1.3 is calculated as a supplement. Numerical results show that the equivalent permeability of inlety⟶outlety is 0.01264 μm2 and 0.01231 μm2 and the equivalent permeability of inletz⟶outletz is 0.09911 μm2 and 0.09993 μm2, respectively. These results validate the jump occurs near the lateral coefficient of 1.0. It means the deviator stress will enhance the seepage once deviating the isotropic spherical stress state but such enhancement is limited and no longer takes effect even further increasing or decreasing lateral coefficients.
4. Conclusions
A series of FEM models representing seepage in coal fracture were built up and studied by means of ANSYS software. The proposed model is a cube with tens of microns side width containing three orthogonal fractures with several microns width. By introducing fluid-solid coupling effect, the seepage law of gas in three-dimensional fractures was simulated and analyzed. The following conclusions can be drawn: (1)Even under the same specified outlet velocity, the pressure difference gradually decreases with the increasing of fracture width. It means the equivalent permeability calculated according traditional Darcy’s seepage law is not a constant for different fracture widths. With the increasing of fracture width, the impact of matrix deformation becomes weaker, and the impacts of nonslip solid-fluid interfaces and the turbulent flow in middle zones also become weaker. Therefore, the flow flux through fractures increases with the increasing of fracture width, but the increasing multiple of flow flux is not equal to the increasing multiple of the inlet or outlet area. The equivalent permeability increases almost linearly with the increasing of fracture width, which indicates fluid seepage becomes easier for larger fractures(2)The linear Darcy’s law is agreed with only when the outlet velocity is relatively large, while the nonlinear Darcy seepage occurs during low velocity zone. Although the differential pressure on each flow channel is basically proportional to the outlet velocity, it must be noticed that their relationship curve deviates from the straight line gradually when the velocity is relatively small. The equivalent permeability increases linearly with the increasing of outlet velocity(3)The influences of the in situ stress under different lateral coefficients on seepage are obviously nonlinear. A jump of the equivalent permeability is noticed near the isotropic spherical stress state, while little variations of the equivalent permeability can be observed for higher or lower lateral coefficients once the lateral coefficient has been increased or decreased more than 20%. Hence, the deviator stress will enhance the seepage, but such enhancement is limited
Data Availability
The data is limited to single case reports and falls under the protected health information. For more information regarding the data availability, please reach out to the corresponding author.
Conflicts of Interest
The authors declare no conflicts of interest.
Acknowledgments
This work was supported by the Second Tibetan Plateau Scientific Expedition and Research Program (STEP) (Grant No. 2019QZKK0904), the National Natural Science Foundation of China (Grant No. 51727807), and the Fundamental Research Funds for the Central Universities (Grant No. 2020YJSMT06).