Application of Extended Finite Element Method in Fracture Propagation Simulation of Plastic Formation


 Current hydraulic fracture models are mainly based on elastic theories, which fail to give accurate prediction of fracture parameters in plasticity formation. This paper proposed a fluid–solid coupling model for fracture propagation in elastoplastic formations. The rock plastic deformation in the model satisfied the Mohr-Coulomb yield criterion and plastic strain increment theory. The coupled model is solved by using extended finite-element method(XFEM) and the cohesive zone method (CZM). The accuracy of the model is verified by comparing the calculated results with existing models. The influences of stress difference, friction angle and dilation angle on fracture shape (length, width), injection pressure, plastic deformation, induced stress and pore pressure are investigated. The results indicate that compared with elastic formation, fracture shape in elastoplastic formation is wider and shorter and fracture propagation is more difficult with greater breakdown pressure and extending pressure. Plastic deformation also cause blunt fracture tip. High stress difference or low friction angle formations tend to occur large plastic deformation area and form wide and short fracture. Compared with friction angle, dilation angle is less sensitive to plastic deformation and fracture parameters and geometry. For the formation with high stress difference and friction angle, plasticity effects on fracture propagation should not be ignored.


Introduction 24
With the increasing energy consumption, efficient development of oil and gas resources has 25 become the focus of attention. Hydraulic fracturing (HF), as a widely used technology, has become 26 increasingly important in oil and gas stimulation (Hu et  propagation mechanism has a great significance for fracture parameter (with, length, height) 28 optimization. The classical fracture propagation models treated formation as an elastic medium 29 without considering the plasticity effects and fluid-solid coupling. (Perkins and Kern, 1961; 30 Geertsma and Klerk, 1969;Nordgren, 1972;Saouma et al, 1987;Dean and Schmidt, 2009;Dahi-31 Taleghani and Olson JE, 2011; Leonhart and Meschke, 2011). However, plastic failure usually 32 occurs in soft rocks of high temperature and pressure formation and unconsolidated formation 33 during fracturing (Brünig, 1999; Abu Al-Rub and Voyiadjis, 2003; Bohloli and de Pater, 2006). Previous studies on the effect of plasticity in HF were just based on the static fracture and 44 constant injection pressure in fracture without considering the coupling between the pressure 45 diffusion and rock deformation. Papanastasiou (Papanastasiou, 1997(Papanastasiou, , 1999(Papanastasiou, , 2012  Where dσij is the stress increment, dεij is the strain increment, dui,j is the medium displacement 87 increment, dfi is the body force increment, D ep ijkl is elastic-plastic coefficient matrix. 88 The displacement and stress at the outer boundary and fluid pressure at the fracture surface are 89 where Γu, Γt, and Γcr are the displacement boundary condition, stress boundary condition, and 91 fracture pressure boundary condition, respectively; u* and t* are the corresponding displacement 92 and stress increment on the boundary; and dp is fluid pressure increment on the fracture. 93 94

The elastic-plastic increment constitutive model 95
When the rock begins to occur plastic failure, the strain increment of the rock is composed of 96 elastic strain increment and plastic strain: 97 The elastic strain increment dε e ij and the stress increment dσij satisfy the Hooke's law: 99 The plastic strain increment follows the flow rule, which is defined as (Vermeer and Borst, 101 In the process of plastic deformation, the yield surface of the rock is a function of strengthening 104 parameters, plastic strain and stress: 105 Where εij is the total strain; ε e Ij is the elastic strain; ε p Ij is the plastic strain; σij is the stress tensor; 108 D e ijkl is the elastic constitutive tensor; F is the load surface equation;  is the strengthening parameter 109 of rock; λ is a quantity related to the slope of the uniaxial curve of stress σ -plastic strain ε p , which 110 can be determined by yield criteria and strengthening theory. 111 The stress increment is obtained by substituting equation (4) into (3). 112 Multiplying both sides of equation (8) by ∂ / ∂ and substituting it into equation (7): The hardening parameter is a function of the equivalent plastic strain and temperature, which 116 is determined as: 117 The equivalent plastic strain increment is defined as: 119 Where ε P is the equivalent plastic strain; T is temperature; ε p x 、ε p y 、ε p z are the plastic strain on 121 the x, y and z axis respectively; γ p xy 、γ p yz 、γ p xz is the shear stress component.

Plastic yield criteria for rocks 132
Rock deformation obeys the Mohr-Coulomb yield criterion, and we stipulates that the tensile 133 stress is positive and the compressive stress is negative: 134 According to the geometry of the mole circle, shear stress τ and normal stress σ can be 136 expressed as following: Where τ is the shear stress; σ is the normal stress; φ is the friction degree; σ1 and σ3 are the 139 maximum and minimum principal stresses, respectively; σ2 are the Intermediate principal stress; 140 Substitute equation (17) into equation (16), we have: 141 Principal stress and principal deviatoric stress satisfy the following equation: 143 Equivalent compressive stress p, Mises equivalent stress q, the first invariant of stress tensor 150 I1 and the second invariant of deviatoric stress J2 satisfy the following equation: 151 Where p is the equivalent pressure stress; q is the Mises equivalent stress; r is the third invariant 153 of deviatoric stress; s is the deviatoric stress.
In the case of equation (23), there will be sharp angles on the yield surface, leading to not 160 unique plastic flow direction, which leads to tedious calculation and slow convergence. In order to 161 solve this problem, Menetrey and Willam (Menetrey and Willam, 1995) proposed a continuous 162 smooth elliptic function as a plastic potential surface. 163 Where 165   cos 6 sin 3 Where  is the dilation angle measured in the p-Rmwq plane at high confining pressure and can 168 depend on temperature and predefined field variables; c0 is the initial cohesion yield stress;  is the 169 eccentricity ratio, here is assumed to be 0.1; e is the "out of roundedness" of the deviatoric section 170 in terms of the ratio between the shear stress along the extension meridian and the shear stress along 171 the compression meridional. 172 173

Fluid flow in fracture and formation 174
Once the fracture is formed, fluid will flow into the fracture at once, fluid flow in the fracture 175 follows the mass conservation law (Yao, 2012 Where qt is the injection flow rate, w is the fracture width, s is the fracture length, t is injection Where ct is the fluid leakoff coefficients, ptop and pbot are pressure in the top and bottom 188 fracture surface, respectively. 189 When the fluid leaks off from the fracture, it will flows through the formation in accordance 190 with the Darcy law. The mass conservation equation is expressed as: 191 bw is the water compression coefficient; μw is the water viscosity; qw are the leakoff fluid; Sw is the 194 water saturation;  is the effective porosity.

Cohesive crack model 197
The singularity of the fracture tip stress field will affect the convergence of the calculation of 198 the model. In order to solve this problem, CZM is adopted to simulate the initiation and propatation 199 of fractures. This method is implemented by embedding an artificially predefined path of cohesive 200 element into the formation. The cohesive element glues two conjugated cohesive surfaces of the 201 reservoir pressure/displacement element. The failure of these element obeys the traction/separation 202 law in the vertical or shear direction (Tvergaard and Hutchinson, 1992;Camacho and Ortiz, 1996). 203 When the cohesive element is completely destroyed, the two bound cohesive surfaces are separated 204 and hydraulic fracture is generated. When the maximum nominal stress criterion is adopted for rock Benzeggagh-Kenane fracture criterion (Benzeggagh and Kenane, 1996). The combined energy 224 dissipated by failure Gc is defined as 225 where Gc is the critical fracture energy release rate, G c Ⅰ and G c Ⅱ are the normal and first shear 227 direction fracture toughness values, GⅠ and GⅡ are the normal and shear direction fracture energy release rate, and η is a material properties parameter, here is assumed to be 2.3. 229 230

3、Model verification 231
In order to simulate HF in elastic-plastic formation, the size of geological model is set as 50 m 232

Effect of stress difference on fracture propagation 282
The effects of horizontal principal stress differences on fracture width, length, injection 283 pressure, fracture width at the wellbore, equivalent plastic strain of fracture, pore pressure and 284 induced stress are simulated. It is assumed that the minimum horizontal principal stress is constant 285 at 30MPa, the horizontal principal stress difference is 5MPa, 8MPa and 10MPa respectively, the 286 pore pressure of the formation is 20MPa, the injection time is 80s, and the friction angle and dilation 287

Poroelastplasticity-20s
Poroelastplasticity-40s Poroelastplasticity-60s Poroelasticity-20s Poroelasticity-40s Poroelasticity-60s angle of the rock are both 15°. The calculation results are presented in figure 6-10. 288 Fig. 6 demonstrates that crack width increases and crack length decreases with the increasing 289 stress difference. Fig. 7 indicates that the crack width at the wellbore and injection pressure increases 290 with the increase of horizontal stress difference. Figure 8 displays that the area of equivalent plastic 291 strain around the crack also increases with the increase of horizontal stress difference. Figure 9-10 292 present that pore pressure and induced stress in x direction surround the fracture also grows with 293 the increasing stress difference. Because high stress difference is easy to generate large plastic zones, 294 and rock compaction caused by the plastic deformation around the fracture leads to a decrease in 295 porosity and permeability, which in turn makes it difficult for the fracturing fluid to flow into the 296 formation, resulting greater pore pressure near the fracture surface. Meanwhile, the induced stress 297 on the fracture surface also makes crack difficult to extent, owing to plastic effect. Therefore, plastic 298 formations with high horizontal stress difference is benefit to form wide and short crack.

Effect of friction angle on fracture propagation 322
In order to study the influence of friction angle on hydraulic fracture propagation, the effects 323 of friction angle on fracture width and length, injection pressure, fracture width at the wellbore, 324 equivalent plastic strain, pore pressure and induced stress are simulated. It is assumed that the 325 minimum horizontal principal stress is constant at 10MPa, the horizontal principal stress difference 326 is 5MPa, and the injection time is 80s. The friction angle of the rock is 10°, 20° and 28°, respectively, 327 and dilation angle is equal to friction angle. The calculation results are presented in figure 11-15. 328 in the early stage of crack expansion, leading to a significant increase of crack width near the 332 wellbore. Figure 12 indicates that the crack width at the wellbore and injection pressure increases 333 with friction angle decrease. When the friction angle decreases to 10°, fracture breakdown pressure 334 and propagation pressure grow significantly, which leads to rapid increasing fracture width near the 335 injection point. Figure 13 displays that the equivalent plastic strain area around the crack increases 336 with the decreasing friction angle. Plastic strain occur obviously at the fracture entrance with friction 337 angle of 10°, which is difficult for crack to extent. Figure 14-15 present that pore pressure and stress 338 in x direction surround the fracture also increases with the decreasing friction angle. 339 Rock with low friction angle are easy to occur plastic deformation and form wide and short 340 crack. 341 342

Effect of dilation angle on fracture propagation 361
The dilation angle reflects the change rate of rock volume, normal displacement and tangential 362 displacement in the shearing process. To study the influence of dilation angle on hydraulic fracture 363 propagation, the effects of dilation angle on fracture width and length, injection pressure, fracture 364 width at the wellbore, equivalent plastic strain, pore pressure and induced stress are simulated. It is 365 assumed that the dilation angle of the rock is 10°, 20° and 28°, respectively, and friction angle is  trend is not obvious. Figure 17 indicates that with injection time more than 20s, dilation angle has 370 little influence on injection pressure and fracture width at the wellbore. Figure 18 displays that the 371 equivalent plastic strain around the crack increases with the decrease of dilation angle and the value 372 of plastic strain is small. Figure 19-20 present that dilation angle has also little effect on pore 373 pressure and stress on x direction. Compared with friction angle, dilation angle is less sensitive to 374 fracture width and length, as well as injection pressure. 375 In comparison to elastic formation, elastic-plastic formation exhibit greater breakdown and 397 propagation pressure, resulting wider and shorter fracture. Plastic deformation also lead to the 398 fracture tip to be blunt. Plastic strain area grows with increasing injection time, and the largest plastic 399 strain appears near the fracture tip. 400 Fracture width increases and fracture length decreases with the increasing stress difference. 401 High stress difference or low friction angle are easy to occur large plastic deformation area, resulting 402 great breakdown pressure, propagation pressure, induced stress in x direction and pore pressure. 403 Soft formations with great stress difference or low friction angle is inclined to form wide and short 404 crack and blunt fracture tip. 405 Compared with friction angle, dilation angle is less sensitive to plastic deformation and fracture 406 parameters. For the formation with high stress difference and friction angle, the influence of 407 plasticity effects on fracture propagation should not be ignored. 408 409