Abstract

In order to investigate the influence of pore pressure on hydraulic fracturing behavior in the local and whole model, the coupled flow-stress-damage (FSD) analysis system RFPA-Flow was used to study the influence of rock heterogeneity, natural stress ratio, double-hole spacing, and water pressure gradient on the stress shadow effect. The numerical results show that the tensile crack induced by pore water pressure is significantly affected by the pore water pressure and water pressure gradient. The larger the pore pressure gradient is, the more asymmetrical the crack development pattern and the smaller the instability pressure of the model. In addition, the shape of hydraulic fracture becomes much more irregular with the increase in rock heterogeneity. The number and shape of tip microcracks under the influence of local water pressure are closely related to the homogeneity of rock. Moreover, when the natural stress difference is large, the hydraulic fracture propagates parallel to the maximum principal stress; when the stress field is close and the spacing between two holes is less than 5 times the diameter, the propagation direction of hydraulic fractures between holes is perpendicular to the maximum principal stress. It is found that no hydraulic fractures occur between the two holes when the distance between holes is greater than 5 times the diameter.

1. Introduction

High-pressure fluid fracturing and mechanical fracturing are the two main approaches to rock cracking, and there are many research studies on the mechanical cracking mechanism of rock. Tests in this laboratory mainly include the uniaxial compression test, triaxial compression test, and tensile test, and the numerical simulation research study is also relatively mature [14]. However, the mechanism of high-pressure fluid fracturing and how the pore water pressure affects the tensile failure mechanism and crack morphology of rock are not well understood. Some field evidence has shown that the influence of pore water pressure on the mechanism of crack initiation and propagation is sometimes unclear. The principle of effective stress is extended to the problem of reservoir fracture, suggesting that the pressure of hydraulic fracture initiation and propagation may decrease with the increase in reservoir pore water pressure. The reopening pressure and propagation pressure of cracks show a decreasing trend when low-rate fluid is injected, which may be because a large amount of fluid has been lost to the reservoir [5]. This conclusion is consistent with the results of laboratory tests. When a large amount of fluid is lost to the rock, the formation fracture pressure will be significantly reduced [6, 7]. Some field tests show that the hydraulic fracture growth pressure increases at high reservoir stress and decreases at low reservoir stress, but this is inconsistent with the principle of effective stress [8, 9]. The reason leading to this phenomenon is that the influence of pore pressure on the crack evolution should be analyzed from both the local and global model. When the local pore water pressure at the crack tip increases, the crack propagation will be enhanced. The overall increase in pore water pressure due to the formation of compressive stress will resist the development of cracks. At the same time, the effect of fluid flow due to pore pressure difference must be taken into account.

Some scholars have discussed the influence of symmetrical or asymmetrical pore pressure distribution on flat crack propagation [1012]. A more complex problem in the process of hydraulic fracturing involves the propagation and evolution of cracks in the field of nonuniform pore water pressure. The propagation of hydraulic fractures is most obviously affected by the internal pore water pressure and stress ratio. Hydraulic fractures in reservoirs generally propagate in the direction parallel to the maximum far-field stress, and the energy driving crack propagation decreases continuously during the propagation process [1, 6]. However, under a relatively uniform natural ground stress field (stress ratio close to 1), the stress field cannot play a dominant role on the propagation of the crack. In this case, the pore pressure difference of the crack tip in the local scope may affect the crack propagation path. The pressure gradient of the stress field in the global scope can influence the crack propagation direction [10]. According to the studies of Geertsma [10] and Cleary [13], the influence of pore water pressure on crack evolution must be studied at local and global scales for hydraulic fracturing nonlinear behavior. The increase in local pore water pressure at the crack tip will drive the continuous development of the crack. However, the overall increase in pore water pressure caused by the increase in in situ stress field may resist the propagation of hydraulic fractures, especially tensile cracks. In addition, Detournay et al. [14] studied the initiation and propagation processes of the tensile crack in double-hole and three-hole SLATE samples under asymmetric pore pressure through experiments, and the results showed that the local pore water pressure played a direct role in tensile crack propagation. Bruno and Nakagawa [15] conducted a numerical simulation of shale fracturing test of Bruno and Nakagawa by means of numerical methods and analyzed the progressive crack propagation-based system under asymmetric pore pressure distribution. According to Wang et al. [16], crack migration is mainly controlled by two parameters, one is the difference between far-site stress deviation and characteristic porous elastic stress, and the other is the stress disturbance caused by stress field at the crack tip and characteristic elastic stress. Different from the stress disturbance in conventional rock mass engineering [1719], the stress disturbance caused by stress ratio is critical to hydraulic fracturing propagation. When the second parameter is small enough, the propagation path of the crack can be completely predicted by the ground stress field. In addition, tensile failure of rock caused by internal water pressure is dominant, and shear failure will be restrained. When rock is subjected to internal water pressure and external load at the same time, it will show different failure modes. Therefore, the internal water pressure and in situ stress of rock are two important test parameters, which determine the initiation and propagation of hydraulic fractures. Warpinski and Teufel [20] found that multiple microcracks or branching cracks are likely to occur near fractured horizontal wells, and the generation and propagation of cracks in the adjacent wellbore will be influenced by each other, affecting the stability of the wellbore wall and the propagation path of cracks. At the same time, it is considered that the distance between adjacent boreholes has an obvious influence on the interaction between cracks. Therefore, the distance between the two phases and the wellbore is another important test parameter. Due to the inherent heterogeneity of rocks, pores, weak surfaces, cinders, and diagenetic minerals lead rocks to be absolutely heterogeneous materials. In the test and numerical calculation of hydraulic fractures, heterogeneity is another important parameter to be considered.

Bruno and Nakagawa [14] demonstrated the effect of pore water pressure on the initiation and propagation of tensile cracks through linear elastic fracture mechanics and laboratory tests. The crack evolution process is also affected by the local pore water pressure at the tip, the distribution of pore water pressure, and the pressure gradient in the whole range. Based on Griffith theory and the maximum strain energy density method, linear elastic fracture mechanics takes into account the influence of fluid pressure at the crack tip, stress intensity factor, and pore water pressure on the strain energy driving crack propagation [20]. In addition, the strength theory can also analyze the crack development to a certain extent. The stress state in rock is determined by the basic equations of elastoplastic mechanics. When the effective stress exceeds the tensile strength at a certain point in the material, the crack will initiate and expand, and the failure point occurs at the maximum stress point [21, 22]. Fracture mechanics and strength theory can help us to know and understand the principle of hydraulic fracturing, but the analytical solution derived from the theory is only applicable to some simple cases. In addition, laboratory test methods bring great inconvenience to control various test parameters affecting the behavior of hydraulic fracture propagation, and it is also very difficult to measure the pore water pressure caused by crack evolution in the process of heterogeneous rock fracturing [14].

For some complex problems, the numerical calculation method can be used to better study the complex mechanical behavior of hydraulic fracturing. In this study, the two-dimensional rock instability program of RFPA-flow, considering the coupled flow-stress-damage was used. Five numerical models are designed to explore the local and whole pore pressure evolution considering the heterogeneity of the rock, in situ stress ratio, double-hole spacing, and pore water pressure gradient. To a certain extent, the numerical results fill up the insufficiency of this study on hydraulic fracturing with a double-hole model and provide a basis for laboratory test design and field gas shale development.

2. Numerical Model and Procedures

The RFPA (realistic fracture process analysis) -Flow numerical code in many literature studies has carried on the detailed introduction [23, 24], the program considering the seepage and stress in the process of rock damage to failure features, and it is based on the consolidation theory of Biot and modified effective stress principle of Terzaghi [25]. In addition, the RFPA method considers the effect of rock damage on permeability changes before and after rock failure, and it can describe the nonuniform permeability in the process of hydraulic fracturing to simulate the fracturing evolution behaviors. The reliability of this numerical model has been proved by many scholars [2325], and its detailed theory is not introduced in this study.

The calculation model of two-dimensional double-hole plane strain is shown in Figure 1. The geometric size of the model is 400 mm × 400 mm, and the whole model is divided into 220 × 220 = 48,400 units. Two small holes with a diameter of 30 mm were excavated on the horizontal line through the center of the model. The horizontal stress σH and vertical stress σV are applied on the boundary of the model, and the stress ratio is expressed as K = σH/σV. Hydraulic pressure force was applied inside the edge of the two holes, and fluid pressure increment Δp is 0.5 MPa, until the crack propagation is caused by sample macro damage. In this model, the influences of rock heterogeneity, hole spacing, injection rate, and other factors on the evolution characteristics of hydraulic fractures in the local and global models were comprehensively considered. The mesoscopic input parameters of the model are listed in Table 1.

Using the model established in Figure 1, five numerical simulators are considered in this study, as below:(1)Case 1: Consider the case where the double holes in the model are equal initial water pressures, the increment of water pressure is the same, and the external boundary is a uniform stress field (K = 1). In order to better understand the phenomenon of crack initiation, propagation, and coalescence caused by the distribution of pore water pressure, the model of pore size expansion and deformation is calculated.(2)Case 2: Model stress distribution, pore water pressure distribution, hydraulic fracture morphology, crack propagation path, and fracture pressure were discussed under different rock homogeneity (m = 1.5, 3, 5, 1000).(3)Case 3: Study the evolution formation of hydraulic fracture in the double-hole model under different stress ratios (K = 0.5, 1, 0.8, 1.5), discuss the influence of compressive stress field on pore water pressure in the overall range, and analyze the evolution law of crack.(4)Case 4: The response of hole spacing to hydraulic fracturing is studied. The value of double-hole spacing L is 60 mm, 100 mm, 130 mm, and 170 mm, respectively. The purpose is to explore the influence of double-hole spacing on local pore pressure distribution and the influencing factors driving the development pattern of hydraulic fracture.(5)Case 5: Discuss the propagation law of hydraulic fracture and its influence on the model rupture and instability under the condition of asymmetric water pressure injection in two holes.

The first, second, and fourth models are used to analyze the crack evolution law driven by the local pore water pressure at the tip of hydraulic fracture. The fourth and fifth models focus on the influence of hydraulic pressure gradient on the initiation, propagation, and coalescence of hydraulic fractures within the whole testing model. The third and fifth numerical models focus on the evolution of hydraulic fractures between holes, analyze the interaction mechanism between two holes, and discuss the mutual attraction effect of hydraulic fractures.

3. Numerical Results and Analysis

3.1. Hydraulic Fracture Propagation in a Two-Hole Model

Figure 2 shows the progressive failure process of the two-hole model under the condition of constant water pressure increment, stress σh = σv = 10 MPa, homogeneity m = 2, and pore pressure increment of 0.5 MPa. When the calculation reaches 30 hour step, almost no elements around the two holes are destroyed and there is no sign of crack initiation, which is the stress accumulation stage. When the calculation time step reaches 33, the minimum principal stress accumulates to the tensile strength of the material. Sporadic and scattered microcracks will appear around the two holes. The number of microcracks is increasing, the microcracks are distributed in an umbrella shape, and the microcracks accumulate in the middle area of the two holes. When the calculation time step reaches 35, the model starts to be unstable. When the rupture stress is 27.5 MPa and the water pressure reaches the rupture pressure, the crack continues to expand until the model is destroyed under the condition that the water pressure of the two holes does not increase.

The propagation direction of the hydraulic fracture in the left hole (hole 1) is almost horizontal, and the angle between the propagation direction and the horizontal direction of the hydraulic fracture in the right hole (hole 2) is about 45°. The hydraulic fractures on both sides of hole 1 and hole 2 are, respectively, symmetrical to the pressurized hole. However, due to the interaction between the stress field and the two increasing holes, the propagation length of the hydraulic fractures in the middle of the two holes is larger than that in the outside of the two holes. When it reaches the 35-5 time step, the hydraulic fractures near the left and right pressure holes begin to turn and have a tendency to connect. When it reaches 35 to 12 hours, the hydraulic fractures between the two holes are completely connected in the horizontal direction, which is parallel to the direction of the maximum horizontal principal stress. The calculation model shows that the distribution of pore water pressure has a great influence on the initiation, propagation, and coalescence of cracks, and the coalescence of cracks between double holes is attributed to the distribution of pore water pressure between pores.

For the convenience of comparison and illustration, the hydraulic fracturing calculation of the single-hole model was carried out at the same time. A small hole with a diameter of 30 mm was excavated in the center of the model. The strength distribution, boundary conditions, and pore pressure increment of the meso-unit were consistent with those in the double-hole model. Before the 34 h step of the single-hole model, there is no damage within the element. When the calculation reaches the 39 h step (the rupture pressure is 29.5 MPa), the model becomes unstable (Figure 3). The fracture pressure of the single-hole model is greater than that of the double-hole model, which is more prone to instability and failure due to the influence of pore water pressure between pores.

In the case of increasing water pressure, there is still an important characteristic pressure, namely, the initiation pressure, before the model reaches the rupture pressure. Crack initiation pressure corresponds to the inflection point from the straight line to the curve on the pressure-time curve. Andreev et al. [25] determined the crack initiation stress by observing the expansion of pore size in different directions. Cracking initiation stress can be monitored by acoustic emission (AE) [23]. Some scholars also believe that the fracture initiation stress corresponds to the pressure value when the rate of increase of borehole water pressure reaches the maximum value. Here, the aperture expansion and acoustic emission monitoring are used to determine the value. Figures 4 and 5 show the relationship between acoustic emission and aperture increment and pressure in the process of increasing borehole pressure. Due to the influence of the pore water pressure between the two pores, the initiation pressure of the single-hole model is greater than that of the double-hole model. At the same time, the crack initiation pressure of hole 1 is less than that of hole 2 due to the influence of material heterogeneity.

3.2. Effect of Rock Homogeneity on Hydraulic Fracture Propagation

Figure 6 shows the calculation result of case 2. The values of homogeneity were 1.5, 3, 5, and 1 000, the stress ratio was K = 1, and the double-hole spacing was 130 mm. The increment of water pressure in the hole is the same as in case 1. For strong heterogeneity of rock material (m = 1.5), the expansion path of crack is the longest, the cracks at both sides of hole 1 and 2 are asymmetric, and a large number of branch cracks form at the hydraulic fracture tip. This eventually leads to the formation of two main hydraulic fractures in a horizontal well. With the decrease in heterogeneity (i.e., m value increases), the number of microcracks distributed around the hole becomes less and less. For homogeneous rock materials (m = 1,000), almost no microcracks were observed, and the damage and fracture were confined to two hydraulic fractures. The hydraulic fractures extended little and were flat and smooth, developing in a symmetrical shape.

In conclusion, due to the existence of rock heterogeneity, hydraulic fractures always choose the path of least resistance to propagate, which depends on the statistical distribution of strength characteristics, resulting in irregularities of hydraulic fracture path. According to the theory of fracture mechanics, it is considered that a large number of microcracks are developed at the crack tip for heterogeneous rock. For homogeneous materials, the tip forms a plastic zone, which is again confirmed by the calculation results in Figure 6.

Driven by the local pore water pressure at the tip, the hydraulic fractures continue to initiate, expand, and connect. When the rock homogeneity is different, the distribution of pore water pressure along the horizontal line through the double-hole core is shown in Figure 7. The pore water pressure is the largest around the two holes and decreases when it is far away from the hole edge. Moreover, it can be found that the decrease in pore water pressure between holes is less than that on both sides of hole 1 and hole 2, and the higher pore water pressure is distributed in the middle of the double holes, so the propagation of hydraulic fracture between holes is stronger than that on the left side of hole 1 and right side of hole 2.

The distribution of pore water pressure is independent of rock homogeneity because the water pressure is pumped into the inner wall of the hole in all the calculation models. The pore pressure distribution of several homogeneity models is basically the same. Figure 8 shows the distribution of effective stress along with the cross-hole core profile when the homogeneity coefficient is 1.5, 5, and 1,000. It can be seen that when the heterogeneity of rock is strong (m = 1.5), due to the great mechanical difference of mesoscopic elements in the model, the stress curve fluctuations and jumps are more obvious. With the increase in homogeneity, the fluctuation decreases gradually. When the homogeneity is 1,000, the distribution of the effective maximum and minimum stresses is smooth and linear. In addition, the stress of the heterogeneous rock fluctuates on both sides of the stress curve of the homogeneous sample, and the average stress of the heterogeneous rock is basically the same as that of the homogeneous sample. It can be seen from Figure 8 that the greater the degree of homogeneity, the closer the effective minimum principal stress around the two holes, and the effective stress of homogeneous material is completely symmetrical to the two injection holes. In the center of the connection between two holes, the local effective minimum principal stress is the largest and gradually decreases as it approaches the hole wall, which is also the reason why the tensile hydraulic fracture runs through the area between holes and is parallel to the horizontal stress.

3.3. Effect of Stress Ratio on Hydraulic Fracture Propagation

Figure 9 shows the calculation results of Scenario 3. The homogeneity of the model is 2, and the spacing between double holes is 130 mm. The evolution characteristics of hydraulic fractures are simulated when K is 0.5, 0.8, 1.5, and 2.0. When K = 0.5, the stress difference is the largest, and the vertical tensile crack develops along the direction of σV. Due to the strong heterogeneity of this model (m = 2), the crack development of the two holes is not symmetrical, the crack extension length of the left hole is larger than that of the right hole, and there is no horizontal crack at this time. When K = 0.8, horizontal cracks gradually appeared in the right hole, and the two vertical cracks on the outside of the right hole deviated upward to the right, showing a tendency to develop in the horizontal direction. The results show that the development of hydraulic fracture is not only affected by the natural stress field but also restricted by the stress field and hole spacing.

When K < 1, the natural stress field and the distance between holes have similar effects on the crack propagation, and the crack propagation direction may be ultimately affected by the homogeneity of the rock medium. When K = 1.5, both horizontal crack and inclined crack appear near hole 1. However, due to the influence of hole 2, the crack finally propagates along the horizontal direction, and the crack propagation along the horizontal direction is dominant. When K = 2.0, the maximum principal stress is in the horizontal direction, and only the horizontal hydraulic fracture initiation, propagation, and connection occur in hole 1 and hole 2. This situation also shows that in the strong anisotropic stress field, the stress field controls the growth of the hydraulic fracture, and the interaction between the two holes cannot be ignored.

3.4. Effect of Hole Spacing on Hydraulic Fracture Propagation

Figure 10 shows the calculation results of Scenario 4. The homogeneity of the model m = 2 and stress ratio K = 1 were calculated to investigate the influence of different double-hole spacing on crack development. The initial pore pressure of the two holes is 10 MPa, and the pore pressure increment is 0.5 MPa. When the distance between the two holes is relatively close (L = 60 mm), three hydraulic fractures develop around hole 1 and hole 2, respectively, and the direction of hydraulic fractures is neither horizontal nor vertical, which indicates that the initiation and propagation of hydraulic fractures are affected by stress field and pore pressure at the same time when the two holes are close together.

With the increase in hole spacing, the interaction between the two holes becomes weaker. When the hole spacing is 90 mm, two vertical cracks appear in hole 1, while three hydraulic fractures appear in hole 2. Under the influence of hole 1, one horizontal crack develops. When the hole spacing increases to 130 mm, two vertical cracks develop in both holes, and no horizontal cracks appear. When the hole spacing reaches a maximum of 160 mm, only two vertical cracks develop in hole 2, and there is no hydraulic fracture initiation around hole 1. The results show that with the increase in hole spacing, the influence on the hydraulic fracture becomes weak. When the hole spacing is greater than 5 times the diameter, the development of hydraulic fracture between holes will not be affected by the pore water pressure of double holes. Geertsma and Detournay believed that the increase in compressive stress field would increase the pore water pressure and thus resist the generation of hydraulic fracture. It can be seen from the calculation results in this section that for the double-hole model, the factors that resist the formation of hydraulic fracture also include the distance between the two holes.

3.5. Effect of Initial Hydraulic Pressure on Hydraulic Fracture Propagation

Figure 11 plots the calculation results of case 5. The model considers the effect of pore pressure gradient on crack growth, where the homogeneity m = 2 and the stress ratio K = 1. The initial pore pressure of hole 1 is 10 MPa, and the initial pore pressure of hole 2 is 13 MPa, 15 MPa, and 17 MPa, respectively. The pressure increment in the two holes is 0.5 MPa, so the pressure difference between the two holes is guaranteed to be equal in each calculation time step. The results show that the crack starts from the hole with high pore pressure (hole 2) and propagates to the hole with low pore pressure (hole 1) with an increase in pore pressure difference. When the model fracture stress is reached, no crack initiation is observed around hole 1 due to the small initial pressure. The simulation results are consistent with the experimental results of Bruno et al., indicating the tendency of hydraulic fracture propagation from high pore pressure to low pore pressure under hydraulic gradient conditions.

The relationship between the critical pressure and the initial pressure difference of the double hole is shown in Figure 12. Both the crack initiation pressure and the fracture pressure in the model have a good linear relationship with the pressure difference. The fracture pressure and initiation pressure of the model decrease gradually with the increase in pressure difference.

In order to further explore the influence of the effect of the global pore water pressure on the crack evolution of the model, the distribution of effective stress and pore water pressure in the cross-hole profile is shown in Figures 13 and 14, respectively. Due to the existence of the pressure gradient in the double hole, the distribution of the minimum principal stress in the two holes is not symmetrical. In the middle region of the double hole, the accumulation rate of the effective minimum principal stress in hole 2 is higher than that in hole 1. When the tensile stress exceeds rock tensile strength the hydraulic fracture initiates at hole 2 and expands toward hole 1. It can also be seen from the distribution curve of pore water pressure that with the increase in pore pressure gradient, the pore water pressure around pore 2 is obviously larger than that around pore 1, and it gradually decreases with the increase in the distance from the pore 2. The existence of high pore pressure zone between the two holes leads to the development of hydraulic fractures more dramatically than that on both sides of the pore.

4. Discussion

In this study, the evolution of hydraulic fracture in a two-hole model was numerically investigated using a coupled flow-stress-damage model. The effects of rock homogeneity, stress ratio, hole spacing, and hydraulic pressure gradient were investigated. The characteristics of hydraulic fracture propagation in the two-hole model are the classical issue of stress shadow. The results of this study are different from that of pre-existing studies, and the homogeneity of rock was first considered in the two-hole model fracturing. Although the stress disturbance on rock fracturing has been revealed [18, 26], however, a numerical model considering disturbed stress on hydraulic fracture propagation is not common. It is accepted by many scholars that the shale formation is disturbed by the frequent injection of high-pressure fluids; this is to say, the pore pressure in the shale formation is impacted by the dynamic stress disturbance. In this study, the disturbance caused by the stress ratio is considered, and in further studies, the influence of disturbed loads, such as cyclic injection loads and vibration loads, should be added to the model.

5. Conclusions

In this study, the evolution of hydraulic fractures under local and global pore water pressures in a two-hole model is investigated numerically. Five kinds of numerical simulators are designed, considering the factors such as rock heterogeneity, natural stress ratio, double-hole spacing, and pore pressure gradient. The fracture evolution morphology, stress distribution, pore water pressure distribution, and critical pressure fraction of the five models were analyzed, and the main conclusions are summarized as below:(1)Under the condition of constant initial water pressure, compared with the single-hole model, the two-hole model is more prone to instability failure due to the influence of pore water pressure. Under the action of local pore water pressure at the crack tip, the interaction between the two holes is more obvious, and the fracture propagation length between the two holes is larger than that outside the two holes. When the pore pressure gradient is considered, the hydraulic fracture starts in the well with high initial pore pressure, and the crack mainly propagates in the middle region of the two holes under the influence of the “attraction effect” of the two holes.(2)Due to the differential heterogeneity of the rock model, the effective stress curves of the two-hole model have different shapes. With the increase in homogeneity, the stress curve becomes flatter, continuous, and smooth. Under the action of local pore pressure, the number and morphology of microcracks at the hydraulic fracture tip decrease with the increase in rock homogeneity. When the rock is homogeneous, the plastic zone is developed at the fracture tip, and no microcracks are observed.(3)It is shown that the overall stress ratio of different natural pore pressure impacts the performance in the formation of crack development. The number of cracks decreases with the increase in pressure ratio, and fracture propagation path also occurs. The numerical results suggest that the increase in pore water pressure caused by the in situ stress field in an overall scope will resist the propagation of hydraulic fracture.(4)For the double-hole model, the double-hole spacing is another important factor affecting the crack evolution. The results show that the pore water pressure between the two pores has an obvious effect on the fracture development when the distance is less than 5 times the hole diameter, and the crack is perpendicular to the maximum principal stress. When the distance is greater than 5 times the diameter, the pore water pressure between the two pores has almost no effect on the fracture evolution. However, it should be pointed out that the determination of the critical distance of double-hole influence may also be related to rock homogeneity. The 5-time hole spacing is obtained when rock homogeneity is 2, and the relationship between the homogeneity and the critical influence distance should be studied in further studies.

Data Availability

The experimental data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest.

Acknowledgments

The authors would like to thank Mechsoft Technology (Dalian) Co., Ltd. to provide the RPFA-Flow CODE. This research was funded by the Beijing Science and Technology Planning Project (Z181100005118012) and the China Coal Science and Industry Group Science and Technology Innovation Venture Capital Special Key Project (2018-2-ZD007).