## Abstract

A 3D non-linear fluid–solid coupling model for horizontal fracture of vertical well was established with the ABAQUS code. The wellbore, cement casing, perforation, pay layer and barriers were included in the model. Fluid–solid coupling elements were used to describe the behavior of formation stress–seepage flow coupling; pore pressure cohesive elements were employed to simulate the process of fracture initiation and propagation in formation. A typical horizontal fracturing process of a vertical well of Daqing Oilfield, China was simulated with the model. All the concerned parameters in simulation were taken from the field measurements.

#### Authors:

##### J. Zhang, F. J. Biao, S. C. Zhang and X. Wang

School of Petroleum Engineering, China University of Petroleum, Beijing 102200, China. Department of Modern Mechanics, University of Science and Technology of China.

Received: 4 November 2010, Accepted: 23 November 2011, Published online: 21 December 2011

© The Author(s) 2011.

The simulated bottom-hole pressure evolution is consistent with the data measured from the field. The configurations of the fracture and porous pressure distributions in the fracture are presented and discussed.

## Introduction

Hydraulic fracturing has been one of the most frequently implemented techniques for stimulating production of oil/gas reservoirs for several decades (Sneddon 1946).

Common water-based fracturing fluid with chemical additives is pumped at high rate and pressure into the formation during hydraulic fracturing process. When the pumping pressure exceeds the strength of the formation rock, fractures are induced and propagated into the formation, and then the propping agent is pumped into the fractures to keep them from closing after pumping pressure is released. Therefore, a man-made passage with high conductivity is constructed and hydrocarbon can flow into the well from the low-permeability formation (Economides and Nolte 2000).

Horizontal or vertical fractures of vertical wells may be developed depending on the factors such as depth of the well, distributions of in situ stresses and lithological parameters in formation. Fractures usually propagate perpendicular to the minimum principal stress. At shallow depths, the minimum principal stresses correspond to the overburden; therefore, horizontal fractures are constructed (Settari and Michael 1986). Numerical simulations of vertical fractures have been studied by many authors (Morales and Abou-Sayed 1989; Boone and Detournay 1990) and reviewed (Fragachan et al. 1993).

Most of the models fall into pseudo 3-dimension (p-3D) and are successful in design and performances of hydraulic fractures. Results of microseismic imaging show that fracture orientation is related to thrust fault near the wellbore. A horizontal fracture could be formed below the fault (Maxwell et al. 2009).

Engineering practices of Daqing Oilfield, China show that many of the hydraulic fractures are horizontal in areas of remained reserves. Significant re-stimulation could be gained from horizontal hydraulic fracturing (Schneider et al. 2007), but studies on mechanism and numerical simulation of horizontal fractures are rarely founded in the literature to our knowledge. A quasi static radial model of horizontal fracture was proposed (Sneddon 1946).

The relationship among width, radius and pressure in the fracture was presented, but flow and variation of pressure in the fracture as well as leak off across the fracture walls were not considered in the model. Interference among horizontal fractures of multi-pay layer reservoirs was studied (Zhang and Zhang 2004). However in the model plan, strain assumption was made, the length of the fracture was pointed to be a fixed value and the loading applied on the fracture walls distributed force but not fluid pressure.

In the present paper, a 3D non-linear fluid–solid coupling model for horizontal fracture was established with the ABAQUS code. The wellbore, cement casing, pay layer, barriers and perforations were included in the model. Fluid–solid coupling elements were used to describe the behavior of formation stress–seepage flow coupling; pore pressure cohesive elements based on damage mechanics were employed to simulate the process of fracture initiation and propagation.

A typical horizontal fracturing process of a vertical well in the Daqing Oilfield, China was simulated with the model. The simulated results show that the evolution of the bottom-hole pressure matches the data measured from the field very well. The correctness and reliability of the proposed model is validated. The horizontal fracture/configurations as well as porous pressure distributions in the fracture are also presented and discussed.

## Mathematical physics model

Hydraulic fracturing can be regarded as two physical processes: one is the fluid field in the well and fractures; the other is the stress–seepage flow coupling field in formation. The two processes are connected by fluid exchange and interaction of pressure on the fracture walls.

### Governing equations of stress–seepage flow coupling field in formation

In the current configuration, the equilibrium equation of porous formation can be written as (Zienkiewicz and Taylor 2005)

where σ, p_{w}, δε, δv and **f** are the effective stress matrix, pore pressure, virtual strain rate matrix, virtual velocity vector, surface force vector and volume force vector, respectively. The effective stress couples solid deformation with fluid flow by the following formulation (Economides and Nolte 2000)

where * σ* is the stress caused by solid deformation. Continuity equation of seepage flow considering large deformation of porous formation is (Malvern 1969

where J and x are the rate of volume dilatation and space vector, respectively; ρ_{w}, n_{w} and *v*_{w} are fluid density, porosity and seepage velocity, respectively.The kinetics equation is the Darcy’s law in the following form (Marino and Luthin 1982)

where * k* and

*are permeability matrix and gravity acceleration vector, respectively.*

**g**### Flow equations of fluid in fractures

The flow of proppant-laden fluid in the fracture is assumed to be incompressible Newtonian fluid and can be resolved as a tangential component along the cohesive element walls and a normal component across the cohesive element walls, respectively. The tangential one can be expressed as (Dean and Schmidt 2008)

where * q*,

*t*,

*μ*and

*p*are the vector of volume flow rate along the cohesive element walls per tangential unit length, the opening thickness of the cohesive elements, the coefficient of viscosity of fracturing liquid in the cohesive elements and fluid pressure in the cohesive elements, respectively. The normal component represents speed of fluid flowing in or out the walls of the cohesive elements. That is the filtration rate of fracturing liquid in engineering and expressed as (Hagoort et al. 1978)

where *q*_{t} and *q*_{b} are the volume flow rate across the top and bottom walls of the cohesive elements, respectively; *c*_{t} and *c*_{b} are the filtration coefficients of the top and bottom walls of the cohesive elements, respectively; *p*_{t} and *p*_{b} are the pore pressures at the top and bottom walls, respectively; and *p*_{i} is the fluid pressure along the middle plan of cohesive elements.

### Damage model of cohesive element

Initiation and extension of fracture are simulated with the cohesive element. The damage initiation of the formation can be expressed by the following criterion formula in cohesive element (Camanho and Dávila 2002)

where *σ*_{n} is the normal stress, *σ*_{s} and *σ*_{t} are the tangential stresses, and σ^{0}n is the intensity of tension of formation. σ^{0}s and σ^{0}tare the threshold stresses of tangential damage. The symbol < > denotes that only tensile stress could make the cohesive element damaged.For simulating degeneration of elastic modulus of cohesive element after damage initiated as expressed in Eq. (7), the following linear damage evolution rule is adopted.

where *E*^{0} and *E* are the initial element elastic modulus (without damage) and damaged element elastic modulus, respectively, and *d* is the damage factor which can be calculated with the following formula (Turon et al. 2006)

where *σ*_{eff} and *δ* are the effective tensile stress and corresponding strain, respectively. *G*_{0} and *G*^{C} are the elastic strain energy at damage initiation and at breaking of the element, respectively.

## Finite element discretization and the simulation model

The deformation, strain, stress and the fluid seepage in the formation during hydraulic fracturing are described by Eqs. (1)–(5) which couple each other and are highly non-linear. A corresponding incremental finite element formula was derived in detail (Zhang et al. <2010). The final expression of the formula and the meanings of the variables in the formula are listed in the following form

where Δ(* u*)

^{n+1}and Δ(

**p**_{w})

^{n+1}are incremental displacements and incremental pore pressure, respectively. The submatrices in the left side of the formula and load vector in the right side are complicatedly related to the current deformation and pressure.

The incremental finite element formula is solved with the well-known Newton–Raphson iteration scheme. The tolerance value of iteration residual in each increment is set as 0.5% of the force. The eight node hexahedron elements are used for describing deformation and seepage in formation. There are four nodal unknowns at each node, i.e., incremental displacements in three directions of the Cartesian coordinate system and incremental pore pressure. The plan cohesive elements with six nodal pairs are adopted to simulate the fracture initiation and propagation.

A typical vertical oil well in Daqing Oilfield, China is simulated with the proposed model. The wellbore, cement casing, pay layer, barriers and perforations were included in the model as shown in Fig. 1. The depth of the artificial bottom hole is 1,230.2m. A fracturing pipe is used to produce pipe string with 62 mm interior diameter. The outside diameter and thickness of the wellbore are 140 and 7.72 mm, respectively. The outside diameter of the cement casing is 200mm; the outside diameter and length of the perforation are 11 and 786 mm, respectively. The thickness of the pay layer is 5.9 m. The pay layer is interlinked by the top and bottom barriers, respectively.

**Fig. 1 **Local schematic of the model near the well hole

In situ stresses along *X*, *Y* and *Z* axes are 23, 26 and 20 MPa, respectively. The density of the reservoir is 2,200 kg/m^{3}. The initial porosity and initial porous pressure of formation are 0.20 and 10.8 MPa, respectively. The permeability of the reservoir is 0.069 μm^{2}. The normal and tangential damage threshold stresses of cohesive elements are 2.5 MPa. The proppant is quartz sand. Other related parameters are listed in Table 1.

**Table 1 **Formation lithological parameters

In simulation, the slurry is equivalent to a liquid with viscosity related to proppant concentration as (Adachi et al. 2007)

The sizes of the simulation model along *X*, *Y* and *Z* axes are 300, 300 and 100 m, respectively. The *XY* plane is located in the middle of the pay layer and *Z* axis is toward up. All normal displacements are restricted and pore pressure kept at 10.86 MPa on the model boundaries during the whole fracturing process. The finite element mesh includes 72,800 elements and 116,688 nodes.

Figure 2 depicts the sketch of the finite element mesh.

**Fig. 2 **Sketch of the finite element mesh.

The simulation of the fracturing process was carried out with the field operation curves including liquid injection rate, proppant concentration and injection duration as shown in Fig. 3. The major controllable quantities in hydraulic fracture practice are the injected pad liquid rate and proppant, which are pumped into the wellbore from their storage tanks. The quantities in operation are measured and recorded by instruments connected in a computer van. It can be seen that the measured wellhead pressure increases sharply before the break point due to the hardness of fracture initiation. The rock is damaged and the fracture propagates steadily after that point, the pressure decreases and then almost keeps unchanged.

Additionally as the proppant is injected into the fracture, the hydrostatic head increases while the perforation friction loss decreases. These factors also cause the wellhead pressure to decrease. The pad fluid without proppant is injected in the first several minutes of fracturing and then the proppant is added; the concentration increases in a step form. Therefore, a time lag between proppant concentration and wellhead pressure increase is induced as shown in Fig. 3.

**Fig. 3 **Field operational curves during fracturing.