The bedding layers in shale reservoirs are highly developed with numerous bedding joints. In the hydraulic fracturing process, the propagation of hydraulic fractures during shale formation significantly affects the efficient development of reservoirs. In this study, we employed the pulsating cyclic loading technique to develop a model of the interaction between shale hydraulic fractures and bedding layer; in addition, we used a numerical model optimization method for the zero-thickness cohesive element to compare the hydraulic fracture and bedding layer propagation behavior under static and pulsating loading. We studied the propagation behavior of pulsating fracturing bedding layers at different frequencies, amplitudes, and bedding angles. The results revealed that static load fracturing is more conducive to the propagation of primary fracture under the same conditions, and that pulsation fracturing is better for the propagation of bedding layer. The higher frequency and amplitude of the seismic source in the pulsation fracturing process are more favorable to crossing the bedding layers. When the bedding angle is approximately 10°, turning the main fracture along the bedding direction is challenging. The hydraulic fracture along bedding layer propagation was achieved when the bedding angle was approximately 30° and 45°. When the bedding angle was approximately 60° or more, the fracture extended only along one side of the bedding layer. The study results provide guidance for the formation of fracture networks in the field of fracturing construction of shale reservoirs, thus improving the production of shale reservoirs.
Shale gas is widely distributed and has excellent potential for extraction. Efficient extraction of shale gas can effectively relieve the pressure of declining conventional energy sources. Hydraulic fracturing is the primary means of shale gas extraction. Artificial fractures created by hydraulic treatment of shale reservoirs communicate with bedding layers and natural fractures as an effective production method. However, the variability of geological conditions in each reservoir and the strong inhomogeneity of the shale make the fracture pattern of hydraulic fracturing very complex [1–8]. The researchers concluded that the fracture pattern of hydraulic fracturing in shale reservoirs is controlled by natural factors such as natural fracture and bedding distribution, rock characteristics, tectonic stress, and engineering factors such as fracturing fluid viscosity and pumping speed[9–11]. Pulsed fracturing requires much less pressure than standard static load hydraulic fracturing. It does not cause local stress concentration; also, it is difficult for conventional hydraulic fracturing to make the fracturing fluid enter these microstructures. This paper carries out the investigation of the propagation law of hydraulic fractures interact with bedding layer under pulse fracturing.
Hou et al. , Zhang et al. , and Wang  showed that higher horizontal stress differences are favorable for the formation of more extensive prominent hydraulic fractures, but not for the formation of complex fracture networks in the vicinity of the wellbore. However, Jiang et al.  and Guo et al.  suggested that higher stress differences may lead to further extension of the main hydraulic fractures and connectivity with more natural fractures, resulting in a relatively more complex fracture network. Guo et al.  investigates the effects of subsequent injection in parent-well with legacy production on interwell interference using a data set from Eagle Ford Shale. A numerical modeling work flow is presented for the characterization of poroelastic behaviors of multiphase-fluid diffusivity and rock deformation using the finite-element method and multifracture propagation using the displacement discontinuity method. It solves for the spatial-temporal evolutions of pore pressure and in situ stress because of parent-well production and injection and models the fracture propagation during infill-well completion on the basis of updated heterogeneous in situ stresses. Guo et al. proposes a reservoir, geomechanics, fracturing modeling work flow for understanding the interference mechanism and quantifying effects of parent-well fracture geometry, differential stress, and the design of infill-well completion on interwell fracturing interference. Zhang et al.  determined the actual hydraulic fracture morphology within the shale by large-size true triaxial hydraulic fracture simulation combined with industrial high-energy computed tomography (CT). It appears that when the horizontal ground stress difference is too high or too low (such as greater than 12 MPa or less than 3 MPa), it is not conducive to the generation of the volumetric fracture network. In addition, the horizontal in situ stress difference required for fracturing to form a complex fracture system is different under different physical conditions of shale reservoirs. Low viscosity fracturing fluids are believed to activate natural discontinuities and facilitate complex fracture networks, while high viscosity fracturing fluids facilitate the formation of a single primary hydraulic fracture .
Nowadays, most of them use pulsation fracturing technology to CBM development to improve the permeability of CBM by pulsation fracturing and achieve efficient outcomes [19–23]. Among previous researches, Lin et al.  proposed the technique of high-pressure pulsation hydraulic fracturing of coal seams to unload pressure and increase permeability. It concluded that the alternating stress could induce fatigue damage in coal rocks and enhance permeability at the same time it confirmed by field tests that the extraction volume fraction increased by 264.7% on average after pulsation fracturing. Zhai et al.  concluded that the compression-expansion-compression effect caused by pulsation fracturing would lead to fracture penetration within the coal rock. Li et al.  concluded that pulsating stress waves’ reflection, superposition, and energy accumulation are the key to pulsating fracturing to increase penetration. They thus can produce better fracturing results than conventional fracturing with smaller pulsating pressure. Li et al.  found that low frequency and low pressure can fully develop microfractures. At the same time, high frequency and high pressure can promote the rapid propagation of the primary fractures, thus proposing the “dual-frequency-double-pressure” fracturing process. Lin et al.  studied the mechanism of high-pressure pulsation breaking; Li et al.  investigated the effect of pulsation parameters on the fracturing effect of pulsation hydraulic fracturing by indoor experiments. However, most of the current theoretical studies are based on the analytical solution of the fluctuation equation, which cannot simulate the actual stratum. At the same time, the indoor tests often cannot truly reflect the infinite stratum stress wave propagation process due to the limitation of test equipment and scale. The large-scale field tests can only be used as a verification means which is challenging to meet scientific research needs. In contrast, numerical methods can avoid tedious theoretical derivations while considering complex geological features and ground stress environments and are highly adaptable.
It can be seen from the above papers that the development process of shale reservoir and coal reservoir is aimed at forming complex fracture network. In this study, the pulsating cyclic loading technique was employed for the hydraulic fracturing treatment of shale reservoirs. In addition, the effects of various factors on the propagation of bedding layers in shale reservoirs under pulsating cyclic loading were investigated. Moreover, the propagation law of the interaction of hydraulic fractures with bedding layers in shale under static fracturing and pulsating load fracturing was compared. The findings presented in this study will provide guidance in the field of fracturing construction of shale reservoirs.
2. Methods and Models
2.1. Cohesive Element Crack Expansion Principle and Model Optimization Method
The schematic diagram of the cohesive element in the two-dimensional plane is shown in Figure 1, which has six node elements, where nodes 1, 2, 5, and 6 form the upper and lower edges of the cohesive element, and nodes 3 and 4 are the middle nodes of the element. The tangential flow of the fluid is described at the intermediate nodes, while the normal flow is presented at the upper and lower edges of the element.
2.1.1. Cohesive Element in ABAQUS Simulates the Principle of Crack Propagation
2.1.2. Damage Initiation and Evolutionary Guidelines
2.1.3. Fluid Flow Criteria
The flow of fracturing fluid in the cohesive element can be divided into two types. One is the tangential flow within the cohesive element, and the other is the normal flow perpendicular to the upper and lower surfaces of the cohesive element. The difference is that the XFEM method cracks do not depend on the grid elements and flow directly through the grid elements. In contrast, the cohesive element method relies heavily on the grid and flow between the element point and lines, as shown in Figure 2. In the figure, is the fluid injection per unit time, is the fluid outflow per unit time, is the difference between the two, the actual fluid injection acting on the fracture unit, is the fracture top opening velocity, and is the fracture bottom opening velocity.
2.1.4. Model Optimization Method for Embedding Zero-Thickness Cohesive Cohesion Elements 
By embedding the whole model with zero-thickness cohesive cohesion element in the entire domain, it can solve the problem that the cohesive elements between the nodes of the cohesive element cannot be assigned in the process of traditional hydraulic fracture propagation in ABAQUS. The cohesion between rocks in the actual formation cannot be ignored; so, the model is closer to the actual formation after embedding the zero-thickness cohesion element, making the numerical model more accurate. The model optimization method can be specifically divided into the following three steps.
Set up the model mesh type by ABAQUS/CAE, perform meshing, and output the inp file. Get the node information of the solid element from the finite element model inp file by adding a new node and then obtaining solid elements that do not share a common point with each other. The number of solid elements and the order of element nodes are not changed
Insert the cohesive elements into the initial mesh in the input inp file and defining the zero-thickness cohesive element as follows: based on the previous step, element pairs with overlapping surfaces are searched. The zero-thickness cohesive elements are defined at the spatial location of the overlapping surfaces. As shown in Figures 3(a) and 3(b), the red numbers are the element nodes, and the black numbers are the element numbers
The node and element information obtained above is written to an inp file and imported into the ABAQUS user interface to obtain the computational model
2.2. Modeling Shale Bedding Seam Extension under Pulsating Load
Based on the pulsation hydraulic fracturing technology proposed by predecessors in coalbed methane development and research , in this paper, the pulsating fracturing fluctuation function is imported to establish a model for shale bedding fracture extension under pulsating load. It is established based on the following conditions.
Assuming that shale is an anisotropic elastic medium, and many bedding surfaces and bedding fractures are distributed in the shale reservoir. However, the width is much smaller than the stress wavelength caused by pulsation fracturing, and the stress wave will be diffracted when it travels to the bedding layer. The overall stress distribution will not be seriously changed; it can be approximated as an isotropic elastic medium
The displacement in the direction of gravity does not change with the change of azimuth; that is, the stress stays in a plane strain state. The model plane is the vertical two-dimensional cross-section of the formation in horizontal well drilling
In the development of horizontal wells in shale reservoirs, the wellbore size is much smaller than the entire shale reservoir; also, the stress difference along the way is negligible. Therefore, it can be simplified to a point source model with the same phase as shown in Figure 4, where the model midpoint is the wellbore fracture injection point, m; is the perforation length, and the half slit length is the perforation depth, m; is the longitudinal separation distance between shale bedding fractures and the injection point, m; is the angle between the bedding fractures and the main perforation fracture °, and the bedding angle is 90-θ; σH is the maximum principal stress, MPa; σh is the minimum principal stress, MPa
The location of the fractures in the model is grid encrypted, and the cohesive element with zero thickness is embedded in the grid to obtain the model grid schematic is shown in Figure 5. The local grid encryption of the model can improve the efficiency of the calculation and make the crack extension orientation calculation more accurate. At the same time, the embedded zero-thickness cohesive cohesion element makes the simulation results more realistic, where the red diagonal position in the figure is the bedding layers, and the red dot in the center is the injection point.
2.3. Pulsation Fracturing Frequency Amplitude Fluctuation Formula Derivation
3. Model Validation
This paper numerically simulates the vertical expansion behavior of shale bedding layers in a horizontal well fracturing section with a formation of about 1000 m below ground, with fixed boundary conditions and a side length of 20 m, subject to the maximum principal stress and horizontal principal stresses in the formation; for the equilibrium ground stress state, the reliability of the model can be verified by physical simulation experiments and dynamic and mechanical static examples. The relevant physical parameters are shown in Table 1.
3.1. Experimental Verification
The numerical model was verified by hydraulic fracturing experiments. True triaxial loading was used in the experiment, and the loading process is shown in Figure 7, where is the vertical stress, and are horizontal stress, and α is the bedding dip angle. Grouping of similar specimen working conditions is shown in Table 1. The tensile strength of similar shale specimens was determined based on the mechanical parameters of the specimen of Changning shale in Sichuan province. The matrix tensile strength was set as 10.57 MPa, and the bedding tensile strength was 2.89 MPa and 5.36 MPa, respectively.
The experimental results are shown in Figure 8. Through the hydraulic fracturing experiment, it is found that the tensile strength of the bedding layer is weak when the fracture extends along the bedding layers. When the tensile strength of the bedding layer is strong, it will run through the whole fracture. Pulse fracturing is mainly manifested in the process of the pulsating load effect of bedding layer, dimensional bedding plane set weak bedding surface, and the tensile strength of 3 MPa.
3.2. Static Validation
The comparison of the theoretical and simulated results for σxx along the and symmetric axes is calculated and shown in Figure 9. It can be seen that the rest of the region is very close to the theoretical value (maximum error of 3.9%) except for 2 to 3 nodes near the crack tip where the calculated results differ due to the end effect (maximum error of 40%). In fact, the stress singularities at the crack tip cannot be eliminated . Still, the stress singularities can be restricted to a small area by changing the grid step without affecting the overall stress distribution. In view of the fact that the focus of this study is on the propagation and distribution of stress in the formation, not limited to the tip of the fracture, hence, the static response of this model is sufficiently accurate.
3.3. Dynamic Verification
From the above equations, the radial stress values in the fluctuation concentration region at the source point can be derived and compared with the stress values calculated by the numerical model simulation in this paper as shown in Figure 11. As can be seen from the figure, the waveform, amplitude (maximum error of 13.2%), period (maximum error of 4.2%) of the analytic solution, and the simulated solution are basically consistent except for the oscillation of the numerical solution at the initial moment. Only the phase is delayed (maximum error of 38%), and this delay may be caused by the delayed misalignment from the model injection rate oscillation to the radial stress transition. The main focus of this paper is on the frequency and amplitude of the vibration; so, the kinetic response of this model is considered accurate enough.
The comparison of injection pressure with time under pulsating load and static load is shown in Figure 12. From the figure, it can be seen that the injection pressure under pulsating load is basically the same as the trend of pressure change in the fracture. The numerical simulation of the pulsating load injection pressure fluctuates between 12.5 and 44.3 MPa cycle loading. However, the pressure change in the simulation is not obviously same to the perfect state. The pulsating load injection pressure up and down peak is not exactly the same which has 2.8 MPa. The static loading fracture also has certain fluctuations, but the fluctuations are small which is only 0.6 MPa.
4. Numerical Simulation and Discussion
Since the tensile strength of shale is much lower than the applied compressive strength, the brittleness increases with increasing loading. As such, shale failures are considered brittle fractures. Under the action of periodic pulsating stress, when the absolute value of the tension within the fracture is greater than the tensile strength of shale (set to 6 MPa in this study), the fracture undergoes tensile damage and continues to expand with the continuous injection of fracturing fluid. By numerically simulating different stresses (0–6 MPa), frequencies (0–30 Hz), and amplitudes (0–0.006 m3/s), the interaction of shale fractures with the bedding layer was studied.
4.1. Comparison of Static Load Fracturing and Pulsating Fracturing
Figure 13 shows the comparison of dynamic and static fracturing effects under the conditions of maximum principal stress in the formation, minimum principal stress in the formation, injection time , and , and the stress unit in the figure is Pa, where Figure 13(a) shows the simulation results of fracturing and change in stress during fracture propagation under static fracturing with injection rate . Figure 13(b) shows the simulation results of the change in stress during fracture propagation under pulsation fracturing with initial injection rate , frequency , and amplitude .
As can be seen from Figure 13, the main fracture directly penetrated the two bedding layers under static load, and only extremely small cracks formed at the bedding layers; moreover, the propagation of fractures was minimal. In contrast, with the continuous propagation of the fractures under pulsation fracturing, the fractures expanded along the bedding layer when the main crack intersected with the first bedding layer, forming an inclined -shaped crack. However, the fracture direction on one side of the bedding layer was at an angle of 135° with the main fracture propagation, which is almost opposite to the primary fracture direction; thus, fracture propagation is slower, and the hydraulic fractures formed are narrower compared to the fractures formed on the other side of the bedding layer In addition, it can be seen from Figure 14 that under pulsation fracturing, the primary fracture at the bedding layers showed a trend of fine propagation at 10 s, while the main fracture expanded to form a fine fracture also contracted at 18 s. Thus, compared with static fracturing, pulsation fracturing is more favorable for the propagation of bedding layers than the main fractures.
The simulation results revealed that turning the main fracture after contacting the bedding layers is more difficult under static load than under pulsation. In the simulation, the primary fracture did not propagate in the bedding layers. In contrast, hydraulic fractures were more likely to propagate along with the bedding layers under pulsation fracturing even if the fractures propagated in almost opposite directions resulted in narrower hydraulic fractures. Constructing complex fracture networks in shale reservoirs is easier under pulsation fracturing and thus more conducive to the effective development of reservoirs. The fracturing method should be selected according to the demand, or static fracturing should be performed first to propagate the primary fracture fully. Subsequently, pulsating fracturing should be employed to form a more complex fracture network to develop the reservoir efficiently.
4.2. Effect of Frequency
Figure 15 shows the fracture propagation under the conditions of maximum principal stress in the formation, the minimum principal stress in the formation, the injection time , conditions, the frequencies , 10 Hz, 20 Hz, and 30 Hz, and injection rate , and the stress unit in the figure is Pa. Simulation results were obtained for the change in stress during fracture propagation under pulsation fracturing with amplitude .
As can be seen from Figure 15, after the main crack propagated to the first bedding layer at a frequency of 1 Hz, it propagated along one side of the bedding layer. Its direction was at an angle of 45° with hydraulic fracture propagation, while there was no expansion on the other side of the bedding layer at an angle of 135°. At 10 Hz, after theprimary fracture intersected with the first bedding layer, the fracture propagated along the bedding layer and formed a tilted -shaped fracture. The bedding layer on the other side at an angle of 135° with the expansion direction of the main crack also created a narrow hydraulic fracture. At 20 Hz, after the fractures propagated to form inclined-shaped cracks, longer and broader hydraulic fractures also developed on the other side with an angle of 135° from the primary fracture propagation direction. At 30 Hz, with the continuous injection of fracturing fluid, the main fracture also formed an inclined -shaped fracture at the bedding layer. An obvious hydraulic fracture formed on the other side at an angle of 135° from the main fracture propagation direction. Although it was still not as thorough as the extension of the bedding layer with the angle of 45°, its fracturing effect was extremely pronounced.
From the simulation results for pulsation fracturing, it can be seen that hydraulic fractures are more likely to propagate along with the bedding layer when other factors are identical. The higher the pulsation frequency, the better the fracture propagation effect, making the shale reservoir more likely to form a complex fracture network and more conducive to the effective development of the reservoir.
4.3. Effect of Amplitude
Figure 16 shows the simulation results of the change in stress under pulsation fracturing at maximum principal stress , minimum principal stress , injection time , , frequency , injection rate , and amplitude , 0.001m3/s, 0.002m3/s, and 0.003m3/s, and the stress unit in the figure is Pa, respectively. The simulation results of stress change during fracture propagation under pulsation fracturing were plotted.
It can be seen from Figure 16 that when the amplitude was 0.0005 m3/s, the primary fracture propagated and interacted with the first bedding layer. The fracture propagated along the side of the bedding layer, but the fracture was narrower than the one in Figure 15(a). Its direction was at an angle of 45° with hydraulic fracture propagation. At the amplitude of 0.001 m3/s, after the main fracture intersects with the first bedding layer, the crack expands along the bedding layer and forms an inclined -shaped fracture. The bedding layer on the other side with an angle of 135° with the direction of propagation fracture also generates a narrow hydraulic fracture. At the amplitude of 0.002 m3/s, after fracture expansion resulted in an inclined -shaped fracture, a wider hydraulic fracture formed on the side of the bedding layer fracture at an angle of 135° with the main fracture expansion direction. At the amplitude of 0.003 m3/s, with the continuous pulsation injection fracturing fluid, fracture propagation resulted in a wider and longer hydraulic fracture on the other side of the bedding layer fracture at an angle of 135° with the main fracture propagation direction. Although it was not as thorough as the fracture propagation on the side of the bedding layer fracture at an angle of 45° with the main fracture expansion direction, its fracturing effect was extremely obvious.
The simulation results show that when other factors remain unchanged in pulsation fracturing, the main fracture propagates along the bedding layer. The higher the amplitude of pulsation fracturing, the greater the effect of bedding layers in fracture propagation. The higher the amplitude, the more conducive the shale reservoir to form a complex fracture network under pulsation fracturing, which enables the efficient development of shale reservoirs.
4.4. Effect of Bedding Angle
Figure 17 shows the simulation results of stress change during fracture expansion under pulsation fracturing with the maximum principal stress in the formation, minimum principal stress in the formation, injection time , , frequency , injection rate , and amplitude of 0.001 m3/s, and the stress unit in the figure is Pa, respectively. The simulated results of stress change during fracture expansion under pulsation fracturing when the bedding angles are 10°, 30°, 45° and 60°, corresponding to θ equal to 80°, 60°, 45°, and 30°, respectively, were plotted.
From Figure 17, it can be seen that when the bedding angle was 10°, the angle between the main fracture propagation direction and the bedding layer on both sides was more extensive. Turning the fracture direction was difficult, and the main fracture directly penetrated the bedding layer and did not extend to the bedding layer. At the bedding angle of 30°, an inclined -shaped hydraulic fracture is formed with the continuous pulsating injection of fracturing fluid. A wider and longer hydraulic fracture developed on the other side of the bedding layer at an angle of 135° with the primary fracture propagation direction. The bedding layer at an angle of 45° with the main fracture propagation direction propagated almost to the same extent. When the bedding angle is 45°, after the primary fracture intersected with the first bedding layer, the crack extended along the bedding fracture and formed an inclined -shaped fracture. The bedding layer on the other side with an angle of 135° with the main fracture propagation direction also created a narrow hydraulic fracture. When the bedding layer angle was 60°, the main fracture propagated to the first bedding layer and then propagated along the side of the bedding layer. The fracture was wider, and its direction was at an angle of 45°, while the fracture at the other side of the bedding layer at an angle of 135° did not propagate.
The simulation results revealed that when other factors remain unchanged, the main fracture expansion was difficult to turn along the bedding layer when the bedding angle was 10°; moreover, it was more difficult to form a complex fracture network by using the fracturing treatment. When the bedding angle was 30o and 45°, the main fracture expanded along the bedding layer. When it interacted with the bedding layer, hydraulic fractures formed on both sides of the bedding layer were more conducive for forming a complex fracture network during the fracturing process, enabling efficient reservoir development. When the bedding angle was 60°, the main fracture propagated only along one side of the bedding layer and the hydraulic fracture was wider. Thus, during the fracturing treatment, the larger the bedding angle, the easier the propagation of the primary fracture along the side, where the fracture trend of the bedding layers is similar to the propagation direction.
In this study, the model optimization method of zero-thickness cohesive element was used to import the pulsating fracture fluctuation function, the shale bedding fracture propagation model was established under pulsating load, two fracture loading methods (static fracturing and pulsating fracturing) were compared, and the fracturing effect were simulated and analyzed under static load. Moreover, the effect of different frequencies, amplitudes and bedding layer angles on the pulsating fracturing effect was investigated. The conclusions are as follows.
Turning the primary fracture after contacting the bedding layers under static load than under pulsation. The hydraulic fracture is more likely to expand along with the bedding layers under pulsation fracturing, and it is easier to form a complex fracture network in the shale reservoir during fracturing treatment, compared with static fracturing, pulsation fracturing is not conducive to the propagation of the primary fracture
Under the same conditions, hydraulic fractures are more likely to propagate along with the bedding layers under pulsation fracturing. The higher the frequency and amplitude of the pulsation source, the better the fracture propagation effect, thus, making it easier to form a complex fracture network in shale reservoirs
The primary fracture propagation was difficult along the bedding layers when the bedding angle was 10°. The main fracture propagated along the bedding layers when the bedding angles were 30° and 45°, and hydraulic fractures formed on both sides of the bedding layers, which is more conducive to the formation of complex fracture networks during the fracturing treatment. The main fracture propagated only along one side of the bedding layers when the bedding angle was 60°, and the hydraulic fractures formed were wider. When the bedding angle was larger, the primary fracture propagated along the side of the bedding layers, where the fracture trend was similar to the propagation direction
The data used to support the findings of this study are all shown in uploaded manuscript.
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
This work was supported by the Natural Science Foundation of Heilongjiang Province of China (YQ2021E005), National Foundation Cultivation Fund of Northeast Petroleum University (2021GPL-02), Shale Oil Research Project of Northeast Petroleum University (YYYZX202102), and the National Natural Science Foundation of China (No. 51774094).