Numerical Simulation for Fracture Propagation in Elastoplastic Formations

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 ﬂ uid – solid coupling model for fracture propagation in elastoplastic formations. The rock plastic deformation in the model satis ﬁ es the Mohr-Coulomb yield criterion and plastic strain increment theory. The extended ﬁ nite-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 e ﬀ ects of stress di ﬀ erence, 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 di ﬃ cult, 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 di ﬀ erence 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 di ﬀ erence and friction angle, the e ﬀ ect of plasticity deformation on fracture propagation should not be ignored.


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][4][5][6][7][8]. However, plastic failure usually occurs in soft layer or unconsolidated formation with high temperature and pressure during fracturing [9][10][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][13][14][15][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][18][19][20][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][23][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.

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 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]: where dσ ij is the stress increment, dε ij is the strain increment, du i,j is the medium displacement increment, df i is the body force increment, and D ep ijkl is the elastic-plastic coefficient matrix. The displacement and stress at the outer boundary and fluid pressure at the fracture surface are [31] 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.

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]: The elastic strain increment dεe ij and the stress increment dσ ij satisfy the Hooke's law: The plastic strain increment follows the flow rule, which is defined as [33]: In the process of plastic deformation, the yield surface of the rock is a function of strengthening parameters, plastic strain, and stress [32]: where ε ij is the total strain, εe Ij is the elastic strain, εp Ij is the plastic strain, σ ij is the stress tensor, De ijkl 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).
Multiplying both sides of equation (8) by ∂F/∂σ ij and substituting it into equation (7): The hardening parameter is a function of the equivalent plastic strain and temperature, which is determined as [32]: The equivalent plastic strain increment is defined as follows: where − ε P is the equivalent plastic strain; T is the temperature; εp x, εp y, and εp z are the plastic strain on the x, y and z axis, respectively; and γp xy, γp yz, and γp xz are the shear stress component. Substituting equation (10) and equation (5) into equation (9). where Substituting equations (12) and (13) into equation (3) and then substituting equation (5) into equation (8), the elastic-plastic matrix can be obtained.   3 Geofluids the tensile stress is positive and the compressive stress is negative [22]: According to the geometry of the mole circle, shear stress τ and normal stress σ can be expressed as follows [34]: 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 Principal stress and principal deviatoric stress satisfy the following equation [32]: where σ 1 and σ 3 are the maximum and minimum principal stresses, respectively, σ 2 is the intermediate principal stress, I 1 is the first invariant of stress tensor, J 2 is the second invariant of deviatoric stress, and Θ is the deviatoric polar angle. Substitute equation (19) into equation (18): Equivalent compressive stress p, Mises equivalent stress q, the first invariant of stress tensor I 1 , and the : 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): where 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] where where ψ is the dilation angle measured in the p − R mw q plane at high confining pressure and can depend on temperature and predefined field variables; c 0 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.

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]: where q t is the injection flow rate, w is the fracture width, s is the fracture length, t is injection time, and v t and v b are the fracturing fluid leakoff rate through the top and bottom fracture surfaces, respectively.

Geofluids
Both tangential flow along and normal flow across the fracture surfaces can be calculated using the Poiseuille equation [37,38]: 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]: where c t is the fluid leakoff coefficients and p top and p bot 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]: where ρ w is the water permeability, k w is the water permeability, p w is water pressure in formation, b w is the water compression coefficient, μ w is the water viscosity, q w is the leakoff fluid, S w 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  9 Geofluids 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][44][45]: where δ n and δ s are, respectively, the normal and tangential displacement jump vector, δ eq are effective displacement, igure 13: Effect of stress difference on wellbore fracture width and injection pressure. 10 Geofluids T 0 n and T 0 s are, respectively, the normal and tangential stress component of the stress element under a linear elastic condition, < > denotes the Macaulay bracket, f c is a parameter evaluating the fracture failure, and f t 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][45][46] where d is the scalar damage variable between fractures, 0 < d < 1, and T n and T s 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 deter-mined by the Benzeggagh-Kenane fracture criterion [47]. The combined energy dissipated by failure G c is defined as where G c is the critical fracture energy release rate, GcI and GcII are the normal and first shear direction fracture toughness values, G I and G II are the normal and shear direction fracture energy release rate, and η is a material properties parameter; here is assumed to be 2.3.

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 as p n+1 = ωp n−1 + ð1 − ωÞp n ; the pressure distribution and  11 Geofluids 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.

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 as 50 × 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.
To 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 (w o ) and injection pressure solution (p w ) of KGD model can be given [48]: 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   Figures 2 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. Figure 4 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. Figure 5 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.

Geofluids
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.  14 Geofluids Figure 7 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.

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 Figures 12-16. Figure 12 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.

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 Section 4.2. The calculation results are presented in Figures 17-21. Figure 17 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. igure 20: Effect of dilation angle on pore pressure.

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 in x 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