Numerical model of fracture growth in hydraulic re-fracturing

Simulation of fracture propagation with FEM method requires re-meshing to provide more accurate results. This raises a question about the determination of the direction and criterion for mesh modification. In the case of general-purpose CAE-packages, we deal with a stationary mesh, and the fracture path is usually represented as a chain of elements with degraded properties. The algorithm proposed in this paper is based on the ANSYS Mechanical APDL language for stepwise geometry reconstruction and mesh modification in accordance with the current configuration of a growing fracture and provides a more accurate description of its shape. The fracture propagation process is divided into stages. Each subsequent stage differs from the previous one by the fracture shape modified due to the crack length increment in the calculated direction. To check the adequacy of the model, an experiment on fracture propagation in glass specimens with an initial notching under uniaxial compression was performed. The laboratory experiments were carried out to determine the fracture toughness of rocks. The developed numerical model has been used to solve the problem of refracturing for different stress anisotropy in the oil-bearing rock formation.


INTRODUCTION
n the last few years the effect of the directional hydraulic re-fracturing (HF) on the stress-strain state of oil-bearing rock formations has been extensively discussed in the scientific and engineering literature. Thus, as it was stated and verified by the specialists of OJSC "NK ROSNEFT' in [1,2], the repeated HF leads to the formation of fracture, the direction of which differs from the direction of the primary HF. The fact of reorientation of the secondary hydraulic fracture with respect to the primary HF has been also supported by foreign specialists [3,4]. The application of the directional hydraulic re-fracturing method offers much promise from the viewpoint of increasing the recovery of oil pools especially in the oil fields, in which the majority of wellbores has been already subject to primary hydraulic fracturing. As applied to these cases, it is of practical importance to analyze thoroughly an expected change in I the direction of the secondary HF and the prospects of using its positive effects for increasing the wellbore capacity. The practical implementation of the secondary HF method suggests the development of a new direction in the hydraulic fracturing technique. As a consequence, an essential improvement of the efficiency of oil wells and accordingly the growth of the oil recovery coefficient of the entire oil fields, which demonstrate a notable decrease in the production rate, is expected. The mathematical modeling of the process of fracture growth relies both on the analytical [5] and numerical [6] methods. The study, which is close to the present work with respect to the set of hypotheses used for problem consideration, is reported in paper [7] where the problem of HF growth was solved using author's version of the hyper-singular boundary element method. At present, there is a series of customized commercial software packages for numerical simulation of fracture propagation allowing for mesh adaptation including the case of three-dimensional problem formulation. Thus, work [8] offers a 3Dmodel to describe the growth of a crack in machine elements, which was implemented numerically as a finite element ZENCRACK program (developed by the company Zentech International Limited) applied in fracture mechanics. Another example of adaptive remeshing to trace crack propagation is the Franc3D program (FRacture ANalysis Code for 3D), which is based on the boundary element method. The possibilities of this package are demonstrated in work [9] in the context of the problem of crack growth in the gear wheel. Both packages are intended for needs of mechanical engineering. The main constraint in finding the finite-element solution to the problem with the aid of generally recognized universal finite-element packages is the necessity of reconstructing the mesh to provide the desired accuracy in the evaluation of the fracture growth direction and fracture growth criterion. In the universal CAE software, in which the mesh remains invariable, the fracture path is traditionally represented as a chain composed of the elements with degraded properties (see, for example [10]). This necessitates a search for an optimal mesh topology and essentially decreases the solution accuracy. The algorithm proposed in our work uses the facilities of ANSYS Parametric Design Language in the framework of ANSYS Mechanical Package. It allows a step-wise reconstruction of the mesh according to the current configuration of the computational domain and provides the most accurate description of the shape of the growing fracture and the stress-strain state in its neighborhood.

STRESS FIELD AFTER FORMATION OF THE PRIMARY HF FRACTURE AND INITIATION OF THE SECONDARY HF FRACTURE IN THE ORTHOGONAL DIRECTION
fter completion of the primary HF the properties of the rock formation in the neighborhood of the oil well change. The formation of the fracture improves the general picture of rock permeability causing a redistribution of internal loads, which essentially changes the stress-strain state (SSS) of the rock formation. In view of this fact a series of computations have been done to get answers to the two main questions: how does the fracture pressure change under the conditions of the secondary HF and in what direction does the secondary hydraulic fracture propagate? The fracture geometry was described by the Perkins-Kern model [11,12]. It was based on the assumption that the fracture cross section in the plane normal to the axis of the well is an ellipse and its maximum width is detected in the nearwellbore region. This is the area adjacent to the well where the filtering characteristics of the productive strata have changed as a result of physicochemical processes initiated by processing method and operating conditions. The general computational scheme and the system of coordinates are shown in Fig. 1. Then, we get a fracture of half-length 50 f x  m, the direction of which coincides with the direction of the highest horizontal stresses H  The stress state of the wellbore as a function of the fracture width w was analyzed at different ratios of maximum stress ( H  ) to minimum stress values( h  ). The problem was solved by the finite element based on the ANSYS package under the assumption that the mechanical properties of the proppant (the material used to secure the main hydraulic fracture) are similar to those of the rock formation. Simulation of the existing hydraulic fracture involves a preset of fracture wing displacements according to a maximum opening width w and variation of w with a horizontal distance x according to the laws [13]: 1 4 ,0 3,57 , A where  is the fracture liquid viscosity; q the rate of liquid pumping;  is the plane strain modulus of the rock ( E -is Young's modulus;  is Poisson ratio). As one might expect the FEM computations demonstrate the growth of tangential stresses   on the well walls with increase of the opening width of the primary HF. Fig. 2 presents the results of six calculations for different values of maximum fracture opening.  The maximum of incremental stresses   is located on the walls of the well. Then, the stresses in the rock formation rapidly decrease in the radial direction until they are stabilized at a distance of 2-3 m from the wellbore, namely, at 20 30 c r    , where  is the current radius and c r is the radius of the wellbore (Fig. 3). Based on the characteristic features of the incremental stress distribution   we can conclude that in the isotropic field of horizontal stresses H h    it is more probable that re-fracturing will occur in the direction perpendicular to the primary HF, i.e., the secondary HF in the isotropic field of horizontal stresses is orthogonal to the primary fracture. In this case, the pressure increase can be observed at the moment of fracture initiation. Thus, at the pressure of the primary fracture of 80.0 MPa and the width of the primary fracture opening near the wellbore equal to 20 mm the secondary fracture pressure is 92.2 MPa. Essentially, the situation changes in the case when the minimum stress h  is about 70-95% of the maximum stress 47 H   MPa. Irrespective of the appearance of incremental stresses   , in all examined cases the most likely direction of the secondary HF is the direction parallel to the primary fracture path. As in the first case, the pressure begins to grow at the time of initiation of the secondary fracture and the primary fracture opening increases. The degree of pressure growth depends on the anisotropy of the stress field. Thus, the width of the primary fracture opening 20 w  mm at 47

 
MPa by a factor of 1.37. The deformation of porous liquid-saturated rocks is a complicated process, during which the deformation of the mineral matrix of rocks (under the action of varying effective stresses and fluid pressure gradients) occurs simultaneously with fluid filtration in pores (under the action of the fluid pressure gradients and bulk deformation of the skeleton). In the following, a drop of oil pressure due to extraction of fluid from the well with the primary HF causes a change in the stress state of the rock formation in the vicinity of the well. Therefore, an adequate treatment of the processes accompanied the deformation of a porous liquid-saturated medium requires that the system of differential equations describing the deformation of the rock skeleton and fluid filtration should be considered simultaneously. The solution of this coupled problem was found numerically using the ANSYS software package. To simulate the process of liquid removal from the wellbore, we performed coupled calculations for SSS and liquid filtration at the prescribed bottom-hole pressure. The finite element mesh was built using the CPT213 element, which makes it possible to implement numerically the filtration consolidation model [6]. By virtue of the problem symmetry we performed computations only on a quarter of the computational domain shown in Fig. 1. The computational 500 500  m domain reproduced the wellbore and the existing primary HF filled with proppant. The permeability of the proppant was taken to be equal to 200 D, the permeability of collector was 100 µD. These values are typical for the local oil fields. The problem was solved in two steps. At the first step we investigated the stress state of the rock formation after the occurrence of the primary hydraulic fracturing, which implies that the problem was solved at initial stresses 47 H h     MPa and initial oil reservoir pressure of 20 MPa. The values for the problem under consideration were taken from field data. At the second stage a bottom-hole pressure of 12 MPa was applied and the operation of the well under such conditions was simulated within the period of three years. The time factor is accounted for explicitly in the coupled simulation in ANSYS. The calculated distribution of the oil bed pressure in the vicinity of the well was represented in the form of ellipse due to the influence of the primary HF. The simulation results show that the well exploitation is accompanied by a redistribution of the pore pressure, which in its turn leads to a change in the stress-strain state of the surrounding formation. Near the fracture and well the full stresses are reduced to 44.5 MPa , that is, become smaller by 2.5 MPa against the initial 47 MPa. Thus, pumping of liquid from the well with the primary HF leads to a reduction of horizontal stresses in the near-well region. However, in the present work the horizontal stresses were calculated taking no account of the above-mentioned incremental stresses, which are generated by the primary HF. The results of computation taking into account both of these factors are shown in Fig. 4 in the form of curves of horizontal stresses x  in the vertical cross-section over the well bore (the maximal opening of the primary HF near the well was set equal to 20 mm). As is evident from the plots, consideration of incremental stresses caused by opening of primary fracture leads to a growth of horizontal compressive stresses. The observed growth turns to be much greater than a reduction of stresses initiated by the drop of oil pore pressure caused by liquid pumping out. The growth of stresses is most pronounced in the vicinity of fracture but as the distance from the fracture increases the stresses decrease drastically and at a distance of 65 m from the fracture the stress curves practically merge together. The calculations show that in the vicinity of fracture the horizontal stresses reach 60 MPa, but at a distance of 2 m their value does not exceed 48 MPa, which is by only 1 MPa higher than the initial values (47 MPa dashed line in Fig. 4). These results show that after carrying out the directional hydraulic fracturing the main effect on the redistribution of stresses in the vicinity of well is produced by the deformation of rock formation caused by generation of the primary fracture rather than by the operation of the well. In this case, the main conclusion will be as follows: the growth of the secondary HF in the direction perpendicular to the path of the primary fracture is possible only in the presence of isotropic or sub-isotropic horizontal stresses. The stress-strain state of the Priobskiy oilfield, in which the effect of reorientation of the secondary HF was recorded, is characterized by a rather small difference between the minimum and maximum horizontal stresses. In particular, the acoustic anisotropy does not exceed 5%. Under the conditions of anisotropic horizontal stresses, the secondary HF is initiated by application of technical means allowing engineers to create a weakening surface at a certain depth [1] and in the desired direction.

EVOLUTION OF RE-FRACTURING IN THE ANISOTROPIC STRESS FIELD
he computational scheme for this case is similar to that of the previous problem (see Fig. 1). The only difference is that in the direction of the lowest horizontal stresses h  the rock formation is subject to the condition of preweakening by the initial fracture of half -length  [11][12][13]. Upon substituting the opening width Eqn. (1) into the flow law with allowance made for the elliptical shape of the fracture and integrating the resultant expression we get the law of pressure distribution in the secondary fracture: where 0 (0, ) In formula (2) the time is taken into account implicitly, since the half-length f X changes with time. Hence, by solving the flow equation at each instant of time we find the growing length of the secondary fracture and pressure distribution inside it. Here we used the following assumptions:1) leakage of the fracture liquid into the oil bed is ignored; 2) the fracture remains rectilinear, since the PKN model provides for the existence of a straight fracture. In other words, the pressure distribution was taken from a model for a straight fracture although the fracture in this case was not straight. The version of the PKN model (originally 3D) used by the authors is based on the assumption that the profile of the average horizontal cross section remains invariable throughout the height of the fracture. This introduces a certain error in the computations, which decreases with increasing oil bed thickness. Opening of the fracture (at most to a width of 20 mm) is significantly smaller than the height of the fractured oil bed (more than 2 m) and the length of the fracture (extending as long as 20 m) In this case, in the first approximation the computation error for 2D case can be considered acceptably small compared to the error in the 3D formulation. The obtained equations were solved for average parameters typical of the oilfields developed in the Perm region (Tab. 1). Since the existence of the primary fracture changes the picture of the initial SSS, the wellbore pressure is expressed as

Algorithm for numerical implementation of the secondary fracture growth model
The proposed algorithm for numerical computation of the secondary fracture growth relies on the following hypotheses: Hypothesis I. Fracture initiation necessitates the fulfillment of the following criterion: where Par is the criterion of the fracture stability, where [ ] Par is its critical value. For criteria we can take the principal stress at the fracture top 1  the stress intensity coefficient I K and the angle of fracture opening cr  . Hypothesis II. The direction of fracture propagation coincides with the vector of the normal to the first principal stress at the fracture top directed at the smallest (in the absolute magnitude) angle (energy-optimal) to the current direction of the fracture. The process of the secondary HF growth is simulated in a step-wise manner. Each subsequent step differs from the previous one by a changed topology of the secondary fracture caused by its extension in a specified direction to a specified length. The solution is sought for the problem of linear elasticity under the assumption of small strains The behavior of the fracture is described at each step by the system of standard equations of the stationary problem of elasticity (div 0 is the total strain tensor, 4Ĉ is the tensor of elastic constants, , u S S  are the parts of the boundary with prescribed displacements and loads, respectively. For numerical implantation of the finite element model we used the ANSYS Mechanical APDL Consider the algorithm, which allows us to analyze the secondary fracture growth based on the limiting value of its opening: The simulation of fracture propagation at the k-th time step consists of the following steps: 1. Calculate the stress-strain state for the current configuration of the computational domain, bearing in mind that the length-wise distribution of pressure is governed by law (2) 1 4 where L is the current length of the fracture; w ; p  is a variable quantity. Verify the possibility of fracture growth by controlling the fulfillment of inequality (3). Note that at the first step the initial value of 0 k p  is calculated as 0 0.2488 1 6.1873 p L   and at subsequent steps -as 0 . By the opening angle ( ) cr L  we mean the angle in the polar coordinate system with its origin lying at the top of the fracture, which is generated as a result of fracture opening under pressure between the lines connecting the fracture top with two adjacent nodes of the finite-element mesh on its wings (Fig. 5). In the figure, the geometry of the fracture near the top in the unloaded state is shown by solid lines and dashed line denote its geometry in the loaded state; the circles denote the position of the mesh nodes. Vector cr x determines the current direction of the fracture. The opening angle due to smallness of the displacements is calculated by the formula   (3) is not fulfilled to a specified accuracy, then the variable parameter k p  is calculated by the chord method from expression (2) on the invariable mesh at the k-th time step according to the following integral scheme: 2.1. Scale the value of the pressure increment for the next iteration using formula

If condition
where crij x is the j -th component of the vector crj In the example given in fig. 6 the fracture will grow in the direction 2 cr cr  x x . 4. Reconstruct the geometry of the sub-region around the fracture. Here, it is necessary to exclude the region itself and its mesh (by means of commands ACLEAR, ADELE), leaving only the boundary of the sub-region. Extend the part of the boundary adjacent to the fracture top. The increment of the fracture length at the k-th step is equal to k l  and in the general case it can be a variable quantity. Construct a sub-region within a newly generated closed contour, generate a mesh on its surface and set new boundary conditions. 5. Calculate the stress-strain state of the rock formation, in which the fracture has a new configuration 6. Repeat the algorithm item 1-4 N times The calculated total length of the growing fracture is In the case when the critical value of the stress intensity coefficient К I , is taken to be the fracture growth criterion, the generation of the finite element mesh suggests the usage of a special type of element arrangement-singular finite element meshing (Fig. 7).

Experimental verification of the proposed algorithm
To verify the adequacy of the accepted hypotheses we performed a fracture experiment to investigate the growth of fractures in glass specimens under uniaxial tension, by analogy with the experiment described in monograph [5]. The specimens were made of 200 100  mm window glass, 4 mm thick and had straight centrally located grooves 2.5 mm thick, which were cut by the abrasive water jet cutting machine at an angle of 30° to the horizontal axis of the specimens. Actually, we repeated the experiment of Hoek and Brown in which the growth of fractures in rock formations under the action of compressive stresses was studied. The compression tests were performed in the LFM-50T test machine with grips 96 mm wide and acceptable load limit of 5 ton-force (Fig. 8). To reduce inhomogeneities of the pressure distribution, we placed the polyurethane inserts 1 mm thick in the zones of sample contact with the grips.   Fig. 9a. The fractures generated in almost all specimens are located both above and below the calculated fractures. In sample № 2 the lower fracture was initiated in the middle of the sample due to the presence of the local stress concentrator. First, in all cases the angle between the primary and secondary fractures remains unchanged and is equal to 85°. As the fracture grows further the opening angle first slightly decreases and then begins to increase. In this case the fracture shows the tendency to change its direction to a vertical one. The fracture paths were digitized and interpolated by the piece-wise linear functions (no less than 30 nodes for one fracture). Owing to the specimen symmetry the fractures initiated at both ends of the initial cut were superimposed. Then at each of 200 nodes of a thicker uniform mesh along the horizontal axis i x a mathematical expectation    was determined (here i is the node number and j is the test number). As is evident from Fig. 9b, the relative discrepancy between the calculated and actual fracture paths does not exceed 5%, which substantiates the adequacy of the proposed algorithm for computation of the growing fracture path.

Experimental determination of fracture toughness coefficients of rock formations
To determine the value of the intensity coefficient К 1С we tested the specimens 30 mm in diameter and 60 mm in length made from the core material extracted from the wellbore of the 118 Enapaevskiy oilfield of "LUKOIL -Perm" CoLtd from the terrigenous interval of 1538-1541 m. Prior to tests of determining the critical coefficient of stress intensity we evaluated the static and dynamic elastic properties for some specimens in reservoir conditions. Schematic diagram of testing cylindrical specimens with o-ring groove are presented in Fig. 10. This testing scheme was chosen from [15] and adapted was rock samples. This type of test is convenient because core samples of standard size can be used. The size of a notch was 7 mm. Tensile tests of rock specimens were carried out in the Instron ElectroPuls E10000 electrodynamic two-axis test system, which is designed for compression, tension, bending and torsion static tests; for dynamic fatigue tests at the frequency of 100 Hz; for two-axis (tension-compression) static and dynamic tests under loads up to 10 kN/100 Nm. Testing specimens of this type required the development and manufacture of special tooling for securing cylindrical specimens in the grips of testing machine during tensile tests. Tensile tests were carried out at room temperature and at the prescribed speed of grips displacement, which was 0.2 mm/min. a b Figure  The stress intensity factor К 1С was calculated by the formulas presented in work [15] for a cylindrical specimen with a surface o-ring fracture subject to tensile stress. Here the notation corresponds to that of Fig. 10. We tested 25 specimens and based on the obtained results calculated the mean value of К 1С , which was equal to 0.108 MPa·m 1/2 at the maximum and minimum values of 0.044 MPa·m 1/2 and 0.168 MPa·m 1/2 , respectively. The tensile strength of the specimens  р was found to vary from 0.55 MPa to 2.34 MPa, which are consistent with the values typical for the rocks of this kind. Based on the results of our tests we determined the correlation relationship between the critical stress intensity coefficient and the static modulus of elasticity under the oil bed conditions (Fig. 11). It is clear that further testing is needed to improve the correlation since the amount of tested samples was quite low. The obtained experimental data did not reveal the dependence of К1C and р on the porosity, velocity of longitudinal wave and dynamic elastic properties of specimens.

THE RESULTS OF SIMULATION OF THE SECONDARY HF GROWTH IN THE ANISOTROPIC STRESS FIELD
he developed numerical model and the results of experimental investigation of К 1C were used to solve the problem of the secondary fracture growth at the following values of the coefficient of the stress field anisotropy a K and the width of the primary fracture opening ,0 w w . The calculations were done at IC K =0.127, which is the mean of the tested rocks. The results of computation are shown in Fig. 12. As it is seen from the figure, with increasing width of the primary fracture opening, the secondary HF practically does not change its direction towards the highest compressive stress even at large stress anisotropy and grows in the direction perpendicular to the primary fracture. The main reason of such behavior is that at large width of the primary fracture opening the level of the initial stress field anisotropy in the T vicinity of fracture decreases. This causes an increase in the pressure minimum, at which the conditions for the primary fracture growth are generated. The direction of the secondary HF depends not only on the width of opening of the primary fracture, but also on the coefficient of stress anisotropy a K . The larger is the opening width of the primary fracture and the closer is the anisotropy coefficient to 1.0 the nearer is the path of the secondary HF to the normal to the primary fracture. Thus, at the width of the primary fracture opening of 5 mm the direction of the secondary HF is close to the normal at 0.85 a K  , and at the opening width of 20 mm this happens at 0.7 a K  .

CONCLUSIONS
he applications of the developed numerical model of the fracture growth under the conditions of hydraulic refracturing in the oil bed have shown that the following factors are responsible for the generation of the secondary fracture orthogonal to the primary fracture. 1. Anisotropy of the medium, in which the fracture is propagating (the coefficient of stress anisotropy a K should be higher than 0.8) 2. An increase the pumping pressure of fracture liquid 3. An increase of the opening width of the primary fracture In the presence of all these factors the effect of the secondary HF will be maximal because its realization opens up the maximum of non-depleted zones of the oil pool. Our calculations were carried out based on the assumption that during hydraulic re-fracturing the fracture liquid does not penetrate into the primary fracture, which means that it is necessary to provide isolation of the primary fracture while executing re-fracturing. Otherwise, it is probable that its opening proceeds -the conclusion that has been also arrived at by foreign researchers investigating the growth of the secondary HF in the framework of physical models [4]. Furthermore, as was already mentioned, to resist the additional circumferential stresses caused by generation of the primary HF, it is necessary to create weakening surfaces at certain depth and in the certain direction. However, at equal components of the horizontal stresses the development of the secondary HF in the direction perpendicular to the primary HF is possible without artificial technical interference.