# Numerical Simulation for Fracture Propagation in Elastoplastic Formations

Authors: Yafei Hu; Jin Zhao; Lihu Cao; Jinzhou Zhao; Junshi Li; Zhiying Wu; Jianfeng Hou
Journal: Geofluids
 (2021)

---

## Abstract

Current hydraulic fracture models are mainly based on elastic theories, which fail to give accurate prediction of fracture parameters in plasticity formation. This paper proposes a fluid–solid coupling model for fracture propagation in elastoplastic formations. The rock plastic deformation in the model satisfies the Mohr-Coulomb yield criterion and plastic strain increment theory. The extended finite-element method (XFEM) combined with the cohesive zone method (CZM) is used to solve the coupled model. The accuracy of the model is validated against existing models. The effects of stress difference, friction angle, and dilation angle on fracture shape (length, width), injection pressure, plastic deformation, induced stress, and pore pressure are investigated through the model. The results indicate that compared with elastic formation, fracture propagation in elastoplastic formation is more difficult, the breakdown pressure and extending pressure are greater, and fracture shape is wider and shorter. The plastic deformation causes the fracture tip to become blunt. Under the condition of high stress difference or low friction angle formation, it is prone to occur large plastic deformation zones and form wide and short fracture. Compared with friction angle, dilation angle is less sensitive to plastic deformation, fracture parameters, and fracture geometry. For the formation with high stress difference and friction angle, the effect of plasticity deformation on fracture propagation should not be ignored.


---

## Body

## 1. Introduction

With the increasing energy consumption, efficient development of oil and gas resources has become the focus of attention. Hydraulic fracturing (HF), as a widely used technology, has become increasingly important in oil and gas stimulation [1, 2]. Studying fracture propagation mechanism has a great significance for fracture parameter (width, length, height) optimization. The classical fracture propagation models treat formation as an elastic medium without considering the plasticity effects and fluid-solid coupling [3–8]. However, plastic failure usually occurs in soft layer or unconsolidated formation with high temperature and pressure during fracturing [9–11].Various experiments have been conducted to study fracture propagation in soft rocks. The results indicate that with the increase of confining pressure, the rock property gradually appears plastic character, which has great influence on fracture pressure and fracture shape [12–16]. Complex fractures and fluid lag near the fracture can be also observed in soft rock during fracturing. Initiation pressure and fracture width are also affected by the confining stresses in unconsolidated formation [17–21]. Owing to rock plasticity and strong fluid-solid coupling, traditional fracture propagation models fail to accurately predict the fracture parameters in soft formation.Previous studies on the effect of plasticity in HF were just based on static fracture and constant injection pressure in fracture without considering the coupling between the pressure diffusion and rock deformation. Papanastasiou [22–24] proposed that plastic yielding was affected by rock dilation, stress difference, and the cohesive zone and utilized the Mohr-Coulomb criterion to study the effect of pore pressure and permeability on fracture parameters without considering the fluid–solid coupling in porous formation. Zaki et al. [25] used the cam-clay model to assess the influence of rock plasticity on fracture volumes based on the simple models that would rely in using lower modulus or higher fracture toughness to account for plasticity. Wang et al. [26] proposed the poroelastic HF model in brittle and ductile formations. Liu et al. [27] established fracture propagation model based on the Drucker-Prager yield criterion. Lin et al. [28] investigated the effect of Young’s model, Poisson’s ratio, and injection rate on fracture propagation based on the modified cam-clay model. However, the influence of rock plasticity on the pore pressure and induced stress distribution and the effect of rock friction angle, dilation angle, and principal stress difference on fracture propagation based on fluid–solid coupling have not been considered in the above studies. Therefore, it is necessary to establish a fluid-solid coupling elastoplastic model for fracture propagation.In this study, a fluid–solid coupling elastoplasticity HF model is proposed. The rock plastic deformation conforms to the Mohr-Coulomb criterion and plastic strain increment theory. Fluid flow within the fracture follows the mass conservation law. The fracture propagation criterion is based on the traction/separation law of cohesive element. By combining the extended finite-element method (XFEM) and the cohesive zone method (CZM), the Abaqus 6.14 software has been utilized to simulate fracture initiation and growth in elastoplastic formation. The simulation results are also compared with existing fracture expansion models. Meanwhile, the influences of principal stress difference, friction angle, and dilation angle on fracture length, width, injection pressure, plastic deformation, pore pressure, and induced stress have been further investigated.


## 2. Mathematical Model

Fracture propagation involves rock deformation, rock failure, fluid flowing in fracture and formation, the dynamic variation of pore pressure, and induced stress, which are a very complex process and involve multiple coexisting and interdependent subprocesses [29]. In addition, formation heterogeneity and plasticity make modeling more complex. In order to simplify the process of crack propagation, the following assumptions are made:(1)
Formation is a homogeneous, porous elastic-plastic rock medium(2)
Fracture extending without proppant transport(3)
Fluid flows within the fracture existing fluid leakoff(4)
Fluid flows through the formation in accordance with the Darcy law


### 2.1. Fluid-Solid Coupling Equation of Rock Deformation

Based on the above assumptions, rock deformation conforms to pore-elastoplastic theory. According to the stress balance equation, the stress equilibrium, strain displacement, and strain–stress equations can be written as follows [30]:(1)dσij,j+dfi=0,dεij=12dui,j+duj,i,dσij=Dijklepdεkl,where dσij is the stress increment, dεij is the strain increment, dui,j is the medium displacement increment, dfi is the body force increment, and Depijkl is the elastic-plastic coefficient matrix.The displacement and stress at the outer boundary and fluid pressure at the fracture surface are [31](2)du=du∗,onΓu,dσ•n=dt∗,onΓt,dnTσ•n=−dp,onΓcr,where Γu, Γt, and Γcr are the displacement boundary condition, stress boundary condition, and fracture pressure boundary condition, respectively; u∗ and t∗ are the corresponding displacement and stress increment on the boundary; and dp is fluid pressure increment on the fracture.


### 2.2. The Elastic-Plastic Increment Constitutive Model

When the rock begins to occur plastic failure, the strain increment of the rock is composed of elastic strain increment and plastic strain [32]:(3)dεij=dεije+dεijp.The elastic strain incrementdεeij and the stress increment dσij satisfy the Hooke’s law:(4)dσij=Dijkledεkle.The plastic strain increment follows the flow rule, which is defined as [33]:(5)dεijp=λ∂F∂σij.In the process of plastic deformation, the yield surface of the rock is a function of strengthening parameters, plastic strain, and stress [32]:(6)F=Fσ,ijεijp,κ,(7)dF=∂F∂σijdσij+∂F∂κdκ+∂F∂εijpdεijp=0,where εij is the total strain, εeIj is the elastic strain, εpIj is the plastic strain, σij is the stress tensor, Deijkl is the elastic constitutive tensor, F is the load surface equation, κ is the strengthening parameter of rock, and λ is a quantity related to the slope of the uniaxial curve of stress σ-plastic strain εp, which can be determined by yield criteria and strengthening theory.The stress increment is obtained by substituting equation (4) into (3).(8)dσij=Dijkledεkl‐dεklp.Multiplying both sides of equation (8) by ∂F/∂σij and substituting it into equation (7):(9)∂F∂σijDijkledεkl‐dεklp=‐∂F∂κdκ‐∂F∂εijpdεijp.The hardening parameter is a function of the equivalent plastic strain and temperature, which is determined as [32]:(10)∂F∂κdκ=∂F∂κ∂κ∂ε¯pdε¯p+∂κ∂TdT=∂F∂κ∂κ∂ε¯p+∂κ∂T∂T∂ε¯pdε¯p=∂F∂κ∂κ∂ε¯p+∂κ∂T∂T∂ε¯pdε¯pdε¯ijPdε¯ijP.The equivalent plastic strain increment is defined as follows:(11)dε¯p=23dεijPdεijP1/2=23dεxP−dεyP2+dεyP−dεzP2+dεzP−dεxP2+32dγxyP2+dγyzP2+dγxzP21/2,where −εP is the equivalent plastic strain; T is the temperature; εpx, εpy, and εpz are the plastic strain on the x, y and z axis, respectively; and γpxy, γpyz, and γpxz are the shear stress component.Substituting equation (10) and equation (5) into equation (9).(12)λ=∂F/∂σijDijkledεklA,where(13)A=∂F∂σijDijkle−∂F∂κ∂κ∂ε¯p+∂κ∂T∂T∂ε¯pdε¯pdε¯ijP‐∂F∂εklp∂F∂σkl.Substituting equations (12) and (13) into equation (3) and then substituting equation (5) into equation (8), the elastic-plastic matrix can be obtained.(14)dσij=Dijklepdεkl,(15)Dijklep=Dijkle‐Dijkle∂F/∂σpq∂F/∂σrsDpqrseA.


### 2.3. Plastic Yield Criteria for Rocks

Rock deformation obeys the Mohr-Coulomb yield criterion, and we stipulate that the tensile stress is positive and the compressive stress is negative [22]:(16)τ=c−σtanϕ.According to the geometry of the mole circle, shear stressτ and normal stress σ can be expressed as follows [34]:(17)τ=σ1−σ32cosϕ,σ=σ1+σ32+σ1−σ32sinϕ,where τ is the shear stress, σ is the normal stress, φ is the friction degree, σ1 and σ3 are the maximum and minimum principal stresses, respectively, and σ2 is the intermediate principal stress.Substitute equation (17) into equation (16), we have(18)σ1−σ32cosϕ=c−σ1+σ32+σ1−σ32sinϕtanϕ.Principal stress and principal deviatoric stress satisfy the following equation [32]:(19)σ1=23J2cosΘ+13I1,σ2=23J2cos2π3‐Θ+13I1,σ3=23J2cos2π3+Θ+13I1,where σ1 and σ3 are the maximum and minimum principal stresses, respectively, σ2 is the intermediate principal stress, I1 is the first invariant of stress tensor, J2 is the second invariant of deviatoric stress, and Θ is the deviatoric polar angle.Substitute equation (19) into equation (18):(20)F=13I1sinϕ+J2sinΘ+π3+J23cosΘ+π3sinϕ−ccosϕ=0.Equivalent compressive stressp, Mises equivalent stress q, the first invariant of stress tensor I1, and the second invariant of deviatoric stress J2 satisfy the following equation [25]:(21)p=‐13traceσ==‐I13,q=32sijsij=3J2,cos3Θ=r3q3,where p is the equivalent pressure stress, q is the Mises equivalent stress, r is the third invariant of deviatoric stress, and s is the deviatoric stress.Substitute equation (21) into equation (20):(22)F=fp,q,Θ=Rmcq−ptanϕ−c=0,where(23)Rmc=13cosϕsinΘ+π3+13cosΘ+π3tanϕ.In the case of equation (23), there will be sharp angles on the yield surface, leading to not unique plastic flow direction, which leads to tedious calculation and slow convergence. In order to solve this problem, Menetrey and Willam [35] proposed a continuous smooth elliptic function as a plastic potential surface.(24)G=εc0tanψ2+Rmwq2−ptanψ,where(25)Rmw=41−e2cos2Θ+2e−1221−e2cosΘ+2e−141−e2cos2Θ+5e2−4e3−sinϕ6cosϕ,(26)e=3−sinϕ3+sinϕ,where ψ is the dilation angle measured in the p−Rmwq plane at high confining pressure and can depend on temperature and predefined field variables; c0 is the initial cohesion yield stress; ε is the eccentricity ratio, and here is assumed to be 0.1; e is the “out of roundedness” of the deviatoric section in terms of the ratio between the shear stress along the extension meridian and the shear stress along the compression meridional.


### 2.4. Fluid Flow in Fracture and Formation

Once the fracture is formed, fluid will flow into the fracture at once. Fluid flow in the fracture follows the mass conservation law [36]:(27)∂w∂t+∂qt∂s+vt+vb=0,where qt is the injection flow rate, w is the fracture width, s is the fracture length, t is injection time, and vt and vb are the fracturing fluid leakoff rate through the top and bottom fracture surfaces, respectively.Both tangential flow along and normal flow across the fracture surfaces can be calculated using the Poiseuille equation [37, 38]:(28)qt=w312μ∂p∂s,where μ is the fluid viscosity and p is the fluid pressure.Fluid leakoff rate is mainly affected by pressure and leakoff coefficient, which satisfy the following equation [39, 40]:(29)vt=ctp−ptop,vb=ctp−pbot,where ct is the fluid leakoff coefficients and ptop and pbot are the pressure in the top and bottom fracture surface, respectively.When the fluid leaks off from the fracture, it will flow through the formation in accordance with the Darcy law. The mass conservation equation is expressed as [34]:(30)∇ρwkwbwμw∇pw+qw=∂∂tφρwswbw,where ρw is the water permeability, kw is the water permeability, pw is water pressure in formation, bw is the water compression coefficient, μw is the water viscosity, qw is the leakoff fluid, Sw is the water saturation, and ϕ is the effective porosity.


### 2.5. Cohesive Crack Model

The singularity of the fracture tip stress field will affect the convergence of the calculation results. In order to solve this problem, CZM is adopted to simulate the initiation and propatation of fractures. This method is implemented by embedding an artificially predefined path of cohesive element into the formation. The cohesive element glues two conjugated cohesive surfaces of the reservoir pressure/displacement element. The failure of these element obeys the traction/separation law in the vertical or shear direction [41, 42]. When the cohesive element is completely destroyed, the two bound cohesive surfaces are separated, and hydraulic fracture is generated. When the maximum nominal stress criterion is adopted for rock failure, initiation propagation occurs as the maximum nominal stress ratio reaches a critical value [43–45]:(31)fc=maxTnTn0,TsTs0,(32)1≤fc≤1+ft,where δn and δs are, respectively, the normal and tangential displacement jump vector, δeq are effective displacement, Tn0 and Ts0 are, respectively, the normal and tangential stress component of the stress element under a linear elastic condition, < > denotes the Macaulay bracket, fc is a parameter evaluating the fracture failure, and ft is assured to be 0.05.The fracture damage evolution is determined by damage variables, and the expressions of damage effects on the normal and tangential stress components of the element can be obtained as [44–46](33)Tn=1−dTn0,Tn0≥0Tn0,Tn0<0,(34)Ts=1−dTs0,where d is the scalar damage variable between fractures, 0<d<1, and Tn and Ts are, respectively, the values of the normal stress component and the tangential stress component of the traction/separation model.When the critical fracture energy of rock material along the first and second shear directions is similar, the damage evolution during fracture propagation can also be determined by the Benzeggagh–Kenane fracture criterion [47]. The combined energy dissipated by failure Gc is defined as(35)Gc=GIc+GIIc−GIcGIIc−GIcη,where Gc is the critical fracture energy release rate, GcI and GcII are the normal and first shear direction fracture toughness values, GI and GII are the normal and shear direction fracture energy release rate, and η is a material properties parameter; here is assumed to be 2.3.


### 2.6. Model Solution

The iterative coupling approach is used to solve the coupled model. Fluid flow and solid deformation are solved separately and sequentially, and the coupling terms are iterated until convergence. The fluid flow model is solved first based on the initial guess of injection pressure; then, the solid deformation is solved by the XFEM stiffness matrix, which is based on the solution of the previous iteration. The pressure distribution is calculated until convergence is reached. If the convergence criterion is not satisfied, the guessing value is modified aspn+1=ωpn−1+1−ωpn; the pressure distribution and rock deformation will be recalculated with the modified value. When the rock failure condition is met, the fracture will propagate and add a new element for the next time step.


## 2.1. Fluid-Solid Coupling Equation of Rock Deformation

Based on the above assumptions, rock deformation conforms to pore-elastoplastic theory. According to the stress balance equation, the stress equilibrium, strain displacement, and strain–stress equations can be written as follows [30]:(1)dσij,j+dfi=0,dεij=12dui,j+duj,i,dσij=Dijklepdεkl,where dσij is the stress increment, dεij is the strain increment, dui,j is the medium displacement increment, dfi is the body force increment, and Depijkl is the elastic-plastic coefficient matrix.The displacement and stress at the outer boundary and fluid pressure at the fracture surface are [31](2)du=du∗,onΓu,dσ•n=dt∗,onΓt,dnTσ•n=−dp,onΓcr,where Γu, Γt, and Γcr are the displacement boundary condition, stress boundary condition, and fracture pressure boundary condition, respectively; u∗ and t∗ are the corresponding displacement and stress increment on the boundary; and dp is fluid pressure increment on the fracture.


## 2.2. The Elastic-Plastic Increment Constitutive Model

When the rock begins to occur plastic failure, the strain increment of the rock is composed of elastic strain increment and plastic strain [32]:(3)dεij=dεije+dεijp.The elastic strain incrementdεeij and the stress increment dσij satisfy the Hooke’s law:(4)dσij=Dijkledεkle.The plastic strain increment follows the flow rule, which is defined as [33]:(5)dεijp=λ∂F∂σij.In the process of plastic deformation, the yield surface of the rock is a function of strengthening parameters, plastic strain, and stress [32]:(6)F=Fσ,ijεijp,κ,(7)dF=∂F∂σijdσij+∂F∂κdκ+∂F∂εijpdεijp=0,where εij is the total strain, εeIj is the elastic strain, εpIj is the plastic strain, σij is the stress tensor, Deijkl is the elastic constitutive tensor, F is the load surface equation, κ is the strengthening parameter of rock, and λ is a quantity related to the slope of the uniaxial curve of stress σ-plastic strain εp, which can be determined by yield criteria and strengthening theory.The stress increment is obtained by substituting equation (4) into (3).(8)dσij=Dijkledεkl‐dεklp.Multiplying both sides of equation (8) by ∂F/∂σij and substituting it into equation (7):(9)∂F∂σijDijkledεkl‐dεklp=‐∂F∂κdκ‐∂F∂εijpdεijp.The hardening parameter is a function of the equivalent plastic strain and temperature, which is determined as [32]:(10)∂F∂κdκ=∂F∂κ∂κ∂ε¯pdε¯p+∂κ∂TdT=∂F∂κ∂κ∂ε¯p+∂κ∂T∂T∂ε¯pdε¯p=∂F∂κ∂κ∂ε¯p+∂κ∂T∂T∂ε¯pdε¯pdε¯ijPdε¯ijP.The equivalent plastic strain increment is defined as follows:(11)dε¯p=23dεijPdεijP1/2=23dεxP−dεyP2+dεyP−dεzP2+dεzP−dεxP2+32dγxyP2+dγyzP2+dγxzP21/2,where −εP is the equivalent plastic strain; T is the temperature; εpx, εpy, and εpz are the plastic strain on the x, y and z axis, respectively; and γpxy, γpyz, and γpxz are the shear stress component.Substituting equation (10) and equation (5) into equation (9).(12)λ=∂F/∂σijDijkledεklA,where(13)A=∂F∂σijDijkle−∂F∂κ∂κ∂ε¯p+∂κ∂T∂T∂ε¯pdε¯pdε¯ijP‐∂F∂εklp∂F∂σkl.Substituting equations (12) and (13) into equation (3) and then substituting equation (5) into equation (8), the elastic-plastic matrix can be obtained.(14)dσij=Dijklepdεkl,(15)Dijklep=Dijkle‐Dijkle∂F/∂σpq∂F/∂σrsDpqrseA.


## 2.3. Plastic Yield Criteria for Rocks

Rock deformation obeys the Mohr-Coulomb yield criterion, and we stipulate that the tensile stress is positive and the compressive stress is negative [22]:(16)τ=c−σtanϕ.According to the geometry of the mole circle, shear stressτ and normal stress σ can be expressed as follows [34]:(17)τ=σ1−σ32cosϕ,σ=σ1+σ32+σ1−σ32sinϕ,where τ is the shear stress, σ is the normal stress, φ is the friction degree, σ1 and σ3 are the maximum and minimum principal stresses, respectively, and σ2 is the intermediate principal stress.Substitute equation (17) into equation (16), we have(18)σ1−σ32cosϕ=c−σ1+σ32+σ1−σ32sinϕtanϕ.Principal stress and principal deviatoric stress satisfy the following equation [32]:(19)σ1=23J2cosΘ+13I1,σ2=23J2cos2π3‐Θ+13I1,σ3=23J2cos2π3+Θ+13I1,where σ1 and σ3 are the maximum and minimum principal stresses, respectively, σ2 is the intermediate principal stress, I1 is the first invariant of stress tensor, J2 is the second invariant of deviatoric stress, and Θ is the deviatoric polar angle.Substitute equation (19) into equation (18):(20)F=13I1sinϕ+J2sinΘ+π3+J23cosΘ+π3sinϕ−ccosϕ=0.Equivalent compressive stressp, Mises equivalent stress q, the first invariant of stress tensor I1, and the second invariant of deviatoric stress J2 satisfy the following equation [25]:(21)p=‐13traceσ==‐I13,q=32sijsij=3J2,cos3Θ=r3q3,where p is the equivalent pressure stress, q is the Mises equivalent stress, r is the third invariant of deviatoric stress, and s is the deviatoric stress.Substitute equation (21) into equation (20):(22)F=fp,q,Θ=Rmcq−ptanϕ−c=0,where(23)Rmc=13cosϕsinΘ+π3+13cosΘ+π3tanϕ.In the case of equation (23), there will be sharp angles on the yield surface, leading to not unique plastic flow direction, which leads to tedious calculation and slow convergence. In order to solve this problem, Menetrey and Willam [35] proposed a continuous smooth elliptic function as a plastic potential surface.(24)G=εc0tanψ2+Rmwq2−ptanψ,where(25)Rmw=41−e2cos2Θ+2e−1221−e2cosΘ+2e−141−e2cos2Θ+5e2−4e3−sinϕ6cosϕ,(26)e=3−sinϕ3+sinϕ,where ψ is the dilation angle measured in the p−Rmwq plane at high confining pressure and can depend on temperature and predefined field variables; c0 is the initial cohesion yield stress; ε is the eccentricity ratio, and here is assumed to be 0.1; e is the “out of roundedness” of the deviatoric section in terms of the ratio between the shear stress along the extension meridian and the shear stress along the compression meridional.


## 2.4. Fluid Flow in Fracture and Formation

Once the fracture is formed, fluid will flow into the fracture at once. Fluid flow in the fracture follows the mass conservation law [36]:(27)∂w∂t+∂qt∂s+vt+vb=0,where qt is the injection flow rate, w is the fracture width, s is the fracture length, t is injection time, and vt and vb are the fracturing fluid leakoff rate through the top and bottom fracture surfaces, respectively.Both tangential flow along and normal flow across the fracture surfaces can be calculated using the Poiseuille equation [37, 38]:(28)qt=w312μ∂p∂s,where μ is the fluid viscosity and p is the fluid pressure.Fluid leakoff rate is mainly affected by pressure and leakoff coefficient, which satisfy the following equation [39, 40]:(29)vt=ctp−ptop,vb=ctp−pbot,where ct is the fluid leakoff coefficients and ptop and pbot are the pressure in the top and bottom fracture surface, respectively.When the fluid leaks off from the fracture, it will flow through the formation in accordance with the Darcy law. The mass conservation equation is expressed as [34]:(30)∇ρwkwbwμw∇pw+qw=∂∂tφρwswbw,where ρw is the water permeability, kw is the water permeability, pw is water pressure in formation, bw is the water compression coefficient, μw is the water viscosity, qw is the leakoff fluid, Sw is the water saturation, and ϕ is the effective porosity.


## 2.5. Cohesive Crack Model

The singularity of the fracture tip stress field will affect the convergence of the calculation results. In order to solve this problem, CZM is adopted to simulate the initiation and propatation of fractures. This method is implemented by embedding an artificially predefined path of cohesive element into the formation. The cohesive element glues two conjugated cohesive surfaces of the reservoir pressure/displacement element. The failure of these element obeys the traction/separation law in the vertical or shear direction [41, 42]. When the cohesive element is completely destroyed, the two bound cohesive surfaces are separated, and hydraulic fracture is generated. When the maximum nominal stress criterion is adopted for rock failure, initiation propagation occurs as the maximum nominal stress ratio reaches a critical value [43–45]:(31)fc=maxTnTn0,TsTs0,(32)1≤fc≤1+ft,where δn and δs are, respectively, the normal and tangential displacement jump vector, δeq are effective displacement, Tn0 and Ts0 are, respectively, the normal and tangential stress component of the stress element under a linear elastic condition, < > denotes the Macaulay bracket, fc is a parameter evaluating the fracture failure, and ft is assured to be 0.05.The fracture damage evolution is determined by damage variables, and the expressions of damage effects on the normal and tangential stress components of the element can be obtained as [44–46](33)Tn=1−dTn0,Tn0≥0Tn0,Tn0<0,(34)Ts=1−dTs0,where d is the scalar damage variable between fractures, 0<d<1, and Tn and Ts are, respectively, the values of the normal stress component and the tangential stress component of the traction/separation model.When the critical fracture energy of rock material along the first and second shear directions is similar, the damage evolution during fracture propagation can also be determined by the Benzeggagh–Kenane fracture criterion [47]. The combined energy dissipated by failure Gc is defined as(35)Gc=GIc+GIIc−GIcGIIc−GIcη,where Gc is the critical fracture energy release rate, GcI and GcII are the normal and first shear direction fracture toughness values, GI and GII are the normal and shear direction fracture energy release rate, and η is a material properties parameter; here is assumed to be 2.3.


## 2.6. Model Solution

The iterative coupling approach is used to solve the coupled model. Fluid flow and solid deformation are solved separately and sequentially, and the coupling terms are iterated until convergence. The fluid flow model is solved first based on the initial guess of injection pressure; then, the solid deformation is solved by the XFEM stiffness matrix, which is based on the solution of the previous iteration. The pressure distribution is calculated until convergence is reached. If the convergence criterion is not satisfied, the guessing value is modified aspn+1=ωpn−1+1−ωpn; the pressure distribution and rock deformation will be recalculated with the modified value. When the rock failure condition is met, the fracture will propagate and add a new element for the next time step.


## 3. Model Verification

In order to simulate HF in elastic-plastic formation, the Abaqus 6.14 software has been utilized to solve this coupling model. The size of geological model is set as50×50 m with a cohesive element embedded in the central of the model. The injection point is in the center of the cohesive zone. The minimum principal stress is in y direction. All the external boundary displacements are fixed in normal directions. The geological model is shown in Figure 1, and the input parameters of simulation model are presented in Table 1. The simulation process consists of two steps: geostatic step and HF step. The purpose of geostatic step is to balance the initial formation stress and fluid pressure. HF step is to simulate the change of fluid injection, rock deformation, and fracture initiation and expansion during fracturing.Figure 1
Geological model for hydraulic fracture simulation.Table 1
Input parameters for hydraulic fracture model.
PropertiesValueYoung modulus (MPa)15000Poisson ratio0.25Material cohesion (MPa)1.5Material friction angle (°)28Material dilation angle (°)28Formation permeability coefficients (m/s)1×10−7Maximum horizontal stress (MPa)20Minimum horizontal stress (MPa)15Porosity0.1Fluid viscosity (cp)1Specific weight of fluid (KN/m3)9.8Initial pore pressure (MPa)0.1Injection rate (m3/s)5×10−3Pressure dependent leakoff coefficient (m3/MPa.s)1×10−14Critical fracture energy (KN/m3)28Damage initiation stress (MPa)0.32To verify the accuracy of the model, we compared the calculated results with KGD model. KGD model is a simple analytical model. The wellbore fracture width (wo) and injection pressure solution (pw) of KGD model can be given [48]:(36)wo=2.36q3μ1−v2Eh31/6t1/3,(37)pw=1.09μE21−v221/3t−1/3+σmin.We used the calculated dimensionless injection pressure (ratio of injection pressure to maximum injection pressure) and dimensionless wellbore fracture width (ratio of wellbore fracture width to maximum wellbore fracture width) to valid the model; the simulated results are presented in Figures2 and 3. After the model is validated, the changes of injection pressure, fracture width, fracture length, and the equivalent plastic strain zone with time in elastic and elastic-plastic formation are also investigated by the model. The calculated results are demonstrated in Figures 4–6.Figure 2
Comparison of calculated dimensionless injection pressure with analytical model.Figure 3
Comparison of calculated dimensionless wellbore fracture width with analytical model.Figure 4
Comparison of injection pressure and wellbore fracture width calculated by poroelastic and poroelastoplastic models.Figure 5
Comparison of fracture parameters with different time calculated in poroelastic and poroelastoplastic formation.Figure 6
Equivalent plastic strain zone of fracture varies with injection time.
(a)
Injection time at 20 s(b)
Injection time at 40 s(c)
Injection time at 60 s(d)
Injection time at 80 sFigures2 and 3 exhibit that the simulated dimensionless injection pressure and dimensionless wellbore fracture width are in good agreement with KGD model. The change trend of injection pressure and fracture width at the wellbore is also in accordance with KGD model. Injection pressure in both models is first increased and then decreased to a constant value. Wellbore fracture width is also gradually increasing with injection time. In early injection time, there is a slight deviation in between dimensionless wellbore fracture width calculated by our model and that calculated by KGD model. This may be caused by fluid leakoff. On the whole, the agreement between the two models is encouraging.Figure4 shows that the injection pressure and fracture width at the wellbore in elastic-plastic formation are greater than that in elastic formation, and the breakdown pressure and propagation pressure are also greater in elastic-plastic formation. Injection pressure is first increased and then decreased to a constant value. Wellbore fracture width is also gradually increasing with injection time.Figure5 presents that fracture in elastic-plastic formation is incline to be wider and shorter than elastic formation. Because plastic damage requires more energy, plastic deformation of rock leads the fracture tip to be blunt, resulting fracture more difficult to extend. Figure 6 displays that the equivalent plastic strain area increases with injection time, and the largest plastic strain appears near the fracture tip. The plastic deformation also causes rock compaction, which results that the porosity and permeability decrease around the fracture with less fluid leakoff and higher injection pressure.


## 4. Numerical Results and Discussion

### 4.1. Effect of Stress Difference on Fracture Propagation

The effects of horizontal principal stress differences on fracture width, length, injection pressure, fracture width at the wellbore, equivalent plastic strain of fracture, pore pressure, and induced stress are simulated. It is assumed that the minimum horizontal principal stress is constant at 30 MPa; the horizontal principal stress difference is 5 MPa, 8 MPa, and 10 MPa, respectively; the pore pressure of the formation is 20 MPa; the injection time is 80s; and the friction angle and dilation angle of the rock are both 15°. The calculation results are presented in Figures7–11.Figure 7
Effect of stress difference on fracture width and length.Figure 8
Effect of stress difference on wellbore fracture width and injection pressure.Figure 9
Effect of stress difference on equivalent plastic strain zone of fracture.
(a)
Stress difference 5 MPa(b)
Stress difference 8 MPa(c)
Stress difference 10 MPaFigure 10
Effect of stress difference on pore pressure.
(a)
Stress difference 5 MPa(b)
Stress difference 8 MPa(c)
Stress difference 10 MPaFigure 11
Effect of stress difference on stress onx direction.
(a)
Stress difference 5 MPa(b)
Stress difference 8 MPa(c)
Stress difference 10 MPaFigure7 demonstrates that crack width increases and crack length decreases with the increasing stress difference. Figure 8 indicates that the crack width at the wellbore and injection pressure increases with the increase of horizontal stress difference. Figure 9 displays that the area of equivalent plastic strain around the crack also increases with the increase of horizontal stress difference. Figures 10 and 11 present that pore pressure and induced stress in x direction surround the fracture also grows with the increasing stress difference. Because high stress difference is easy to generate large plastic zones, and rock compaction caused by the plastic deformation around the fracture leads to a decrease in porosity and permeability, which in turn makes it difficult for the fracturing fluid to flow into the formation, resulting greater pore pressure near the fracture surface. Meanwhile, the induced stress on the fracture surface also makes crack difficult to extent, because of the plastic effect. Therefore, plastic formations with high horizontal stress difference are benefit to form wide and short crack.


### 4.2. Effect of Friction Angle on Fracture Propagation

In order to study the influence of friction angle on hydraulic fracture propagation, the effects of friction angle on fracture width and length, injection pressure, fracture width at the wellbore, equivalent plastic strain, pore pressure, and induced stress are simulated. It is assumed that the minimum horizontal principal stress is constant at 10 MPa, the horizontal principal stress difference is 5 MPa, and the injection time is 80 s. The friction angle of the rock is 10°, 20°, and 28°, respectively, and dilation angle is equal to friction angle. The calculation results are presented in Figures12–16.Figure 12
Effect of stress difference on fracture width and length.Figure 13
Effect of stress difference on wellbore fracture width and injection pressure.Figure 14
Effect of friction angle on equivalent plastic strain zone of fracture.
(a)
Friction angle 10°(b)
Friction angle 20°(c)
Friction angle 28°Figure 15
Effect of friction angle on pore pressure.
(a)
Friction angle 10°(b)
Friction angle 20°(c)
Friction angle 28°Figure 16
Effect of friction angle on stress onx direction.
(a)
Friction angle 10°(b)
Friction angle 20°(c)
Friction angle 28°Figure12 demonstrates that crack width decreases with the increasing friction angle, while crack length increases with the increasing friction angle. Lower friction angle causes the fracture tip to be blunter. When friction angle is 10°, plastic deformation at the injection point has obviously appeared in the early stage of crack expansion, leading to a significant increase of crack width near the wellbore. Figure 13 indicates that the crack width at the wellbore and injection pressure increases with friction angle decrease. When the friction angle decreases to 10°, fracture breakdown pressure and propagation pressure grow significantly, which leads to rapid increasing fracture width near the injection point. Figure 14 displays that the equivalent plastic strain area around the crack increases with the decreasing friction angle. Plastic strain occurs obviously at the fracture entrance with friction angle of 10°, which is difficult for crack to extent. Figures 15 and 16 present that pore pressure and stress in x direction surround the fracture also increases with the decreasing friction angle.Rocks with low friction angle are easy to occur plastic deformation and form wide and short crack.


### 4.3. Effect of Dilation Angle on Fracture Propagation

The dilation angle reflects the change rate of rock volume, normal displacement, and tangential displacement in the shearing process. To study the influence of dilation angle on hydraulic fracture propagation, the effects of dilation angle on fracture width and length, injection pressure, fracture width at the wellbore, equivalent plastic strain, pore pressure, and induced stress are simulated. It is assumed that the dilation angle of the rock is 10°, 20°, and 28°, respectively, and friction angle is 28°. The other input parameters are the same as Section4.2. The calculation results are presented in Figures 17–21.Figure 17
Effect of dilation angle on fracture width and length.Figure 18
Effect of dilation angle on wellbore fracture width and injection pressure.Figure 19
Effect of dilation angle on equivalent plastic strain zone of fracture.
(a)
Dilation angle 10°(b)
Dilation angle 20°(c)
Dilation angle 28°Figure 20
Effect of dilation angle on pore pressure.
(a)
Dilation angle 10°(b)
Dilation angle 20°(c)
Dilation angle 28°Figure 21
Effect of dilation angle on stress onx direction.
(a)
Dilation angle 10°(b)
Dilation angle 20°(c)
Dilation angle 28°Figure17 demonstrates that crack width decreases with the increasing dilation angle, and this trend is not obvious. Figure 18 indicates that with injection time more than 20 s, dilation angle has little influence on injection pressure and fracture width at the wellbore. Figure 19 displays that the equivalent plastic strain around the crack increases with the decrease of dilation angle, and the value of plastic strain is small. Figures 20 and 21 present that dilation angle has also little effect on pore pressure and stress on x direction. Compared with friction angle, dilation angle is less sensitive to fracture width and length, as well as injection pressure.


## 4.1. Effect of Stress Difference on Fracture Propagation

The effects of horizontal principal stress differences on fracture width, length, injection pressure, fracture width at the wellbore, equivalent plastic strain of fracture, pore pressure, and induced stress are simulated. It is assumed that the minimum horizontal principal stress is constant at 30 MPa; the horizontal principal stress difference is 5 MPa, 8 MPa, and 10 MPa, respectively; the pore pressure of the formation is 20 MPa; the injection time is 80s; and the friction angle and dilation angle of the rock are both 15°. The calculation results are presented in Figures7–11.Figure 7
Effect of stress difference on fracture width and length.Figure 8
Effect of stress difference on wellbore fracture width and injection pressure.Figure 9
Effect of stress difference on equivalent plastic strain zone of fracture.
(a)
Stress difference 5 MPa(b)
Stress difference 8 MPa(c)
Stress difference 10 MPaFigure 10
Effect of stress difference on pore pressure.
(a)
Stress difference 5 MPa(b)
Stress difference 8 MPa(c)
Stress difference 10 MPaFigure 11
Effect of stress difference on stress onx direction.
(a)
Stress difference 5 MPa(b)
Stress difference 8 MPa(c)
Stress difference 10 MPaFigure7 demonstrates that crack width increases and crack length decreases with the increasing stress difference. Figure 8 indicates that the crack width at the wellbore and injection pressure increases with the increase of horizontal stress difference. Figure 9 displays that the area of equivalent plastic strain around the crack also increases with the increase of horizontal stress difference. Figures 10 and 11 present that pore pressure and induced stress in x direction surround the fracture also grows with the increasing stress difference. Because high stress difference is easy to generate large plastic zones, and rock compaction caused by the plastic deformation around the fracture leads to a decrease in porosity and permeability, which in turn makes it difficult for the fracturing fluid to flow into the formation, resulting greater pore pressure near the fracture surface. Meanwhile, the induced stress on the fracture surface also makes crack difficult to extent, because of the plastic effect. Therefore, plastic formations with high horizontal stress difference are benefit to form wide and short crack.


## 4.2. Effect of Friction Angle on Fracture Propagation

In order to study the influence of friction angle on hydraulic fracture propagation, the effects of friction angle on fracture width and length, injection pressure, fracture width at the wellbore, equivalent plastic strain, pore pressure, and induced stress are simulated. It is assumed that the minimum horizontal principal stress is constant at 10 MPa, the horizontal principal stress difference is 5 MPa, and the injection time is 80 s. The friction angle of the rock is 10°, 20°, and 28°, respectively, and dilation angle is equal to friction angle. The calculation results are presented in Figures12–16.Figure 12
Effect of stress difference on fracture width and length.Figure 13
Effect of stress difference on wellbore fracture width and injection pressure.Figure 14
Effect of friction angle on equivalent plastic strain zone of fracture.
(a)
Friction angle 10°(b)
Friction angle 20°(c)
Friction angle 28°Figure 15
Effect of friction angle on pore pressure.
(a)
Friction angle 10°(b)
Friction angle 20°(c)
Friction angle 28°Figure 16
Effect of friction angle on stress onx direction.
(a)
Friction angle 10°(b)
Friction angle 20°(c)
Friction angle 28°Figure12 demonstrates that crack width decreases with the increasing friction angle, while crack length increases with the increasing friction angle. Lower friction angle causes the fracture tip to be blunter. When friction angle is 10°, plastic deformation at the injection point has obviously appeared in the early stage of crack expansion, leading to a significant increase of crack width near the wellbore. Figure 13 indicates that the crack width at the wellbore and injection pressure increases with friction angle decrease. When the friction angle decreases to 10°, fracture breakdown pressure and propagation pressure grow significantly, which leads to rapid increasing fracture width near the injection point. Figure 14 displays that the equivalent plastic strain area around the crack increases with the decreasing friction angle. Plastic strain occurs obviously at the fracture entrance with friction angle of 10°, which is difficult for crack to extent. Figures 15 and 16 present that pore pressure and stress in x direction surround the fracture also increases with the decreasing friction angle.Rocks with low friction angle are easy to occur plastic deformation and form wide and short crack.


## 4.3. Effect of Dilation Angle on Fracture Propagation

The dilation angle reflects the change rate of rock volume, normal displacement, and tangential displacement in the shearing process. To study the influence of dilation angle on hydraulic fracture propagation, the effects of dilation angle on fracture width and length, injection pressure, fracture width at the wellbore, equivalent plastic strain, pore pressure, and induced stress are simulated. It is assumed that the dilation angle of the rock is 10°, 20°, and 28°, respectively, and friction angle is 28°. The other input parameters are the same as Section4.2. The calculation results are presented in Figures 17–21.Figure 17
Effect of dilation angle on fracture width and length.Figure 18
Effect of dilation angle on wellbore fracture width and injection pressure.Figure 19
Effect of dilation angle on equivalent plastic strain zone of fracture.
(a)
Dilation angle 10°(b)
Dilation angle 20°(c)
Dilation angle 28°Figure 20
Effect of dilation angle on pore pressure.
(a)
Dilation angle 10°(b)
Dilation angle 20°(c)
Dilation angle 28°Figure 21
Effect of dilation angle on stress onx direction.
(a)
Dilation angle 10°(b)
Dilation angle 20°(c)
Dilation angle 28°Figure17 demonstrates that crack width decreases with the increasing dilation angle, and this trend is not obvious. Figure 18 indicates that with injection time more than 20 s, dilation angle has little influence on injection pressure and fracture width at the wellbore. Figure 19 displays that the equivalent plastic strain around the crack increases with the decrease of dilation angle, and the value of plastic strain is small. Figures 20 and 21 present that dilation angle has also little effect on pore pressure and stress on x direction. Compared with friction angle, dilation angle is less sensitive to fracture width and length, as well as injection pressure.


## 5. Conclusions

A fluid-solid coupling elastoplasticity fracture propagation model has been established. Effects of plasticity and damage with have been accounted. The extended finite element methods combined with cohesive zone are utilized to approach to model the fracture process in permeable rocks. The Abaqus 6.14 software has been used to investigate fracture propagation in elastoplastic formation. According to the above analysis, it can be concluded as follows:(1)
Calculation results by the model are compared with KGD model, and they have good agreements. It is suitable to simulate fracture propagation in fracturing treatment(2)
Compared with elastic formation, fracture propagation in elastoplastic formation is more difficult, the breakdown pressure and extending pressure are greater, and fracture shape is wider and shorter. Plastic deformation also leads to the fracture tip to be blunt. Plastic strain area grows with the increasing injection time, and the largest plastic strain appears near the fracture tip(3)
Fracture width increases, and fracture length decreases with the increasing stress difference. Under the condition of high stress difference or low friction angle, it is easy to occur large plastic deformation area, result in great breakdown pressure, propagation pressure, induced stress inx direction, and pore pressure. Soft formations with great stress difference or low friction angle are inclined to form wide and short crack with blunt fracture tip(4)
Compared with friction angle, dilation angle is less sensitive to plastic deformation and fracture parameters. For the formation with high stress difference and friction angle, the influence of plasticity deformation on fracture propagation should not be ignored



---

*Source: 6680023-2021-08-26.xml*