# Deformation Behavior between Hydraulic and Natural Fractures Using Fully Coupled Hydromechanical Model with XFEM

Authors: Fei Liu; Zhifeng Luo; Yu Sang; Liqiang Zhao; Changlin Zhou
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

There has been a growing consensus that preexisting natural fractures play an important role during stimulation. A novel fully coupled hydromechanical model using extended finite element method is proposed. This directly coupled scheme avoids the cumbersome process during calculating the fluid pressure in complicated fracture networks and translating into an equivalent nodal force. Numerical examples are presented to simulate the hydraulic fracture propagation paths for simultaneous multifracture treatments with properly using the stress shadow effects for horizontal wells and to reveal the deformation response and interaction mechanism between hydraulic induced fracture and nonintersected natural fractures at orthotropic and nonorthotropic angles. With the stress shadow effects, the induced hydraulic flexural fracture deflecting to wellbore rather than transverse fracture would be formed during the progress of simultaneous fracturing for a horizontal well. The coupled hydromechanical simulation reveals that the adjacent section to the intersection is opened and the others are closed for orthogonal natural fracture, while the nonorthogonal natural fracture is activated near the intersection firstly and along the whole section with increasing perturbed stresses. The results imply that the induced hydraulic fracture tends to cross orthotropic natural fracture, while it is prior to being arrested by the nonorthotropic natural fracture.


---

## Body

## 1. Introduction

Unconventional reservoirs, such as shale gas, have an extremely tight matrix, but well-developed natural fractures [1]. Volume stimulation drives these natural fractures to expand continuously and introduce shear slip in brittle rocks, thus forming complicated networks of intersecting natural fractures and artificial fractures that further increase the stimulated volume [2]. Therefore, the volume stimulation method has become one of the key technologies in unconventional oil and gas development [3]. However, simulation hydraulic fracturing in the presence of natural fracture network is a challenging task, owing to the complex interactions between fracturing fluid, rock matrix, and rock interfaces, as well as the mechanical behaviors before, during, and after the intersections between hydraulic fractures and preexisting natural fractures [4–6].The morphology of complicated fracture networks created by volume fracturing can be estimated by determining the mechanical mechanism and propagation behavior of hydraulically induced fractures before, during, and after intersecting with natural fractures. In laboratory experiments, a true triaxial fracturing simulation system can be used in combination with high-energy CT scanning, acoustic emission, and microscopic X-ray chromatography to characterize the intersection behavior between hydraulically induced fractures and natural fractures during the fracturing of artificial rock or reservoir outcrop rock samples. These methods also reveal the influences of in situ stress state, horizontal principle stress difference, approach angle, natural fracture interface properties, injection fluid amount, and fracturing fluid viscosity [7–11]. Scholars have developed a series of criteria to describe the intersection behavior between hydraulic fractures and natural fractures. For example, Renshaw and Pollard proposed a criterion for the intersection between hydraulic fractures and orthogonal fractures (R-P criterion) [12]. Gu and Weng expanded the R-P criterion to encompass nonorthogonal conditions (G-W criterion) by considering that the approach angles of hydraulic fractures and natural fractures are between 0° and 90° [13]. Meanwhile, Chen et al. considered the dip angle of natural fractures and expanded the R-P criterion to three-dimensional conditions [14]. Dahi-Taleghani and Olson investigated the tensile property and shearing behavior of natural fractures while hydraulic fractures are approaching but before they come into contact [15]. All of these criteria only consider the mechanical behavior caused by the interaction between hydraulic fractures and natural fractures, without accounting for the influence of the hydraulic effect induced by the injection of fracturing fluid.Volume fracturing presents a complicated hydromechanical coupling problem, with two special conditions. Firstly, there would be fluid flow within the fracture network and matrix rock and fluid exchange or fluid leak between the fracture and the surrounding porous media, because of the pressure difference between the fracture network and the matrix rock during fracturing fluid pumping. Secondly, although the net pressure drives the initiation and propagation of fractures, the increased pore fluid pressure in the matrix surrounding the fractures would impact fracture propagation. There were considerable differences between piston-like displacement, single phase and multiphase leakage behaviors for nonintersected fractures, and matrix rock using extended finite element method (XFEM) [16]. Based on the vital important role of leakage, this paper promoted the multiphase fluid flow model to fully coupled hydromechanical model for fractured porous media to provide a numerical tool for simulating hydraulic fracturing of fractured porous media. This decision is enlightened on the strong discontinuity of displacement at either side of fractures and the weak discontinuity of pore fluid pressure, as XFEM represents the most commonly used method in simulating discontinuous media. This paper also discusses the deformation features and mechanical mechanisms before the intersection between hydraulically induced fractures and orthogonal and nonorthogonal natural fractures.


## 2. Fully Coupled Hydromechanical Model

The fractured porous media were composed of matrix rock and macro fractures as a whole. Although there were plenty of micro fractures, which show performances similar to those of macro fractures, developed in matrix rock, the fractured porous media were considered to be constructed by an intact matrix and macro fracture networks for simplicity, as displayed in Figure1.Figure 1
Schematic diagram of hydromechanical coupling for fractured porous media.The strong form, weak form, and the associated discretization form of governing equations for hydromechanical coupling inside both matrix domain and fracture domain are demonstrated in this section.


### 2.1. Strong Form

Figure2 shows a bounded control body Ω⊂Rn (n=2,3) with an outer boundary of Γ.Figure 2
Control body and boundary conditions of fractured porous media.The effective stress equilibrium equation for porous media is written as(1)∇·σ+ρg=∇·σ′-αIp+ρg=0.According to the law of mass conservation, the continuity equations for the wetting and nonwetting phase fluid flow in the matrix domain are expressed as(2)α-ϕKsSw,mSw,m+pc,m∂Sw,m∂pc,m+ϕSw,mKw-ϕ∂Sw,m∂pc,m∂pw,m∂t+α-ϕKsSw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pnw,m∂t+αSw,mIT∂ε∂t+∇·-kmkm,rwμw∇pw,m+QwstcBw=0,α-ϕKsSnw,mSw,m+pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pw,m∂t+α-ϕKsSnw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕSnw,mKnw-ϕ∂Sw,m∂pc,m∂pnw,m∂t+αSnw,mIT∂ε∂t+∇·-kmkm,rnwμnw∇pnw,m=0.Similarly, the continuity equations for the wetting and nonwetting phase fluid flow inside the fracture domain are given by(3)Sw,fKw-∂Sw,f∂pc,f∂pw,f∂t+∂Sw,f∂pc,f∂pnw,f∂t+∇·-kfkf,rwμw∇pw,f+Sw,f∇·u˙=0,∂Sw,f∂pc,f∂pw,f∂t+Snw,fKnw-∂Sw,f∂pc,f∂pnw,f∂t+∇·-kfkf,rnwμnw∇pnw,f+Snw,f∇·u˙=0.The displacement and pressure of the wetting and nonwetting phase initially are provided according to the initial conditions(4)u=u0,pw,m=pw,m0,pnw,m=pnw,m0.The displacement and pressure of the wetting and nonwetting phase at the boundary are provided by the forced boundary conditions(5)ux=u~∀x∈Γuu,pw,mx=p~w,m∀x∈Γpw,pnw,mx=p~nw,m∀x∈Γpnw.The stress and mass flow rate of the wetting and nonwetting phase at the boundary are provided by the natural boundary conditions(6)σx·n=t~∀x∈Γut,-kmkm,rwμw∇pw,m·n=q~wρw∀x∈Γqw,-kmkm,rnwμnw∇pnw,m·n=q~nwρnw∀x∈Γqnw.The fluid pressure is applied on the fracture surface(7)σx·nΓf+=Sw,fpw,f+Snw,fpnw,f·nΓf∀x∈Γf+,σx·nΓf-=-Sw,fpw,f+Snw,fpnw,f·nΓf∀x∈Γf-.The pore fluid pressure at the fracture boundary is continuous, with fluid exchange occurring between fractures and the surrounding matrix:(8)pw,m=pw,f∀x∈Γf,pnw,m=pnw,f∀x∈Γf,-kmkm,rlμl∇pl,m·nΓf+=q~l+ρll=w,nw;∀x∈Γf+,-kmkm,rlμl∇pl,m·nΓf-=q~l-ρll=w,nw;∀x∈Γf-,-kfkf,rlμl∇pl,f·nΓf+=-q~l+ρll=w,nw;∀x∈Γf+,-kfkf,rlμl∇pl,f·nΓf-=-q~l-ρll=w,nw;∀x∈Γf-.


### 2.2. Weak Form

The weak form of the governing equations can be obtained by combining the weighted residual method with the associated boundary conditions. With regard to fluid flow in fractures, since the fracture width is much smaller than its length, the fluid pressure variation along its width can be neglected (see Figure3). The areal integral inside the fracture domain Ωf can be transformed into the linear integral along the fracture surface Γf. The fluid exchange items or leak between the fractures and the surrounding matrix rock can then be vanished by adding the equivalent integrals of the matrix rock and fractures, which have the same absolute value but opposite signs. Thus, it becomes unnecessary to explicitly express the fluid exchange items or leak [16–18]:(9)∫ΩδεTσ′dΩ-∫ΩδεTαΙSw,mpw,mdΩ-∫ΩδεTαΙSnw,mpnw,mdΩ-∫ΓfδuTRTIfnSw,fpw,fdΓ-∫ΓfδuTRTIfnSnw,fpnw,fdΓ=∫ΓutδuTt~dΓ,(10)∫ΩwpTα-ϕKsSw,mSw,m+pc,m∂Sw,m∂pc,m+ϕSw,mKw-ϕ∂Sw,m∂pc,m∂pw,m∂tdΩ+∫ΩwpTα-ϕKsSw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pnw,m∂tdΩ+∫ΩwpTαSw,mIT∂ε∂tdΩ+∫Ω∇wpTkmkm,rwμw∇pw,mdΩ+∫ΓqwwpTq~wρwdΓ+∫ΓfbwpTSw,fKw-∂Sw,f∂pc,f∂pw,f∂tdΓ+∫ΓfbwpT∂Sw,f∂pc,f∂pnw,f∂tdΓ+∫ΓfwpTSw,fu˙y′dΓ+∫Γfbkfkf,rwμw∂wpT∂x′∂pw,f∂x′dΓ=0,(11)∫ΩwpTα-ϕKsSnw,mSw,m+pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pw,m∂tdΩ+∫ΩwpTα-ϕKsSnw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕSnw,mKnw-ϕ∂Sw,m∂pc,m∂pnw,m∂tdΩ+∫ΩwpTαSnw,mIT∂ε∂tdΩ+∫Ω∇wpTkmkm,rnwμnw∇pnw,mdΩ+∫ΓqnwwpTq~nwρnwdΓ+∫ΓfbwpT∂Sw,f∂pc,f∂pw,f∂tdΓ+∫ΓfbwpTSnw,fKnw-∂Sw,f∂pc,f∂pnw,f∂tdΓ+∫ΓfwpTSnw,fu˙y′dΓ+∫Γfbkfkf,rnwμnw∂wpT∂x′∂pnw,f∂x′dΓ=0.Figure 3
Schematic diagram of fluid flows in a fracture.


### 2.3. Discretization Form

XFEM introduces Heaviside step functions and branch functions to enrich the shape function, which make it easy to deal with typical strong discontinuous problems like fractures without grid matching and remeshing. The XFEM displacement field for multiple fractures is [19](12)uhx=∑i∈SNixuai+∑j=1nc∑i∈SHjubi,jNixHjx-Hjxi︸enrichment  of  bisected  support  for  fracturej+∑j=1nt∑i∈SCj∗NixRjx∑l=14uci,jlψjlx-ψjlxi︸enrichment  of  slit  support  for  fracturej+∑j=1nx∑i∈SJjudi,jNixJjx︸enrichment  of  intersecting  support  for  cross  pointj.In fractured porous media, the natural fractures are the main channels for seepage. The pore fluid pressure field is continuous at the fracture surface, while the derivative (i.e., flow rate) is discontinuous. This presents a typical weak discontinuous problem. Similar to the interpolation of displacement field, the pressure field can be interpolated as [16, 18](13)phx=∑i∈SNixpai+∑j=1nc∑i∈SHjpbi,jNixfjx-fjxi︸enrichment  of  bisected  support  for  fracturej+∑j=1nt∑i∈SCj∗NixRjx∑l=13pci,jlBjlx-Bjlxi︸enrichment  of  slit  support  for  fracturej+∑j=1nx∑i∈SJjpdi,jNixJjpx-Jjpxi︸enrichment  of  intersecting  support  for  cross  pointj.The discretization form of the fully coupled hydromechanical model for fractured porous media can be obtained by regarding the trial functionsδu and wp in (9)–(11) as the variation of displacement and pressure interpolation functions, as shown in(14)0000Hw,m+Hw,f000Hnw,m+Hnw,fu¯p¯wp¯nw+-KLw,m+Lw,fLnw,m+Lnw,fLw,m+Lw,fTCw,m+Cw,fSw,m+Sw,fLnw,m+Lnw,fTSnw,m+Snw,fCnw,m+Cnw,fddt⁡u¯p¯wp¯nw=-dF¯dtq¯wq¯nw.The stiffness matrixK, coupling matrix L, flow matrix H, compressibility matrix C, saturation matrix S, equivalent nodal force vector F¯, and equivalent node flow rate vector q¯ are assembled matrixes for each element matrix labeled by the degrees of freedom, where each element matrix is expressed as(15)Kαβ=∫ΩBαTDTBβdΩα,β=a,b,c,d,Ll,m,αβ=∫ΩBαTαISl,mNβpdΩl=w,nw;α,β=a,b,c,d,Ll,f,αβ=∫ΓfNαuTRTIfnSl,fNβpdΓl=w,nw;α,β=a,b,c,d,Hl,m,αβ=∫Ω∇NαpTkmkm,rlμl∇NβpdΩl=w,nw;α,β=a,b,c,d,Hl,f,αβ=∫Γf∇NαpTRTIfτbkfkf,rlμlIfτTR∇NβpdΓl=w,nw;α,β=a,b,c,d,Cw,m,αβ=∫ΩNαpTα-ϕKsSw,mSw,m+pc,m∂Sw,m∂pc,m+ϕSw,mKw-ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Cnw,m,αβ=∫ΩNαpTα-ϕKsSnw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕSnw,mKnw-ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Cl,f,αβ=∫ΓfNαpTbSl,fKl-∂Sw,f∂pc,fNβpdΓl=w,nw;α,β=a,b,c,d,Sw,m,αβ=∫ΩNαpTα-ϕKsSw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Snw,m,αβ=∫ΩNαpTα-ϕKsSnw,mSw,m+pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Sl,f,αβ=∫ΓfNαpTb∂Sw,f∂pc,fNβpdΓl=w,nw;α,β=a,b,c,d,F¯α=∫ΓutNαuTt~dΓα=a,b,c,d,q¯l,α=∫ΓqlNαpTq~lρldΓl=w,nw;α=a,b,c,d.Equation (14) is implicitly solved in the time domain using the generalized trapezoidal method. The implicit pressure explicit saturation method (IMPES) is used to treat the multiphase fluid flow. This strong nonlinear multiphase fluid hydromechanical coupling problem can be solved by incorporating the Newmark scheme [20].


## 2.1. Strong Form

Figure2 shows a bounded control body Ω⊂Rn (n=2,3) with an outer boundary of Γ.Figure 2
Control body and boundary conditions of fractured porous media.The effective stress equilibrium equation for porous media is written as(1)∇·σ+ρg=∇·σ′-αIp+ρg=0.According to the law of mass conservation, the continuity equations for the wetting and nonwetting phase fluid flow in the matrix domain are expressed as(2)α-ϕKsSw,mSw,m+pc,m∂Sw,m∂pc,m+ϕSw,mKw-ϕ∂Sw,m∂pc,m∂pw,m∂t+α-ϕKsSw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pnw,m∂t+αSw,mIT∂ε∂t+∇·-kmkm,rwμw∇pw,m+QwstcBw=0,α-ϕKsSnw,mSw,m+pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pw,m∂t+α-ϕKsSnw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕSnw,mKnw-ϕ∂Sw,m∂pc,m∂pnw,m∂t+αSnw,mIT∂ε∂t+∇·-kmkm,rnwμnw∇pnw,m=0.Similarly, the continuity equations for the wetting and nonwetting phase fluid flow inside the fracture domain are given by(3)Sw,fKw-∂Sw,f∂pc,f∂pw,f∂t+∂Sw,f∂pc,f∂pnw,f∂t+∇·-kfkf,rwμw∇pw,f+Sw,f∇·u˙=0,∂Sw,f∂pc,f∂pw,f∂t+Snw,fKnw-∂Sw,f∂pc,f∂pnw,f∂t+∇·-kfkf,rnwμnw∇pnw,f+Snw,f∇·u˙=0.The displacement and pressure of the wetting and nonwetting phase initially are provided according to the initial conditions(4)u=u0,pw,m=pw,m0,pnw,m=pnw,m0.The displacement and pressure of the wetting and nonwetting phase at the boundary are provided by the forced boundary conditions(5)ux=u~∀x∈Γuu,pw,mx=p~w,m∀x∈Γpw,pnw,mx=p~nw,m∀x∈Γpnw.The stress and mass flow rate of the wetting and nonwetting phase at the boundary are provided by the natural boundary conditions(6)σx·n=t~∀x∈Γut,-kmkm,rwμw∇pw,m·n=q~wρw∀x∈Γqw,-kmkm,rnwμnw∇pnw,m·n=q~nwρnw∀x∈Γqnw.The fluid pressure is applied on the fracture surface(7)σx·nΓf+=Sw,fpw,f+Snw,fpnw,f·nΓf∀x∈Γf+,σx·nΓf-=-Sw,fpw,f+Snw,fpnw,f·nΓf∀x∈Γf-.The pore fluid pressure at the fracture boundary is continuous, with fluid exchange occurring between fractures and the surrounding matrix:(8)pw,m=pw,f∀x∈Γf,pnw,m=pnw,f∀x∈Γf,-kmkm,rlμl∇pl,m·nΓf+=q~l+ρll=w,nw;∀x∈Γf+,-kmkm,rlμl∇pl,m·nΓf-=q~l-ρll=w,nw;∀x∈Γf-,-kfkf,rlμl∇pl,f·nΓf+=-q~l+ρll=w,nw;∀x∈Γf+,-kfkf,rlμl∇pl,f·nΓf-=-q~l-ρll=w,nw;∀x∈Γf-.


## 2.2. Weak Form

The weak form of the governing equations can be obtained by combining the weighted residual method with the associated boundary conditions. With regard to fluid flow in fractures, since the fracture width is much smaller than its length, the fluid pressure variation along its width can be neglected (see Figure3). The areal integral inside the fracture domain Ωf can be transformed into the linear integral along the fracture surface Γf. The fluid exchange items or leak between the fractures and the surrounding matrix rock can then be vanished by adding the equivalent integrals of the matrix rock and fractures, which have the same absolute value but opposite signs. Thus, it becomes unnecessary to explicitly express the fluid exchange items or leak [16–18]:(9)∫ΩδεTσ′dΩ-∫ΩδεTαΙSw,mpw,mdΩ-∫ΩδεTαΙSnw,mpnw,mdΩ-∫ΓfδuTRTIfnSw,fpw,fdΓ-∫ΓfδuTRTIfnSnw,fpnw,fdΓ=∫ΓutδuTt~dΓ,(10)∫ΩwpTα-ϕKsSw,mSw,m+pc,m∂Sw,m∂pc,m+ϕSw,mKw-ϕ∂Sw,m∂pc,m∂pw,m∂tdΩ+∫ΩwpTα-ϕKsSw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pnw,m∂tdΩ+∫ΩwpTαSw,mIT∂ε∂tdΩ+∫Ω∇wpTkmkm,rwμw∇pw,mdΩ+∫ΓqwwpTq~wρwdΓ+∫ΓfbwpTSw,fKw-∂Sw,f∂pc,f∂pw,f∂tdΓ+∫ΓfbwpT∂Sw,f∂pc,f∂pnw,f∂tdΓ+∫ΓfwpTSw,fu˙y′dΓ+∫Γfbkfkf,rwμw∂wpT∂x′∂pw,f∂x′dΓ=0,(11)∫ΩwpTα-ϕKsSnw,mSw,m+pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,m∂pw,m∂tdΩ+∫ΩwpTα-ϕKsSnw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕSnw,mKnw-ϕ∂Sw,m∂pc,m∂pnw,m∂tdΩ+∫ΩwpTαSnw,mIT∂ε∂tdΩ+∫Ω∇wpTkmkm,rnwμnw∇pnw,mdΩ+∫ΓqnwwpTq~nwρnwdΓ+∫ΓfbwpT∂Sw,f∂pc,f∂pw,f∂tdΓ+∫ΓfbwpTSnw,fKnw-∂Sw,f∂pc,f∂pnw,f∂tdΓ+∫ΓfwpTSnw,fu˙y′dΓ+∫Γfbkfkf,rnwμnw∂wpT∂x′∂pnw,f∂x′dΓ=0.Figure 3
Schematic diagram of fluid flows in a fracture.


## 2.3. Discretization Form

XFEM introduces Heaviside step functions and branch functions to enrich the shape function, which make it easy to deal with typical strong discontinuous problems like fractures without grid matching and remeshing. The XFEM displacement field for multiple fractures is [19](12)uhx=∑i∈SNixuai+∑j=1nc∑i∈SHjubi,jNixHjx-Hjxi︸enrichment  of  bisected  support  for  fracturej+∑j=1nt∑i∈SCj∗NixRjx∑l=14uci,jlψjlx-ψjlxi︸enrichment  of  slit  support  for  fracturej+∑j=1nx∑i∈SJjudi,jNixJjx︸enrichment  of  intersecting  support  for  cross  pointj.In fractured porous media, the natural fractures are the main channels for seepage. The pore fluid pressure field is continuous at the fracture surface, while the derivative (i.e., flow rate) is discontinuous. This presents a typical weak discontinuous problem. Similar to the interpolation of displacement field, the pressure field can be interpolated as [16, 18](13)phx=∑i∈SNixpai+∑j=1nc∑i∈SHjpbi,jNixfjx-fjxi︸enrichment  of  bisected  support  for  fracturej+∑j=1nt∑i∈SCj∗NixRjx∑l=13pci,jlBjlx-Bjlxi︸enrichment  of  slit  support  for  fracturej+∑j=1nx∑i∈SJjpdi,jNixJjpx-Jjpxi︸enrichment  of  intersecting  support  for  cross  pointj.The discretization form of the fully coupled hydromechanical model for fractured porous media can be obtained by regarding the trial functionsδu and wp in (9)–(11) as the variation of displacement and pressure interpolation functions, as shown in(14)0000Hw,m+Hw,f000Hnw,m+Hnw,fu¯p¯wp¯nw+-KLw,m+Lw,fLnw,m+Lnw,fLw,m+Lw,fTCw,m+Cw,fSw,m+Sw,fLnw,m+Lnw,fTSnw,m+Snw,fCnw,m+Cnw,fddt⁡u¯p¯wp¯nw=-dF¯dtq¯wq¯nw.The stiffness matrixK, coupling matrix L, flow matrix H, compressibility matrix C, saturation matrix S, equivalent nodal force vector F¯, and equivalent node flow rate vector q¯ are assembled matrixes for each element matrix labeled by the degrees of freedom, where each element matrix is expressed as(15)Kαβ=∫ΩBαTDTBβdΩα,β=a,b,c,d,Ll,m,αβ=∫ΩBαTαISl,mNβpdΩl=w,nw;α,β=a,b,c,d,Ll,f,αβ=∫ΓfNαuTRTIfnSl,fNβpdΓl=w,nw;α,β=a,b,c,d,Hl,m,αβ=∫Ω∇NαpTkmkm,rlμl∇NβpdΩl=w,nw;α,β=a,b,c,d,Hl,f,αβ=∫Γf∇NαpTRTIfτbkfkf,rlμlIfτTR∇NβpdΓl=w,nw;α,β=a,b,c,d,Cw,m,αβ=∫ΩNαpTα-ϕKsSw,mSw,m+pc,m∂Sw,m∂pc,m+ϕSw,mKw-ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Cnw,m,αβ=∫ΩNαpTα-ϕKsSnw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕSnw,mKnw-ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Cl,f,αβ=∫ΓfNαpTbSl,fKl-∂Sw,f∂pc,fNβpdΓl=w,nw;α,β=a,b,c,d,Sw,m,αβ=∫ΩNαpTα-ϕKsSw,mSnw,m-pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Snw,m,αβ=∫ΩNαpTα-ϕKsSnw,mSw,m+pc,m∂Sw,m∂pc,m+ϕ∂Sw,m∂pc,mNβpdΩα,β=a,b,c,d,Sl,f,αβ=∫ΓfNαpTb∂Sw,f∂pc,fNβpdΓl=w,nw;α,β=a,b,c,d,F¯α=∫ΓutNαuTt~dΓα=a,b,c,d,q¯l,α=∫ΓqlNαpTq~lρldΓl=w,nw;α=a,b,c,d.Equation (14) is implicitly solved in the time domain using the generalized trapezoidal method. The implicit pressure explicit saturation method (IMPES) is used to treat the multiphase fluid flow. This strong nonlinear multiphase fluid hydromechanical coupling problem can be solved by incorporating the Newmark scheme [20].


## 3. Numerical Simulation

### 3.1. Simultaneous Fracturing for Horizontal Wells

The cluster perforation and simultaneous fracturing are key technologies for volume stimulation of horizontal wells. The propped hydraulic fracture induces additional stresses in the surrounding rock. This stress interference between fractures can complicate the propagation of hydraulic fractures to form complex fracture networks. Crosby et al. [21] and Rahman et al. [22] investigated the propagation behavior of two isolated hydraulic fractures induced by the simultaneous fracturing of a horizontal well.This study simulated Crosby’s DC-8 specimen, as Figure4 shows. Modeling was performed using material properties and injection parameters similar to those employed during the laboratory test DC-8. The specimen was meshed into uniform 80 × 80 structured bilinear quadrilateral elements, with the propagation step set at a constant value of 0.01 m. The coupling process was neglected during this simulation because the synthetic fracture test material was cast by silica flour and acrylic emulsion to provide the specimen with a negligible leak coefficient, 3.8 × 10−8 m/s, and an ultralow permeability, 8 × 10−5 mD. Therefore, only the first item in (14) was used. Equation (9) states that the fluid pressure inside the fracture, pf, should be calculated, which presents a very complicated problem, especially for fractures that are not straight. Moreover, the fluid pressure inside the fracture must be transferred into an equivalent nodal force using the linear integral form, which requires integral points on the fracture surface [23]. We set the net pressure inside the fractures at 3 MPa.Figure 4
Dimensions of the artificial rock specimen.The maximum circumferential stress criterion is selected to determine the propagation condition and crack direction. The stress intensity factors of hydraulic fracture were calculated with the interaction integral(16)I1,2=∫Aσij1ui,12+σij2ui,11-W1,2δ1jq,jdA-∫C+∪C-σi21ui,12+σi22ui,11qm2dΓ.The first item of the right side can be calculated easily. The second item of the right side was derived explicitly:(17)I21,I=-KI2κ+12μ2π∫0r1pfnm2rdΓ+1r1-r2∫r1r2r-r2pfnm2rdΓ,I21,II=-KII2κ+12μ2π∫0r1pfτm2rdΓ+1r1-r2∫r1r2r-r2pfτm2rdΓ.Figure5 shows comparisons between the XFEM simulation results, experimental results, and HYFRANC finite element simulation results.Figure 5
Comparison between the experimental results and XFEM numerical simulation results of the staged fracturing of a horizontal well.The results of XFEM simulation are more consistent with the experimental result in the upper fracture than in the lower fracture. This can be largely attributed to the fact that the lower fracture was initiated preferentially during the experiment. Therefore, the propagation path was relatively straight. However, the two fractures were considered to initiate at the same time in the simulation. Then, the lower fracture was forced to deflect towards the wellbore for stress shadow effect. In comparison, the HYFRANC simulation results [21] and the XFEM simulation results show high consistency with the experimental results. Moreover, the simulation efficiency improved significantly.The stress distribution at the sixth propagation step (Figure6) shows that the x-axis stress, σxx, between the two fractures increased from 5 MPa to about 7.62 MPa (compressive), the y-axis stress, σyy, increased from 7 MPa to about 8.06 MPa (compressive), and the shear stress, σxy, was lower than 0.0002 MPa. In other words, the stress state between the two fractures changed significantly. The maximum horizontal nominal stress direction changed from the x-axis to the y-axis. Perforation fracturing between two fractures tends to cause the hydraulic fracture to propagate along the y-axis, intersect with other fractures, and form a fracture network.Figure 6
The displacement distribution of the sixth propagation step.
(a)
x-axis stress σxx (MPa)
(b)
y-axis stress σyy (MPa)
(c)
x
-
y planar shearing stress σxy (0.1 kPa)


### 3.2. Coupled Hydromechanical Behavior for Hydraulic Fracturing

Figure7 shows a 2D fractured porous media control domain. The physical model has a length and width of 150 m, which contains three fractures. Table 1 lists detailed parameters for the model. The hydraulic fracture, which is parallel to the global x-axis, has a length l1 of 50 m and an initial width of 1 mm. The orthogonal natural fracture, which is parallel to the global y-axis, has a length l2 of 40 m and an initial width of 0.1 mm. The nonorthogonal natural fracture, which lies at an included angle of 45° to the global x-axis, has a length l3 of 30 m and an initial width of 0.1 mm. The initial stress is zero, while the pore fluid pressure inside both the fractures and the matrix rock is 50 MPa. The displacements of the upper, lower, right, and left sides are all restricted (ux=uy=0), and they are therefore impermeable boundaries. In reality, the length, width, height, and net pressure would be increased with increasing injection time, and the stress intensity factors and induced stresses would be strengthened to force the deformation of natural fracture, subsequently. For simplicity, the fractures were presumed as nonpropagating cracks to investigate the deformation behavior between hydraulic and natural fractures before intersection. However, the net pressure and width of fractures were solved with the fully coupled hydromechanical model. Although this assumption was not coinciding with the reality, the numerical simulation can represent the deformation behavior between hydraulic and natural fractures to some extent.Table 1
Model parameters for the seepage-stress coupling analysis.
Item
Value
Young’s modulus,E (GPa)
30
Poisson’s ratio,v
0.25
Biot coefficient,αm
0.85
Rock permeability,km (mD)
0.1
Rock porosity,ϕm
0.08
Bulk modulus of rock,Ks (GPa)
134.6
Coefficient between real and ideal fractures,f
1.3
Bulk modulus of wetting phase,Kw (GPa)
2.2
Viscosity of wetting phase,μw (mPa⋅s)
1
Liquid injection rate,Qw (m3/min)
10
Reservoir thickness,Hf (m)
30
Density of wetting phase,ρw (kg/m3)
1000
Volume factor of wetting phase,Bw
1.05Figure 7
Control domain of fractured porous media for the hydromechanical coupling analysis.A uniform 80 × 79 bilinear quadrilateral mesh was used for calculation, with the intention of locating the injection well point at the element boundary, while the fracture location was independent of the generated mesh.Thex-axis displacement, ux, distribution is shown in Figure 8. The y-axis displacement, uy, distribution is shown in Figure 9. The displacement fields demonstrate that the injected fluid initially flows along the first fracture, forcing the fluid pressure inside the fracture to increase. Meanwhile, the fluid filtrated into the surrounding matrix rock, leading to an increase of pore fluid pressure in the matrix rock. The increase in fluid pressure drives the deformation of fractured porous media, generating displacement along the x-axis and the y-axis. Since the contact condition is not considered here, parts of the other fracture surfaces may come into contact and intersect with each other, while other parts become staggered.Figure 8
x-axis displacement at different moments (mm).
(a)
t = 60 s
(b)
t = 600 sFigure 9
y-axis displacement at different moments (mm).
(a)
t = 60 s
(b)
t = 600 sFracture width represents the key parameter in hydraulic fractures. Equation (12) allows the displacement of the two surfaces of a fracture at any position to be obtained. The fracture width can be calculated by subtracting the positive-side displacement from the negative-side displacement. The width decreases with increasing distance from the injecting well point, and it increases with increasing injection duration, as Figure 10 shows. The fracture width is distributed symmetrically. The orthogonal natural fracture remains closed at first under the effect of stress induced by the left crack tip of a hydraulic fracture. As the fluid pressure inside the hydraulic fracture increases, the crack-tip-induced stress increases and forces the center section to open, while the faraway section remains closed, as shown in Figure 11. This is one of the reasons why a hydraulic fracture intersects a natural fracture at a large intersection angle. In contrast, the stress induced by a hydraulic fracture drives the nonorthogonal natural fracture to open. As Figure 12 shows, only a small part of the initial fracture remains closed, while the whole fracture opens finally.Figure 10
Width section of the hydraulic fracture.Figure 11
Width section of the orthogonal natural fracture.Figure 12
Width section of the nonorthogonal natural fracture.


## 3.1. Simultaneous Fracturing for Horizontal Wells

The cluster perforation and simultaneous fracturing are key technologies for volume stimulation of horizontal wells. The propped hydraulic fracture induces additional stresses in the surrounding rock. This stress interference between fractures can complicate the propagation of hydraulic fractures to form complex fracture networks. Crosby et al. [21] and Rahman et al. [22] investigated the propagation behavior of two isolated hydraulic fractures induced by the simultaneous fracturing of a horizontal well.This study simulated Crosby’s DC-8 specimen, as Figure4 shows. Modeling was performed using material properties and injection parameters similar to those employed during the laboratory test DC-8. The specimen was meshed into uniform 80 × 80 structured bilinear quadrilateral elements, with the propagation step set at a constant value of 0.01 m. The coupling process was neglected during this simulation because the synthetic fracture test material was cast by silica flour and acrylic emulsion to provide the specimen with a negligible leak coefficient, 3.8 × 10−8 m/s, and an ultralow permeability, 8 × 10−5 mD. Therefore, only the first item in (14) was used. Equation (9) states that the fluid pressure inside the fracture, pf, should be calculated, which presents a very complicated problem, especially for fractures that are not straight. Moreover, the fluid pressure inside the fracture must be transferred into an equivalent nodal force using the linear integral form, which requires integral points on the fracture surface [23]. We set the net pressure inside the fractures at 3 MPa.Figure 4
Dimensions of the artificial rock specimen.The maximum circumferential stress criterion is selected to determine the propagation condition and crack direction. The stress intensity factors of hydraulic fracture were calculated with the interaction integral(16)I1,2=∫Aσij1ui,12+σij2ui,11-W1,2δ1jq,jdA-∫C+∪C-σi21ui,12+σi22ui,11qm2dΓ.The first item of the right side can be calculated easily. The second item of the right side was derived explicitly:(17)I21,I=-KI2κ+12μ2π∫0r1pfnm2rdΓ+1r1-r2∫r1r2r-r2pfnm2rdΓ,I21,II=-KII2κ+12μ2π∫0r1pfτm2rdΓ+1r1-r2∫r1r2r-r2pfτm2rdΓ.Figure5 shows comparisons between the XFEM simulation results, experimental results, and HYFRANC finite element simulation results.Figure 5
Comparison between the experimental results and XFEM numerical simulation results of the staged fracturing of a horizontal well.The results of XFEM simulation are more consistent with the experimental result in the upper fracture than in the lower fracture. This can be largely attributed to the fact that the lower fracture was initiated preferentially during the experiment. Therefore, the propagation path was relatively straight. However, the two fractures were considered to initiate at the same time in the simulation. Then, the lower fracture was forced to deflect towards the wellbore for stress shadow effect. In comparison, the HYFRANC simulation results [21] and the XFEM simulation results show high consistency with the experimental results. Moreover, the simulation efficiency improved significantly.The stress distribution at the sixth propagation step (Figure6) shows that the x-axis stress, σxx, between the two fractures increased from 5 MPa to about 7.62 MPa (compressive), the y-axis stress, σyy, increased from 7 MPa to about 8.06 MPa (compressive), and the shear stress, σxy, was lower than 0.0002 MPa. In other words, the stress state between the two fractures changed significantly. The maximum horizontal nominal stress direction changed from the x-axis to the y-axis. Perforation fracturing between two fractures tends to cause the hydraulic fracture to propagate along the y-axis, intersect with other fractures, and form a fracture network.Figure 6
The displacement distribution of the sixth propagation step.
(a)
x-axis stress σxx (MPa)
(b)
y-axis stress σyy (MPa)
(c)
x
-
y planar shearing stress σxy (0.1 kPa)


## 3.2. Coupled Hydromechanical Behavior for Hydraulic Fracturing

Figure7 shows a 2D fractured porous media control domain. The physical model has a length and width of 150 m, which contains three fractures. Table 1 lists detailed parameters for the model. The hydraulic fracture, which is parallel to the global x-axis, has a length l1 of 50 m and an initial width of 1 mm. The orthogonal natural fracture, which is parallel to the global y-axis, has a length l2 of 40 m and an initial width of 0.1 mm. The nonorthogonal natural fracture, which lies at an included angle of 45° to the global x-axis, has a length l3 of 30 m and an initial width of 0.1 mm. The initial stress is zero, while the pore fluid pressure inside both the fractures and the matrix rock is 50 MPa. The displacements of the upper, lower, right, and left sides are all restricted (ux=uy=0), and they are therefore impermeable boundaries. In reality, the length, width, height, and net pressure would be increased with increasing injection time, and the stress intensity factors and induced stresses would be strengthened to force the deformation of natural fracture, subsequently. For simplicity, the fractures were presumed as nonpropagating cracks to investigate the deformation behavior between hydraulic and natural fractures before intersection. However, the net pressure and width of fractures were solved with the fully coupled hydromechanical model. Although this assumption was not coinciding with the reality, the numerical simulation can represent the deformation behavior between hydraulic and natural fractures to some extent.Table 1
Model parameters for the seepage-stress coupling analysis.
Item
Value
Young’s modulus,E (GPa)
30
Poisson’s ratio,v
0.25
Biot coefficient,αm
0.85
Rock permeability,km (mD)
0.1
Rock porosity,ϕm
0.08
Bulk modulus of rock,Ks (GPa)
134.6
Coefficient between real and ideal fractures,f
1.3
Bulk modulus of wetting phase,Kw (GPa)
2.2
Viscosity of wetting phase,μw (mPa⋅s)
1
Liquid injection rate,Qw (m3/min)
10
Reservoir thickness,Hf (m)
30
Density of wetting phase,ρw (kg/m3)
1000
Volume factor of wetting phase,Bw
1.05Figure 7
Control domain of fractured porous media for the hydromechanical coupling analysis.A uniform 80 × 79 bilinear quadrilateral mesh was used for calculation, with the intention of locating the injection well point at the element boundary, while the fracture location was independent of the generated mesh.Thex-axis displacement, ux, distribution is shown in Figure 8. The y-axis displacement, uy, distribution is shown in Figure 9. The displacement fields demonstrate that the injected fluid initially flows along the first fracture, forcing the fluid pressure inside the fracture to increase. Meanwhile, the fluid filtrated into the surrounding matrix rock, leading to an increase of pore fluid pressure in the matrix rock. The increase in fluid pressure drives the deformation of fractured porous media, generating displacement along the x-axis and the y-axis. Since the contact condition is not considered here, parts of the other fracture surfaces may come into contact and intersect with each other, while other parts become staggered.Figure 8
x-axis displacement at different moments (mm).
(a)
t = 60 s
(b)
t = 600 sFigure 9
y-axis displacement at different moments (mm).
(a)
t = 60 s
(b)
t = 600 sFracture width represents the key parameter in hydraulic fractures. Equation (12) allows the displacement of the two surfaces of a fracture at any position to be obtained. The fracture width can be calculated by subtracting the positive-side displacement from the negative-side displacement. The width decreases with increasing distance from the injecting well point, and it increases with increasing injection duration, as Figure 10 shows. The fracture width is distributed symmetrically. The orthogonal natural fracture remains closed at first under the effect of stress induced by the left crack tip of a hydraulic fracture. As the fluid pressure inside the hydraulic fracture increases, the crack-tip-induced stress increases and forces the center section to open, while the faraway section remains closed, as shown in Figure 11. This is one of the reasons why a hydraulic fracture intersects a natural fracture at a large intersection angle. In contrast, the stress induced by a hydraulic fracture drives the nonorthogonal natural fracture to open. As Figure 12 shows, only a small part of the initial fracture remains closed, while the whole fracture opens finally.Figure 10
Width section of the hydraulic fracture.Figure 11
Width section of the orthogonal natural fracture.Figure 12
Width section of the nonorthogonal natural fracture.


## 4. Conclusion

Hydraulic fracturing is an interactive process incorporating both mechanical and hydraulic effects. The hydromechanical coupling model and associated solution scheme proposed under the framework of XFEM can effectively avoid tedious calculations of both the fluid pressure distribution in fractures and the transformation of equivalent nodal forces. This model also considers the influence of pore pressure variation in the rock matrix on fracture propagation.The mechanical module can simulate the fracture propagation pathway of the simultaneous fracturing for a horizontal well. Under the influence of stress perturbation between hydraulic fractures, not the straight transverse fractures but the flexural fractures deflecting toward the wellbore were created for the staged fracturing of a horizontal well. The stress shadow effect introduced by the propped hydraulic fracture reduces the horizontal stress differential between fractures or even leads to stress inversion. Therefore, it is beneficial to enhance the complexity of fracture networks.The induced stresses added to the crack tip of the hydraulic fracture keep most of the orthogonally intersected natural fractures closed, with only a small part opening due to the induced stress and fracturing fluid leakage. The stress concentration at the crack tip of the hydraulic fracture fully opens most of the nonorthogonal intersected natural fractures as time passes. The hydraulic fracture tends to traverse the orthogonal natural fracture and continue to propagate, while it is more easily arrested by the nonorthogonal natural fracture.



---

*Source: 6373957-2017-06-07.xml*