Next Article in Journal
A Numerical Investigation on Kick Control with the Displacement Kill Method during a Well Test in a Deep-Water Gas Reservoir: A Case Study
Next Article in Special Issue
Study on Near-Wellbore Fracture Initiation and Propagation with Fixed-Plane Perforation in Horizontal Well for Unconventional Reservoirs
Previous Article in Journal
Tight Oil Well Productivity Prediction Model Based on Neural Network
Previous Article in Special Issue
Study on the Influence of Perforating Parameters on the Flow Rate and Stress Distribution of Multi-Fracture Competitive Propagation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Use of Extended Finite Element Method to Characterize Stress Interference Caused by Nonuniform Stress Distribution during Hydraulic Fracturing

1
School of Petroleum Engineering, Yangtze University, Wuhan 434000, China
2
Research Institute of CNOOC Shenzhen Branch, Shenzhen 518068, China
3
Downhole Service Company, CNPC Chuanqing Drilling Engineering Company Limited, Chengdu 610051, China
*
Authors to whom correspondence should be addressed.
Processes 2024, 12(10), 2089; https://doi.org/10.3390/pr12102089
Submission received: 5 September 2024 / Revised: 15 September 2024 / Accepted: 17 September 2024 / Published: 26 September 2024

Abstract

:
Stress interference is the main factor affecting hydraulic fracture propagation during multi-well hydraulic fracturing; stress interference is influenced by fracture bending, fracture hits, and asymmetric fracture propagation. To investigate the role of stress interferences among hydraulic fractures with nonuniform stress distribution in an inhomogeneous formation, a hydromechanical coupling extended finite element method was adopted to investigate the fracturing paths that occurred during the first fracturing–fracturing fluid flowback–repeat fracturing process; the asymmetric fracturing that occurred at different child well locations was also studied. The results showed that the area affected by fracturing-induced stress formed a “butterfly type” area. For child wells located within the zone, stress interference resulted in asymmetric fracture propagation; meanwhile, for child wells located outside this zone, stress interference resulted in symmetric fracture geometry. The effect of stress interference on the asymmetry of child well fracture wings was found to be negatively correlated with the distance between the parent well and the child well.

1. Introduction

For low-permeability hydrocarbon reservoirs, hydraulic fracturing increases the contact area with a low-permeability matrix and improves oil and gas production efficiencies [1]. Currently, well pad hydraulic fracturing, which has parent and child wells in adjacent areas, is the main technology that is used to increase the utilization degree of reservoirs and to achieve greater economic benefits [2]. When a child well is successfully fractured, hydrocarbons can be efficiently produced in the infill area; but stress interference between adjacent hydraulic fractures is the main factor contributing to fracture bending, fracture hits, and asymmetric fracture propagation in child wells [3,4]. Therefore, it is crucial to study the parameters that affect the fracture propagation of child wells to support the production of reservoirs [5].
Normally, hydraulic fracturing procedures for parent wells are carried out prior to those for child wells; the induced stress and fluid pressure have a great effect on the magnitude and direction of the principal stresses near the fracture. For this reason, the initial stress state of adjacent child wells is considered to be a superposition of in situ stress and induced stress. During the fracturing, the overlapping initial stressed lead to the development of complex fractures near child wells. The complexity of child well fractures may cause inter-well interference, which can have a negative effect on production performance. In some cases, fracture hits induced by inter-well interference can result in abnormal pressure; this situation can further effect the cementing safety of production wells. Yaich et al. [6] have indicated that the infill fracturing of wells in the Marcellus shale reservoir caused serious interference for the production wells; the formation pressure of the fracturing wells significantly exceeded the original pressure in the work area, and there was a significant decrease in production after the fracturing at the infill wells [7]. Miller et al. [8] investigated the proportion of production wells among shale oil and gas reservoirs in the US that were affected by adjacent fracturing wells; the proportions of affected production wells was found to be as follows: 41% in Eagle Ford, 56% in Niobrara, and 64% in Woodford [9]. The research results obtained by Lindsay et al. indicated that, within a 1000 ft cluster space, 60–80% of fracturing wells can affect the productivity of adjacent production wells [10].
The asymmetric propagation of two wings along a child wellbore has been found to be the main cause of inter-well interference. Walser and Siddiqui [11] found that an asymmetric fracture could significantly reduce recoverable reserves and recovery; additionally, the properties of reservoirs, including reservoir permeability and contrasting permeability between layers, could also affect the long-term recovery of hydrocarbons [12]. Rezaei and Rafiee [13] introduced a novel approach for calculating dynamic fracture propagation and stress changes in depleted regions; they showed that, due to the asymmetry of the stress field, the fracture wing of the infilling well near the parent well changed to longitudinal; meanwhile, the other fracture wing, which was away from the parent well, remained lateral, resulting in highly asymmetric fracture geometry [14]. Based on an elegant microseismic trial conducted by the Athabasca Oil Corporation, Stephenson et al. revealed that, in a complex in situ stress field, the fracture could be asymmetric, and this fracture asymmetry is caused by the well azimuth and the state of in situ stress [15]. Hu et al. employed a 2D fracture model with microseismic observations to capture the fracture hit behavior between adjacent horizontal wells, which showed that stress barriers around horizontally fractured wells could result in asymmetric fracture wings in adjacent wells [16].
Hydraulic fracture models have been shown to be key in analyzing stress interference among hydraulic fractures. Many research efforts have been conducted in the past decade with a focus on quantitatively analyzing the interactions that occur between adjacent fractures [17,18]. Cheng et al. developed a 2D displacement discontinuous method (DDM) to investigate the impacts of the perforation cluster number and spacing on production performance; the results suggested that inadequate cluster spacing could lead to less effective or ineffective fractures, thus resulting in a lower gas rate and ultimate recovery [18,19]. Olson introduced a height-correction formula to calculate fracture morphology when multiple fractures propagate at subcritical speed; the fracture propagation velocity that occurs during stable growth is proportionate to the stress intensity factor at the fracture tips [20]. Wu et al. obtained a more accurate calculation result for multi-fracture stress interference based on the DDM by further improving Olson’s height-correction formula; the modelling results suggested that the mechanism that is key to mitigating stress–shadow effects is to maintain uniform flow resistance within each fracture [21]. Tang et al. constructed a multi-fracture propagation model based on a planar three-dimensional DDM to present complex fracture geometries in a heterogeneous stress field, and the calculation results showed that longitudinal in situ stress distribution could affect the strength of the mutual interference between fractures [22]. Zheng et al. investigated the influence of fracturing sequences on fracture propagation; fracture geometries were analyzed under different fracturing sequences [23].
From the literature review, it can be concluded that the current research mainly focused on multi-fracture propagation and the interaction mechanism between natural fractures and hydraulic fractures; fracture propagation between wells is rarely discussed. In addition, the DDM, FEM, and DEM were the main methods of exploring fracture propagation in numerical simulations. In these simulations, the hydraulic fracture propagation trajectory was predetermined, and it cannot consider the effect of the stress field change on fracture propagation. Comparatively, the XFEM can simulate the effect of the stress field change on fracture propagation; thus, the XFEM was selected for the present study. In this work, we present a numerical investigation of fracture propagation for fractured parent wells and child wells in the Junggar Basin. The first challenge is to describe the inhomogeneous distribution of geostress; then, we must capture the fracture bending process with a consideration of the anisotropic characteristics of the reservoir. Here, we determine the rock parameters using logging data (Section 2.1 and Section 2.2) and estimate the magnitude and direction of 3D principal stresses through inversion analysis based on geostress field test data (Section 2.3 and Section 2.4). Then, we adopted a hydro-mechanical coupling extended finite element method (XFEM) to investigate the fracturing paths of the parent well and the child well, and the relevant parameters that affect the asymmetry of child well fractures are discussed (Section 3).

2. Methods

2.1. Rock Mechanical Tests

In this work, nineteen core samples were collected from the parent well in a sand–mudstone interlayered reservoir in the Junggar Basin for a rock mechanic test. Cylinders with a diameter of 25 mm and a length of about 50 mm were prepared. In the experiment, the loading speed of the samples was 0.005 mm/min, and the samples were cut from a drill core. For the uniaxial test, the cylinders were placed in the loading cell; then, the axial load increased at a constant rate. The stress–strain curves were recorded along this process. The test results on the 10 rock samples are shown in Table 1.

2.2. Determining Rock Parameters with Logging Data

In order to continuously estimate the rock mechanical parameters, the dynamic Young’s modulus (Ed) and dynamic Poisson’s ratio (νd) were calculated using acoustic logging data by Empirical Correlations (1) and (2):
E d = D E N D T S 2 ( 3 D T S 2 4 A C 2 D T S 2 4 A C 2 )
μ d = D T S 2 2 A C 2 2 ( D T S 2 4 A C 2 )
By fitting dynamic rock mechanical parameters to the static rock mechanical parameters, the conversion relations can be presented as follows:
E s = 1.1572 E d 6038.3
μ s = 1.8506 E μ d 0.2502
U C S = 45.3 e x p ( 0.0196 E s )
S T = 6169.32 ( t s / ρ b ) 1.422
In order to verify the usability of the transformation models in Figure 1, we compared the core test results with the rock mechanical curves determined by Equations (3) and (4). Based on these models and the logging data of the parent well, the curves of Young’s modulus and Poisson’s ratio at a depth of 643–645 m are presented in Figure 1; the core test results at the same depth (the black dots) present a good correlation with the curves. This is an indication of the accuracy of these transformation models. We applied these models to the entire well section of the parent well. Figure 2 shows the rock mechanic parameter curves in the depth range of 940–1030 m.

2.3. Estimation of the Pore Pressure

In view of the low mud content of the sand–mudstone alternate formation, we adopt an equivalent depth method to carry out the prediction of pore pressure. For rocks of different depths but with the same physical properties, the equivalent depth method assumes that the effective stress on the rock skeleton is equal [5]:
( σ α p p p ) 1 = ( σ α p p p ) 2
The primary premise of using the equivalent depth method is to predict formation pressure based on logging data; accordingly, one can establish an accurate mudstone normal compaction trend equation. Therefore, according to the logging data, the normal compaction trend equations of the sand and mudstone formations in the study area were established, respectively. Based on the established normal compaction trend line, if the abnormal pressure segment had the same acoustic time difference as a point on the normal compaction trend line, then the two points are indicated to have the same degree of compaction. The pore pressure of the target layer can be calculated as follows:
P A = ρ r A g h A g ( ρ r N ρ w N ) h N
where h A and h N are the depth of the target layer and the equivalent depth point, ρ r A and ρ r N are therock density, and ρ w N is the average pore fluid density.

2.4. Estimation of the Principal Stresses

The in situ geostress at the perforation was obtained by a hydraulic fracturing test. According to the principle of the lowest energy, a hydraulic fracture is always initiated and propagated along the direction that is perpendicular to the minimum principal stress ( σ m i n ); thus, the fracture indicates the direction of the maximum horizontal principal stress ( σ m a x ). The closing pressure of the fracture can be used as an indicator for the magnitude of the horizontal minimum principal stress; the horizontal maximum principal stress could be further calculated from the breaking pressure. Under fully elastic conditions, when drilling the borehole parallel to a certain principal stress direction, the circumferential stress ( σ θ ) and the normal stress ( σ r ) near the wellbore could be presented as follows:
σ θ = 1 2 σ m a x + σ m i n 1 + r a 2 r i 2 + 1 2 σ m a x σ m i n 1 + 3 r a 4 r i 4 c o s 2 θ P w r a 2 r i 2
σ r = 1 2 σ m a x + σ m i n 1 + r a 2 r i 2 + 1 2 σ m a x σ m i n 1 4 r a 4 r i 4 + 3 r a 4 r i 4 c o s θ + P w r a 2 r i 2
where r a is the wellbore radius, r i is the radial distance from the wellbore center, and P w is the fluid pressure in the wellbore. The shear stress could be presented as follows:
τ r θ = 1 2 σ m a x σ m i n 1 + 2 r a 2 r i 2 3 r a 4 r i 4 s i n 2 θ
On the wellbore wall, when → r a and θ = 0:
σ θ = 3 σ m i n σ m a x P w
σ m i n = P s o + ρ h P m P b
σ m a x = σ m i n + 4 ( P e P s o )
σ v = g i = 1 n ρ i D i
where P m is the friction loss along the well section, P b is the formation pore pressure, P e is the fracture extension pressure, and P s o is the instantaneous shut-in pressure. P m , P b , P e , and P s o are from the field fracturing data; the stress state at the perforation can be calculated by Equations (13)–(15). To estimate the principal stresses along the well continuously, here, we adopt the combined spring model. This model assumes that the rock is a homogeneous and isotropic linear elastic body without relative displacement during tectonic movements; the horizontal strain keeps constant. With these assumptions, the combined spring model could comprehensively consider the influence of formation rock mechanics, pore pressure, and tectonic action on principal stresses. Each principal stress component could be presented as follows:
σ m a x = ( μ 1 μ + β 1 ) ( σ v α P p ) + α P w
σ m i n = ( μ 1 μ + β 2 ) ( σ v α P p ) + α P w
The tectonic stress coefficients β 1 and β 2 can be determined by the stress state at the perforation; thus, the geostress along the wellbore can be obtained by Equations (16) and (17). The principal stresses along parent well L25 are shown in Figure 3.

2.5. Three-Dimensional Geostress Field

In situ stress forms under the combined influences of the weight of the rock mass and the tectonic movement, as well as the physicochemical changes in the geological mass, topography, and temperature stress. The geostress field test that was conducted during the operation was the most direct approach for studying deep geostress; however, the cost of in situ stress testing was high, the shapes of the deep geological structures were complex, the distribution of the formation media was uneven and discontinuous, and the rock physical and mechanical properties were variable. The test results only presented the stress conditions in the area near the test point. Therefore, a feasible method was to construct a reasonable geomechanical model on the basis of a fine analysis of the geological structures, and to perform a numerical simulation inversion analysis based on the in situ stress data of limited test points and mathematical and mechanical theories; this approach enabled us to further quantitatively study the stress field in a deep area and analyze the horizontal distribution of the stress field in the deep area.
Taking the tested in situ stress as a basic constraint, the inversion analysis of the 3D geostress could be accomplished by determining the boundary load of the calculation model. The far-field tectonic stress state could be considered as the superposition of the following two basic tectonic states: (a) horizontal extrusion or stretching along the X and Y directions and (b) uniform shear structure deformation. Thus, the displacement boundaries in the model can be expressed as follows:
u x = P x n + P t 1 u y = P y n + P t 2
where u x and u y are the horizontal loading vectors in the X and Y directions; P x n and P y n are the horizontal extrusion or tensile vectors; P t 1 and P t 2 are the horizontal shearing vectors.
The initial stress field inversion was based on an iterative algorithm that gradually corrects the trial calculation of parameters until the error function tended to a minimum. The geostress component can be noted as follows:
σ i j = f P x n , P y n , P t 1 P t 2
The field geostress component was noted as σ i j m k ; at the same spot, the geostress component calculated by the finite element model was σ i j c k . The error function, ψ , can be expressed using the deviation of the calculated stress from the measured stress:
ψ P x n , P y n , P t 1 , P t 2 = σ i j c k σ i j m k / σ i j c k
In this work, numerical inversion was accomplished by a genetic algorithm in MatLab, 2020. The boundary vectors P x n , P y n , P t 1 , a n d P t 2 are noted as input parameters, and the geostress components, σ i j , calculated by the finite element method, are noted as output parameters. The initial error decreased gradually with the genetic generation (Figure 4a); after the inheritance was conducted to the 500th generation, the inversion accuracy rate reached 93.3%. The inversed 3D principal stress fields of an 8.0 × 5.0 km2 area at a depth of 100–1200 m are presented in Figure 4b–d.
The rock mechanical properties and stress fields, calculated according to the experimental tests and logging calculations, are depicted in Table 2.

3. Hydraulic Fracture Modelling

3.1. Modelling Approach

The extended finite element method (XFEM) uses the level set method or the fast advance method to introduce the modification of the shape function within the unit. The continuous function was used to characterize the continuous region, and the enhanced function was used to represent the non-continuous region. The introduced enhancement functions mainly include three functions: the Heaviside function to characterize the discontinuity of the displacement field, the Junction function to characterize the fracture crossing, and the fracture tip increase function to characterize the fracture tip singularity field. The corresponding unit types are the penetration node, the cross node, and the split-tip node (Figure 5).
The approximate form of displacement u r at any point in the finite element model with cracks is as follows:
u r = i M r N α i ( r ) u i + j M r N α j r H r H r j a j                                                                                                                         + h M r N α h r J r J r h b h + k M r N α k r i = 1 4 F l r F l r k c k l
where r is the position vector; M r , M r j , M r h , and M r k are the ordinary node set, step function enhanced node set, connection function node set, and fracture tip enhanced node set of the unit where r is located; u i is the node displacement; N α α = 1 4 is a shape function; α ( i ) is the mapping function from global coordinates to local coordinates;   H ( r ) , J ( r ) , and F l ( r ) are the Heaviside function value, the Junction function value, and the fracture tip enhancement function value of Gaussian point x, respectively. a j , b h , and c k l are the enhanced node degrees of freedom vectors for the penetrating element, the intersection element, and the fracture tip element, respectively. The Heaviside function is a step function using level set function to determine whether a point is located on the upper or lower surface of the crack, thereby indicating the discontinuity of the displacement field:
H r = s g n ( f ( r ) )
where the level set function f ( r ) = 0 represents the fracture surface. If point r was closest to point r* on the fracture surface, the normal vector at point r* was e n , and when e n · r r * > 0 and f r > 0 , point r was located above the fracture surface. The Junction function was a connection function that reflects the effect of the intersection of cracks:
J r = s g n f 2 r s g n f 2 r h   f 2 r · f 2 r 0 > 0 s g n f 1 r s g n f 1 r h   f 2 r · f 2 r 0 < 0
where sgn (r) is the sign function; f 1 r is the level set function of hydraulic fractures; f 2 r is the level set function of natural fractures. The fracture enhancement function is to extract the function reflecting the singularity of the fracture tip according to the progressive displacement field of the fracture tip in linear elastic fracture mechanics, namely:
F r , φ l = 1 4 = r s i n φ / 2 , c o s φ / 2 , s i n φ / 2 s i n φ , c o s φ / 2 s i n φ
where (r, φ ) represents the polar coordinates of the local coordinate system centered on the fracture tip.

3.2. Hydraulic Fracture Modelling at the Parent Well

The targeted reservoir was at a depth of 620 m; the material parameters are calculated by Equations (3)–(6) with logging data. Based on Saint-Venant’s principle, the 3D geostress field near the wellbore was abstracted from the inversed principal stress field in Figure 6a and shown in Figure 6c. The hydraulic fractures extend in a perpendicular direction from the maximum principal stress σ v , and form a vertical fracture in the target reservoir. In the horizontal direction, the range of the model was 1000 m × 1000 m (Figure 6d). According to the fracturing design at parent well L25 provided by the industry, the fracture operation included four stages, and the time was close to the real time.
  • Stage 1: The first fracturing occurred at 50 min; the flow ratio of the fracturing fluid was 2 m3/min.
  • Stage 2: Fracturing fluid flowback occurred for 360 min.
  • Stage 3: Fracturing repeated for 50 min; the flow ratio of the fracturing fluid was 2 m3/min.
After the fracturing operation at the parent well, the industry planned to place a child well at candidate location (Figure 6c). The fracturing operation at the child well occurred at 100 min; the flow ratio of the fracturing fluid was 2 m3/min.

3.2.1. First Fracturing Modelling

The evolution of fluid pressure at the perforation with time is presented in Figure 7. The extrusion stage was not included in this work; at the beginning of the fracturing process, the fluid pressure at the perforation point rises sharply until a breaking pressure is reached. Then, the fluid pressure drops rapidly; following this, it basically maintains a stable value during the subsequent fracture propagation process. The accuracy of the numerical model was verified indirectly using field fracturing data of the parent well. The field fracture pressure was 19.30 MPa and the closing pressure was 14.73 MPa. In our modeling results, the simulated fracture pressure was 18.13 MPa and the closure pressure was 12.33 MPa; the agreement between the simulation results and the field fracturing data was 87.32%, indicating the feasibility of the proposed approach.
At the beginning of the fracturing process, a fracture propagates in the direction of the perforation. After the fracture length reaches about 10 times the length of the borehole’s diameter, it propagates along the direction of the maximum principal stress, and thus forms a turning fracture. The induced stress was tensile at the fracture tip; in the normal direction of the fracture length, it was compressive. This induced compressive stress changes the stress state near the fracture; the direction of σ m a x changed from NNW–SSE to E–W (Figure 7b). The direction of the minimum horizontal principal stress before and after the initial fracturing changes accordingly (Figure 7d). However, the direction of the horizontal principal stress near the fracture tip was not significantly affected; for this reason, the propagation direction of the fracture was not affected by the induced stress.

3.2.2. Fracturing Fluid Flowback and Repeat Fracturing

After the fracturing fluid flowback for 360 min, the pore pressure near the fracture dropped by about 3~5 MPa; the magnitudes of σ m a x and σ m i n changed accordingly, but their directions were not affected. The fracture morphology did not change significantly during the flowback process (Figure 8a,c). Then, during the repeated fracturing process at well L25, the fracture width gradually increased in the first 20 min, but the fracture length did not change significantly. After 30 min, the fracture continued to propagate along the direction of the maximum principal stress; this range is influenced by fracturing-induced stress, which further expands here (Figure 8b,d).

3.3. Modelling Result Discussion

After the first fracturing process, due to the induced tensile stress at the fracture tips and the induced compressive stress in the normal direction of the fracture length, the area affected by fracturing-induced stress forma a “butterfly type” area (shown in Figure 9a). In the E–W direction, the range affected after the first fracturing was 1.5–2 times the length of the fracture. After repeated fracturing, the magnitude of the induced stress was increased, but its influence range did not change significantly (Figure 9b). The fluid pressure variation at the parent well was present in Figure 9c; as shown here, the fluid pressure rose rapidly during the first fracturing process, until it reached the breaking pressure; then, it decreased throughout the fracturing fluid flowback process. Afterward, the repeat fracturing enhanced the fluid pressure again. Here, we investigate the variation in principal stresses at the checkpoint during parent well fracturing. The checkpoint is located 200 m east from the parent well; this is inside the “butterfly type” affected area that is shown in Figure 9a. As shown in Figure 9d, the maximum principal stress turned into the minimum principal stress after the first fracturing process; the repeated fracturing led to an increase in the difference between the maximum and minimum principal stresses; the reorientation of the stresses has a significant impact on the fracturing of the child well here. As indicated in Figure 9a, the checkpoint was the point at which the geostress and the pore pressure were measured during the fracturing of the parent well.

3.4. Fracturing of the Child Well

In the candidate child well location, the induced stress field was superimposed onto the in situ stress field; this stress state change can significantly impact the fracturing of the child well. Here, we present several child wells in different locations to illustrate how inter-well interference affects the fracture wings of the child well. After the parent well repeatedly fractured, a child well fracturing was conducted for d = 200 m, d = 250 m, and d = 300 m, where d represents the distance from the child well to the parent well on the E–W direction. Figure 10 presents the fracture geometry after fracturing at the child well for 100 min, with a flow ratio of 2 m3/min; the fracture width (white curve) was magnified 50 times for visualization. Child well C1 is located within the “butterfly type” zone in Figure 9a, while child wells C2 and C3 are located outside this zone. As shown in Figure 10b, the left fracture wing of the child well propagates in the east–west direction, which is significantly affected by the compressive stress induced by parent well fracturing; meanwhile, the right wing propagates in the northwest–southeast direction and thus present asymmetric fracture propagation. When the child well is located outside the “butterfly type” zone, Figure 10c,d present the symmetric fracture geometry. The effect of stress interference on the asymmetry of child well fracture wings was negatively correlated with d.

4. Conclusions

In this study, we investigated stress interference among hydraulic fractures with nonuniform stress distribution in an inhomogeneous reservoir. We determined the rock parameters with rock mechanical tests and field logging data, and then estimated the magnitude and direction of 3D principal stresses by inversion analysis; this analysis was based on geostress field test data. The hydromechanical coupling XFEM was adopted to investigate the fracturing paths for the fracture wings of the parent well and the child well, and the hydraulic fracturing at different child well locations were also studied. The conclusions can be summarized as follows:
(1)
Due to the induced tensile stress at the fracture tips and the induced compressive stress in the normal direction of the fracture length, the area affected by fracturing-induced stress formed a “butterfly type” area. In the E–W direction, the affected range after the first fracturing process at the parent well was 1.5–2 times the fracture length.
(2)
The fracturing fluid flowback at the parent well changed the magnitude of the maximum and minimum horizontal principal stresses, but their directions were not affected. After repeated fracturing at the parent well, the magnitude of the fracturing-induced stress increased, but its influence range did not increase significantly.
(3)
For a child well located within the “butterfly type” zone, the stress interference results in an asymmetric fracture propagation; meanwhile, for a child well located outside of this zone, a symmetric fracture geometry occurs.
(4)
The non-concurrent parent well and child well completions are the main cause of stress interference, and the effect of stress interference on the asymmetry of the child well fracture wings was negatively correlated with the distance between the parent well and the child well.
(5)
The fracture geometry under different well spacings was analyzed; this provided a method of optimizing the well space during parent well fracturing. In the simulations, the formation properties were assumed to be uniform. This was necessary for investigating the impact of nonuniform stress distribution on fracture propagation.

Author Contributions

Conceptualization, Y.Z., P.W., R.L. and H.Z.; Methodology, Y.Z. and Y.L.; Software, Y.Z., P.W., Y.L., R.L. and H.Z.; Validation, Y.Z. and H.Z.; Formal analysis, R.L.; Investigation, Y.Z., P.W., Y.L. and H.Z.; Resources, R.L. and H.Z.; Writing—original draft, Y.Z., Y.L., R.L. and H.Z. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported by National Natural Science Foundation of China. Grant Number: [62173049].

Data Availability Statement

The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest

Authors Yinghui Zhu, Yi Liao, Ruiquan Liao and Heng Zheng were employed by the company Research Institute of CNOOC Shenzhen Branch. Author Pengxiang Wang was employed by the company Downhole Service Company, CNPC Chuanqing Drilling Engineering Company Limited.

References

  1. Bowers, G.L. Pore pressure estimation from velocity data: Accounting for overpressure mechanisms besides undercompaction. SPE Drill. Complet. 1995, 10, 89–95. [Google Scholar] [CrossRef]
  2. Cheng, Y. Impacts of the number of perforation clusters and cluster spacing on production performance of horizontal shale-gas wells. SPE Reserv. Eval. Eng. 2012, 15, 31–40. [Google Scholar] [CrossRef]
  3. Esquivel, R.; Blasingame, T.A. Optimizing the Development of the Haynesville Shale-Lessons-Learned from Well-to-Well Hydraulic Fracture Interference. In Proceedings of the Unconventional Resources Technology Conference, Austin, TX, USA, 24–26 July 2017; Society of Exploration Geophysicists, American Association of Petroleum: Houston, TX, USA, 2017; pp. 1005–1026. [Google Scholar]
  4. Gala, D.P.; Manchanda, R.; Sharma, M.M. Modeling of fluid injection in depleted parent wells to minimize damage due to frac-hits. In Proceedings of the Unconventional Resources Technology Conference, Houston, TX, USA, 23–25 July 2018; Society of Exploration Geophysicists, American Association of Petroleum: Houston, TX, USA, 2018; pp. 4118–4131. [Google Scholar]
  5. Guo, X.; Ma, J.; Wang, S.; Zhu, T.; Jin, Y. Modeling Interwell Interference: A Study of the Effects of Parent Well Depletion on Asymmetric Fracture Propagation in Child Wells. In Proceedings of the SPE/IATMI Asia Pacific Oil & Gas Conference and Exhibition, Bali, Indonesia, 29–31 October 2019. [Google Scholar]
  6. Hu, X.; Liu, G.; Luo, G.; Ehlig-Economides, C. Model for Asymmetric Hydraulic Fractures with Nonuniform-Stress Distribution. SPE Prod. Oper. 2019, 36, 1–11. [Google Scholar] [CrossRef]
  7. Lawal, H.; Jackson, G.; Abolo, N.; Flores, C. A novel approach to modeling and forecasting frac hits in shale gas wells. In Proceedings of the EAGE Annual Conference & Exhibition Incorporating SPE Europec, London, UK, 10–13 June 2013; Society of Petroleum Engineers: Richardson, TX, USA, 2013. [Google Scholar]
  8. Lindsay, G.; Miller, G.; Xu, T.; Shan, D.; Baihly, J. Production performance of infill horizontal wells vs. pre-existing wells in the major US unconventional basins. In Proceedings of the SPE Hydraulic Fracturing Technology Conference and Exhibition, The Woodlands, TX, USA, 23–25 January 2018; Society of Petroleum Engineers: Richardson, TX, USA, 2018. [Google Scholar]
  9. Miller, G.; Lindsay, G.; Baihly, J.; Xu, T. Parent well refracturing: Economic safety nets in an uneconomic market. In Proceedings of the SPE Low Perm Symposium, Denver, CO, USA, 5–6 May 2016; Society of Petroleum Engineers: Richardson, TX, USA, 2016. [Google Scholar]
  10. Moës, N.; Belytschko, T. Extended finite element method for cohesive crack growth. Eng. Fract. Mech. 2002, 69, 813–833. [Google Scholar] [CrossRef]
  11. Montgomery, C.T.; Smith, M.B. Hydraulic fracturing: History of an enduring technology. J. Pet. Technol. 2010, 62, 26–40. [Google Scholar] [CrossRef]
  12. Olson, J.E. Predicting fracture swarms-The influence of subcritical crack growth and the crack-tip process zone on joint spacing in rock. Geol. Soc. Lond. Spec. Publ. 2004, 231, 73–88. [Google Scholar] [CrossRef]
  13. Rezaei, A.; Rafiee, M.; Bornia, G.; Soliman, M.; Morse, S. Protection Refrac: Analysis of pore pressure and stress change due to refracturing of Legacy Wells. In Proceedings of the Unconventional Resources Technology Conference, Austin, TX, USA, 24–26 July 2017; Society of Exploration Geophysicists, American Association of Petroleum: Houston, TX, USA, 2017; pp. 312–327. [Google Scholar]
  14. Salmachi, A.; Sayyafzadeh, M.; Haghighi, M. Infill well placement optimization in coal bed methane reservoirs using genetic algorithm. Fuel 2013, 111, 248–258. [Google Scholar] [CrossRef]
  15. Soliman, M.Y.; East, L.; Adams, D. Geomechanics aspects of multiple fracturing of horizontal and vertical wells. In Proceedings of the SPE International Thermal Operations and Heavy Oil Symposium and Western Regional Meeting, Bakersfield, CA, USA, 16–18 March 2004; Society of Petroleum Engineers: Richardson, TX, USA, 2004. [Google Scholar]
  16. Stephenson, B.; Galan, E.; Williams, W.; MacDonald, J.; Azad, A.; Carduner, R.; Zimmer, U. Geometry and failure mechanisms from microseismic in the Duvernay Shale to explain changes in well performance with drilling azimuth. In Proceedings of the SPE Hydraulic Fracturing Technology Conference and Exhibition, The Woodlands, TX, USA, 23–25 January 2018; Society of Petroleum Engineers: Richardson, TX, USA, 2018. [Google Scholar]
  17. Surendran, M.; Natarajan, S.; Palani, G.; Bordas, S.P. Linear smoothed extended finite element method for fatigue crack growth simulations. Eng. Fract. Mech. 2019, 206, 551–564. [Google Scholar] [CrossRef]
  18. Tang, H.; Wang, S.; Zhang, R.; Li, S.; Zhang, L.; Wu, Y.-S. Analysis of stress interference among multiple hydraulic fractures using a fully three-dimensional displacement discontinuity method. J. Pet. Sci. Eng. 2019, 179, 378–393. [Google Scholar] [CrossRef]
  19. Tang, H.; Winterfeld, P.H.; Wu, Y.-S.; Huang, Z.-Q.; Di, Y.; Pan, Z.; Zhang, J. Integrated simulation of multi-stage hydraulic fracturing in unconventional reservoirs. J. Nat. Gas Sci. 2016, 36, 875–892. [Google Scholar] [CrossRef]
  20. Walser, D.; Siddiqui, S. Quantifying and mitigating the impact of asymmetric induced hydraulic fracturing from horizontal development wellbores. In Proceedings of the SPE Annual Technical Conference and Exhibition, Dubai, United Arab Emirates, 26–28 September 2016; Society of Petroleum Engineers: Richardson, TX, USA, 2016. [Google Scholar]
  21. Wu, K.; Olson, J.; Balhoff, M.T.; Yu, W. Numerical analysis for promoting uniform development of simultaneous multiple-fracture propagation in horizontal wells. SPE Prod. Oper. 2017, 32, 41–50. [Google Scholar] [CrossRef]
  22. Yaich, E.; Diaz De Souza, O.C.; Foster, R.A.; Abou-Sayed, I.S. A Methodology to Quantify the Impact of Well Interference and Optimize Well Spacing in the Marcellus Shale. In Proceedings of the Spe/csur Unconventional Resources Conference—Canada, Calgary, AB, Canada, 30 September–2 October 2014. [Google Scholar]
  23. Zheng, H.; Pu, C.; Xu, E.; Sun, C. Numerical investigation on the effect of well interference on hydraulic fracture propagation in shale formation. Eng. Fract. Mech. 2020, 228, 106932. [Google Scholar] [CrossRef]
Figure 1. Validation of the applicability of the established transformation models; the black dots represent the core test results.
Figure 1. Validation of the applicability of the established transformation models; the black dots represent the core test results.
Processes 12 02089 g001
Figure 2. Rock mechanic parameters in the depth range of 940–1030 m.
Figure 2. Rock mechanic parameters in the depth range of 940–1030 m.
Processes 12 02089 g002
Figure 3. Principal stresses along parent well L25.
Figure 3. Principal stresses along parent well L25.
Processes 12 02089 g003
Figure 4. Error evolution of the genetic algorithm and the inversed principal stress field.
Figure 4. Error evolution of the genetic algorithm and the inversed principal stress field.
Processes 12 02089 g004
Figure 5. Unit types of the extended finite element method.
Figure 5. Unit types of the extended finite element method.
Processes 12 02089 g005
Figure 6. Horizontal stress field near wellbore; the candidate child well location was 200 m east of parent well L25.
Figure 6. Horizontal stress field near wellbore; the candidate child well location was 200 m east of parent well L25.
Processes 12 02089 g006
Figure 7. Horizontal stress field before and after first fracture.
Figure 7. Horizontal stress field before and after first fracture.
Processes 12 02089 g007
Figure 8. Horizontal stress field before and after repeat fracturing.
Figure 8. Horizontal stress field before and after repeat fracturing.
Processes 12 02089 g008aProcesses 12 02089 g008b
Figure 9. (a) The “butterfly type” area affected by fracturing-induced stress; the checkpoint is located 200 m east of the parent well. (b) The affected range after the first fracturing process and repeat fracturing in the E–W direction. (c) Fluid pressure variation at the parent well. (d) Principal stresses variation at the checkpoint during the fracturing of the parent well.
Figure 9. (a) The “butterfly type” area affected by fracturing-induced stress; the checkpoint is located 200 m east of the parent well. (b) The affected range after the first fracturing process and repeat fracturing in the E–W direction. (c) Fluid pressure variation at the parent well. (d) Principal stresses variation at the checkpoint during the fracturing of the parent well.
Processes 12 02089 g009
Figure 10. The distribution of maximum principal stress and fracture geometry after fracturing at the child well for 100 min with a flow ratio of 2 m3/min for different parent well–child well distances.
Figure 10. The distribution of maximum principal stress and fracture geometry after fracturing at the child well for 100 min with a flow ratio of 2 m3/min for different parent well–child well distances.
Processes 12 02089 g010
Table 1. Test results of the 10 rock samples.
Table 1. Test results of the 10 rock samples.
Sample IDDepth (m)Density (g/cm3)Porosity (%)Confining Pressure (MPa)Compressive Strength (MPa)Poisson’s Ratio (1)Es (MPa)
1-1645.352.591.50.062.5440.30654,760
1-2645.382.601.10.053.2020.05710,544
1-3645.672.631.23.092.2000.26647,382
1-4645.722.611.86.0118.3850.21840,955
1-5645.872.641.43.089.4520.21537,995
1-6645.832.641.89.0107.8140.24137,939
1-7646.022.631.16.0108.8870.21738,042
1-10646.542.621.39.081.3000.19528,079
Table 2. Parameters of the hydraulic fracturing simulation.
Table 2. Parameters of the hydraulic fracturing simulation.
ItemParametersUnitValue
Rock Mechanical propertiesYoung’s ModulusGPa27.6
Poisson’s RateDimensionless0.34
Tensile strengthMPa8.42
Shear strengthMPa8.61
PermeabilitymD0.35
Injection ParametersInjection ratem3/min2.0
Viscosity of fracturing fluidmPa.s20.0
Stress FieldVertical stressMPa15.0
Maximum horizontal principal stressMPa13.6
Minimum horizontal principal stressMPa11.7
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Zhu, Y.; Wang, P.; Liao, Y.; Liao, R.; Zheng, H. Use of Extended Finite Element Method to Characterize Stress Interference Caused by Nonuniform Stress Distribution during Hydraulic Fracturing. Processes 2024, 12, 2089. https://doi.org/10.3390/pr12102089

AMA Style

Zhu Y, Wang P, Liao Y, Liao R, Zheng H. Use of Extended Finite Element Method to Characterize Stress Interference Caused by Nonuniform Stress Distribution during Hydraulic Fracturing. Processes. 2024; 12(10):2089. https://doi.org/10.3390/pr12102089

Chicago/Turabian Style

Zhu, Yinghui, Pengxiang Wang, Yi Liao, Ruiquan Liao, and Heng Zheng. 2024. "Use of Extended Finite Element Method to Characterize Stress Interference Caused by Nonuniform Stress Distribution during Hydraulic Fracturing" Processes 12, no. 10: 2089. https://doi.org/10.3390/pr12102089

APA Style

Zhu, Y., Wang, P., Liao, Y., Liao, R., & Zheng, H. (2024). Use of Extended Finite Element Method to Characterize Stress Interference Caused by Nonuniform Stress Distribution during Hydraulic Fracturing. Processes, 12(10), 2089. https://doi.org/10.3390/pr12102089

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop