This section is on the basis of 1# rock sample. The vertical stress which was exerted on the rock was 20 MPa, and the maximum and minimum horizontal stresses forced on the rock were 8 and 5 MPa, respectively. During the physical simulation process of 1# rock sample, due to that the horizontal stress exerted on the rock sample did not change along the vertical direction, the interlaminar geostress difference was not considered.

After the completion of the hydraulic fracturing experiment, 1# rock sample was taken out of the stress loading system and was split along the propagation direction of the hydraulic fracture so that the fracture face could be revealed. Because the fracturing fluid was colored into red before, the red area on the rock split face was the fracture face area (propagation area), just as shown in the Fig. (5). It can be found that the fracture face was round when the interlaminar geostress difference did not exist, and a symmetrical double-wing fracture whose propagation direction was perpendicular to the minimum horizontal stress σ_h was generated.

Fig. (5) Hydraulic fracture face after the hydraulic fracturing.

Different fracture propagation processes lead to different acoustic emission events [9C. Yashwanth, M. Camilo, S. Carl, and R. Chandra, "An experiment investigation into hydraulic fracture propagation under different applied sresses in tight sands using acostic emissions", J. Petrol. Sci. Eng., vol. 108, pp. 151-161, 2013.
[http://dx.doi.org/10.1016/j.petrol.2013.01.002] ]. Fig. (6) presents the distribution of acoustic emission points at different stages of physical simulation, which can present the dynamic propagation process of the fracture. The dimensions here is the dimension used to process the acoustic emission points. It is larger than the dimension mentioned in the text, because the acoustic emission points appear not only within the rock volume but out of the rock when some noises happened there.

Fig. (6) Acoustic emission points distribution at different stages. (a) 400s; (b) 960s; (c) 2200s; (d)3240s.

It can be found that more acoustic emission points generated with the continuance of the hydraulic fracturing simulation. When the physical simulation time was 400s, because of the relatively short time, acoustic emission points emerged intensively in the location of open hole section of wellbore, just as shown in the Fig. (6a), and it shows the initiation of the hydraulic fracture. With the extension of experimental time to 960 and 2200s, the amount of acoustic emission points intensively increased as shown in Fig. (6b and c), and they present that the fracture gradually propagated. With a further extension of experimental time to 3240s, large numbers of acoustic emission points emerged at the edge of the rock as shown in the Fig. (6d). It could be inferred that the fracture fluid flowed out of the rock sample and the fracture reached the edge of the rock. According to the distribution change of acoustic emission points, it can be found that the lateral and vertical propagation of the fracture occurred, and finally a fracture with round-shape face was formed. It is due to the reason that without the interlaminar geostress difference, no barrier effect existed in the vertical direction, and then the fracture morphology was only controlled by the three-dimension principles stress condition.

In the numerical simulation study, the first group simulation was conducted to investigate the fracture propagation behavior when the interlaminar geostress difference did not exist. In order to simulate the absence of interlaminar geostress difference condition, three layers in the numerical model (Fig. 3) were given a same Young’s modulus value (upper = middle = bottom = 15 GPa). Fig. (7) shows the result of the numerical simulation with presentation of the change in fracture face. It shows from the Fig. (7a) that when the simulation time was 0.5 ms, the fracture face was relatively small and presented a round shape. With the increase of the simulation time from 0.5ms to 1.2 and 1.6ms, the fracture face was obviously enlarged; however, both of their shapes (Fig. 7b and c) were also round despite that the fracture has already reached and crossed the boundaries of the middle layer. Fig. (7d) illustrates that with the extension of the simulation time to 2.5ms, the fracture further propagated into the upper and bottom layers and the fracture face also showed the round shape. The result of the numerical simulation matched very well with that of the physical simulation.

Fig. (7) Fracture propagation in the first group numerical simulation. (a) 0.5ms; (b)1.2ms;(c) 1.6ms; (d)2.5ms.

The 2# rock sample was used to conduct the physical simulation test when the interlaminar geostress difference was considered. The vertical stress exerted on the 2# rock sample was 20 MPa which was same with that of the 1# rock sample. The horizontal stress exerted on the rock sides changed along the vertical direction. The upper and bottom layers were exerted relatively high horizontal stresses (σ_H = 15 MPa, σ_h = 12 MPa), and the middle layer was exerted relatively low stresses (σ_H = 8 MPa, σ_h = 5 MPa).

After the physical simulation was performed, the rocksample was cut into slices along the horizontal direction parallel to the fracture face, and then the cross sections of hydraulic fracture in different locations were obtained, just as shown in the Fig. (8).

Fig. (8) Cross-section of rock sample after the hydraulic fracturing. (a) picture of rock slices and fracture cross-section after physical simulation; (b) schematic of rock slices and fracture morphology.

Fig. (8a) is the picture of rock slices, and the trace lines on different slices reflected the propagation size of the fracture. That is to say, each slice can be perceived to be one cross section of the hydraulic fracture. The length of the line equals the height of the fracture in the certain location of that slice. In order to describe the fracture more specifically, a schematic was drawn as shown in the Fig. (8b). The blue lines represent the cross sections of fracture in different propagation locations along the horizontal direction, and the red line represents the fracture edge. Note that the Fig. (8) just presents one half of the hydraulic fracture. It can be clearly seen that the larger the distance from the wellbore, the smaller the fracture height, and then an elliptical fracture face was finally achieved.

Fig. (9) illustrates the acoustic emission points distribution at several stages of the physical simulation of 2# rock sample. When the simulation time was 440s, points were confined to an area close to the open hole section of wellbore, as shown in the Fig. (9a), and the shape of the area was nearly round. It was due to that the fracture did not reach the boundary of the middle layer at that time. When the physical simulation time was extended to 920s, a large amount of acoustic emission points appeared and located near the boundary of the middle layer (Fig. 9b). It can be inferred that the fracture may reach the boundary. With the increase of the simulation time to 1800s, numbers of points occurred at the right and left parts of the middle layer; however, only a few points appeared at the upper and bottom layers (Fig. 9c). This phenomenon presents that the fracture mainly propagated in the middle layer, and its propagation along the vertical direction was limited by the upper and bottom layers with relatively high stress. From (Fig. 9b and c), more points emerge at the left part than at the right part, and it indicates that the fracture propagation may be not uniform. It may go to the left first then the right, which has no influence on the final shape of the fracture. Fig. (9d) shows the result that at the final time of 2800s, the vast majority of points appeared at the edges of middle layer, and it illustrates that the fracture reached the edge of the rock sample, and the fluid flowed out of the rock. It can be concluded from the above results that with the presence of the interlaminar geostress difference, the fracture height was obviously limited due to the barrier effect of the high stress barriers, and the final shape of the fracture face was elliptical.

Fig. (9) Acoustic emission points distribution at different simulation stages (2# rock). (a) 440s; (b) 920s; (c) 1800s; (d) 2800s.

In the numerical simulation of the fracture propagation in layers with different interlaminar geostresses, three layers were separated and given different Young’s modulus values (upper: 30GPa; middle: 15GPa; bottom: 30GPa) in the numerical model as shown in the Fig. (3). Similarly, the upper layer and bottom layer were given relative high value than that of the middle layer, and then the influence of geostress difference on the fracture propagation was simulated. The simulation result is shown in Fig. (10). It shows that when the simulation time was 1.2ms, the fracture expanded mainly within the middle layer, and the shape of fracture face was round as shown in Fig. (10a and b). The reason is that before the fracture reached the boundary of the middle layer, the stress difference had little influence on the fracture growth, just as shown in the section 4.1. With the increase of the simulation time to 1.6ms, the fracture propagated to the boundaries of the middle layer, and the fracture height was limited due to the high stress of the upper and bottom layers (Fig. 10c). The fracture tended to continually propagate along the horizontal direction. It also proved that the high stress barriers had an obviously barrier effect on the height of the fracture. With a further extension of simulation time from 1.6 to 2.5ms, despite a little degree of the fracture propagation occurred in the upper and bottom layers, the fracture mainly propagated along the horizontal direction in the middle layer due to the barrier effect, and the final shape of the fracture face presented to be elliptical, just as shown in the Fig. (10d). The result of the numerical simulation was consistent with that of the physical simulation in the section 4.2.1. When the height of the fracture is limited, the PKN model is suitable in simulation. However, the fracture height in PKN model does not change in the length direction, which is not consistent with the situation discussed here. According to Figs. (8 and 10), despite the limitation from upper and bottom layers, the fracture grew into these layers.

Fig. (10) Fracture propagation in the second group numerical simulation. (a) 0.5ms; (b) 1.2ms; (c) 1.6ms; (d) 2.5ms.

The consistency between the physical and numerical simulations verifies the feasibility of this numerical simulation method. Therefore, the numerical simulation method which based on the extended finite element method was used to further simulate the fracture propagation behavior when several special layers with different stiffness exist at the same time.

In this section, we simulated a case using the numerical simulation method to investigate the hydraulic fracture propagation behavior when several special layers with different stiffness existed simultaneously. A horizontal well model with a length of about 100 meters was built, as shown in Fig. (11). Three special layers numbered 1, 2 and 3 were placed, and they were given different Young’s modulus values, and the assigned values are shown in Table 2. Among the three layers, the Young’s modulus values of layer 1 and layer 3 were higher than that of layer 2. The layer 2 has a same distance to the wellbore with that of layer 1; however, the distance between layer 3 and wellbore is relatively long compared with that of layer 1 or 2. To form a vertical hydraulic fracture perpendicular to the horizontal wellbore, the minimum horizontal stress σ_h was exerted on the model along the direction which was paralleled to the horizontal wellbore; meanwhile, the maximum horizontal stress σ_H was forced on the model along the perpendicular direction, and a vertical stress was forced along the vertical direction. The parameters and their assigned values are presented in Table. 2.

Fig. (11) Horizontal well model with a hydraulic fracture approaching the pre-existing different layers.

Table 2
Parameters of numerical simulation model for case study.

Fig. (12) Propagation behavior before and after the hydraulic fracture reached the special layers. (a) hydraulic fracture began to extend; (b) hydraulic fracture reached layers 1 and 2; (c) hydraulic fracture reached layer 3.

The result of the numerical simulation of the fracture propagation when several special layers with different stiffness exist is shown in Fig. (12). The result was calculated without remeshing procedure that was an advantage of the simulation method. This figure presents the change of fracture vertical cross-section during propagation at different simulation stages. This view orientation is different from the Figs. (7 and 10). At the first stage (Fig. 12a), the hydraulic fracture began to propagate from the wellbore. At the second simulation stage as shown in Fig. (12b), the fracture crossed the layer 2 and stopped at the layer 1, because the distance between layer 1 and wellbore equals that between layer 2 and wellbore. The crossing phenomenon of fracture is mainly due to the relatively weak stiffness of the layer 2. That is to say, a layer with a weak stiffness may not prevent the fracture propagation. In contrast, the fracture stopped at the layer 1 which has a stronger stiffness. Similar result was obtained at the third simulation stage when the hydraulic fracture reached the layer 3 which also has a relative high stiffness (Fig. 12c). After passing through layer 2, the hydraulic fracture kept growing along the direction toward layer 3 and finally stopped at there. Therefore, the fracture height was limited due to the barrier effect of layers 1 and 3. It can be also found that the hydraulic fracture became wider with the continuance of the simulation. It can be concluded from the above results that the fracture passed through the weak and stopped at the strong. It is due to the reason that the hydraulic fracture propagate when the pressure in the fracture is higher than the rock failure stiffness [20J. Oliver, and A. Huespe, "Continuum approach to material failure in strong discontinuity settings", Comput. Methods Appl. Mech. Eng., vol. 193, pp. 3195-3220, 2004.
[http://dx.doi.org/10.1016/j.cma.2003.07.013] ], and that the pressure is high enough to make the fracture cross the weak layer but not enough to cross the strong layer.