Microsoft Word - numero_49_art_16_2454 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 140 Focused on Russian mechanics contributions for Structural Integrity Numerical model of fracture growth in hydraulic re-fracturing Oleg Yurievich Smetannikov, Yuriy Aleksandrovich Kashnikov, Sergey Gennadievich Ashikhmin, Artem Eduardovich Kukhtinskiy Perm National Research Polytechnic University, 614990, Komsomolskiy av, 29, Perm, Russia sou2009@mail.ru, geotech@pstu.ru, a_s_g_perm@mail.ru, artyom@pstu.ru ABSTRACT. 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 re- fracturing for different stress anisotropy in the oil-bearing rock formation. KEYWORDS. Fracture propagation; FEM; stress anisotropy; fracture path; APDL. Citation: Smetannikov, O. Y., . Kashnikov, Y. A, Ashikhmin, S. G., Kukhtinskiy, A. E., Numerical model of fracture growth in hydraulic re-fracturing, Frattura ed Integrità Strutturale, 49 (2019) 140-155. Received: 03.04.2019 Accepted: 22.05.2019 Published: 01.07.2019 Copyright: © 2019 This is an open access article under the terms of the CC-BY 4.0, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. 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 http://www.gruppofrattura.it/VA/49/2454.mp4 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 141 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 3D- model 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 near- wellbore 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 50fx  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 , f w q x w E       1 4 0 ,0( ) 1w f x w x w x       , (1) A O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 142 where  is the fracture liquid viscosity; q the rate of liquid pumping;  21E E    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. Figure 1: Computational scheme for the problem of secondary oriented HF in rock formation. Figure 2: The incremental growth of incremental stresses  on the well walls depending on the value of the opening angle of the primary hydraulic fracture in the isotropic field of horizontal stresses ( H h   ) at different values of maximum opening w , mm: 3 (symbol ●); 5 (■);7,5 ();10 (♦); 15 (); 20 (□). O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 143 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 30cr   , where  is the current radius and cr is the radius of the wellbore (Fig. 3). Figure 3: Radial distribution of incremental stresses  in the isotropic field of horizontal stresses ( H h   ) at the width of opening of the primary HF 20w  mm. 1 — in the direction of H , 2 — in the direction of h . 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 47H  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 20w  mm at 47H h    MPa is responsible for a maximal increase of the secondary fracture pressure by a factor of 1.15, whereas at 47H  MPa, 32.9h  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 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 144 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 47H 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). Figure 4: Distribution of radial stresses along the trace of the plane crossing the well axis: taking into account liquid removal (solid line) and taking into account both the removal of liquid and opening of the primary fracture (dashed line) in the isotropic field of horizontal stresses: general picture (a) and its magnified fragment (b). 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 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 145 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 pre- weakening by the initial fracture of half -length 2fX  m, whereas in the direction of the highest horizontal stresses H the generation of fracture of half- length 50fx  m has been already completed. The value of the lowest horizontal stresses was prescribed as h a HK   , where aK is the coefficient of the SSS (stress-strain state) anisotropy. The value of the anisotropy coefficient varied from 0.7 to 1.0 and the width of fracture opening in the vicinity of the well varied from 5 to 20 mm. The distribution of the liquid pressure in the primary HF in the course of its evolution was determined based on the Perkins- Kern Nordgren (PKN) model of hydraulic fracturing [11-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: 0 1 1 ,f f x P P A X               (2) where 0 (0, )HP p t    ;  30256 .fA X q hw   In formula (2) the time is taken into account implicitly, since the half-length fX 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 0 H HP p      . The value of correction p is estimated by the independent calculation. By setting different values of pumping periods we can determine the parameters of the secondary HF – its length, opening width, well-head pressure, fracture pressure profile and also the parameters p , A . T O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 146 Parameter Value Wellbore diameter,m 0.216 Elasticity modulus of the rock GPa 40 Poisson’s ratio 0.25 Tensile strength 3.0 Initial oil bed pressure – 0p , MPa 17 Vertical rock formation pressure – V , MPa 47 Lateral rock pressure - H , MPa 47 Rate of liquid feeding in the fracture, m3/min 3.0 Viscosity of the liquid in the fracture, cPs 100 Fracture height – h , m 20 Table 1: Input data for simulation of the secondary HF. 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: [ ]Par Par , 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 IK 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  , Vx ;   T1ˆ 2     u u , Vx ; 4 ˆ ˆˆ C   ) with the following boundary conditions: u U , uSx ; ̂  n P , Sx . Here,  ˆ , t x is the stress tensor,  , tu x is the displacement vector,  ˆ , t x is the total strain tensor, 4Ĉ is the tensor of elastic constants, ,uS 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: ( ) [ ].cr crL   (3) 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 1( , ) 1 1 ,k w x P L x P A L             O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 147 where L is the current length of the fracture; w НР р    (where 47H  MPa); 0.2487 4.4611A L ; 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 0kp is calculated as 0 0.2488 1 6.1873p L  and at subsequent steps - as 0 0 1k kp p    . 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 crx determines the current direction of the fracture. The opening angle due to smallness of the displacements is calculated by the formula cr r nr l nlu l u l   , where ru  , lu  are the circumferential components of the displacement vectors of nodes adjacent to the top and lying on the right and left wings of the fracture; nrl , nll are the distances from the top to the corresponding nodes. Figure 5: The scheme of computation of the fracture opening angle. If condition (3) is not fulfilled to a specified accuracy, then the variable parameter kp 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 1j jk kp a p    with a prescribed coefficient 1.5a  . Here the subscript corresponds to the time step number and superscript - to the iteration number. 2.2. Solve the problem of determining the SSS for modified boundary conditions and find the value of the fracture initiation parameter 1jcr  . 2.3. Using jkp , 1j kp  , jcr , 1j cr  , determined in the previous iterations find the subsequent approximation 2jkp  : 1 2 1 1 1 j j j j jk k crk k j j cr cr p p p p                 , after which evaluate the stress-strain state and the corresponding value of the fracture initiation parameter 2jcr  . 2.4. In the case of non-fulfillment of condition 2 [ ] [ ] j cr cr cr       , where  is the prescribed accuracy, increase the iteration number by 1 and repeat the steps 2.1–2.4. O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 148 Figure 6: Computation of fracture growth direction. 3. Specify the direction of the first principal stress at the fracture top node, at which the intensity of stresses is maximal. As an example, fig. 6 shows its direction at the top of the fracture and also the current orientation of the fracture (dot and dash vector) and possible paths of its growth (vectors 1crx and 2crx ). The ANSYS Parametric Design Language macros contains following steps 3.1. Determine 6 components of the tensor of stresses at the top using the function *GET and load them in the six- element COMP file:  COMP 11 22 33 12 23 13, , , , ,       . 3.2. Determine the direction cosines of the vectors of principal stresses with respect to the global Cartesian coordinate system using the function *VFUN,DIR,DIRCOS, and locate them in the 9-elements DIR file:  DIR 11 12 13 21 22 23 31 32 33, , , , , , , ,c c c c c c c c c , ijc is the cosine of the angle between the i-th principal stress and the j-th coordinate vector. 3.3. Determine possible directions of the fracture propagation vector (Fig. 6):  1 12 11,cr c c x and its inverse  2 12 11,cr c c x . 3.4. Determine the coordinates of the direction vector at the end of the currently existing fracture (Fig. 6)  0 0 0 0 0,cr e b e bx x y y  x 3.5. Calculate the angles 1 , 2 between the vectors 0crx and 1crx , 0crx and 2crx , respectively:  1 01 2 02 0arccosi cri cr cri crx x x x l     ( 1, 2i  ), where crijx is the j -th component of the vector crjx ( 0,1, 2i  ; 1, 2j  ); 2 2 0 01 02( ) ( )cr crl x x  is the length of the vector 0crx . 3.6. According to the hypothesis II, select a directional vector of further advance of the fracture crx , lying at the smallest (in absolute magnitude) angle with the previous vector 0crx : 1 1 2 2 1 2 , ; , . cr cr cr         x x x In the example given in fig. 6 the fracture will grow in the direction 2cr 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 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 149 boundary adjacent to the fracture top. The increment of the fracture length at the k-th step is equal to kl 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 0 1 N k k L L l     meters 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. Figure 7: Types of finite element mesh with a singular block (cloud) for calculation of the secondary HF growth: general view and scaled up fragments. O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 150 Figure 8: Position of the sample during tests. The computations were done for the following glass properties and opening angle criterion 73E  GPa, 0.3  and 0.8  °. The size of the fracture length increment was a variable quantity and was equal to 0.5kl  mm for the first 10 steps of the algorithm and 1mm for all subsequent steps. The calculated and experimental configurations of the fracture are shown in 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 ix a mathematical expectation   8 1 1 8 ji i j y y    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 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 151 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 9: Comparison of experimental fractures in specimens and calculated fractures (thin solid lines) (a); the results of digitalization of fracture paths (b): line 1 – experiment; line 2 – experimental average, line 3 – calculation. 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. 12 ( )Ic P K b f b     , O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 152 b r   , (4) 2 3 4 1 0, 5(1 0, 5 0, 375 0, 363 0, 731 ) 1f            Here the notation corresponds to that of Fig. 10. Figure 10: Diagram of loading of the rock specimen with o-ring fracture for estimating К1C [15] Figure 11: Correlation relationship between К1C and the static modulus of elasticity under the oil bed conditions We tested 25 specimens and based on the obtained results calculated the mean value of К1С, which was equal to 0.108 MPa·m1/2 at the maximum and minimum values of 0.044 MPa·m1/2 and 0.168 MPa·m1/2, respectively. The tensile O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 153 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. Figure 12: The calculated paths of the secondary HF at different widths of the primary HF opening and different values of the anisotropy coefficient of the initial stress field in the vicinity of the wellbore (a) ,0 20ww  mm; (b) ,0 5ww  mm. 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 aK and the width of the primary fracture opening ,0ww . The calculations were done at ICK =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 O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 154 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 aK . 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.85aK  , and at the opening width of 20 mm this happens at 0.7aK  . CONCLUSIONS he applications of the developed numerical model of the fracture growth under the conditions of hydraulic re- fracturing 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 aK 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. REFERENCES [1] Latypov, I.D., Borisov, G.A., Haidar, A.M., Gorin, A.N., Nikitin, A.N., Kardymon, D.V. (2011) Reorientation refracturing on RN-Yuganskneftegaz LLC oilfields, Neftyanoe khozyajstvo – Oil Industry, 6, pp. 34–38. [2] Latypov, I.D., Fedorov, A.I., Nikitin, A.A. (2013). Research of reorientation refracturing, Neftyanoe khozyajstvo – Oil Industry, 10, pp. 74-78. [3] Wright, C.A., Conant, R.A., Stewart, D.E. Emanuel, M.A., Wright, W.W. (1994). Reorientation of propped refracture treatments, SPE paper 28078 presented at the 1994 SPE / ISMR Rock Mechanics in Petroleum Engineering Conference, Delft, Aug. 29-31. DOI: 10.2118/28078-MS. [4] Lan, Zh., Zhang, G., Spe, Hou, F., He, X., Liu X. (2008). Evaluation of refracture reorientation in both laboratory and field scales, SPE International Symposium and Exhibition on Formation Damage Control, 13-15 February, Lafayette, Louisiana, USA, SPE-112445-MS. DOI: 10.2118/112445-MS. [5] Parton, V.Z., Morozov, E.M. (1985). Mekhanika uprugoplasticheskogo razrusheniya [Elastic-plastic fracture mechanics], Nauka, Moscow. [6] Fadeev, A.B. (1987). Metod konechnykh ehlementov v geomekhanike [Finite element method in geomechanics], Nedra, Moscow. [7] Zubkov, V.V., Koshelev, V.F., Lin'kov, A.M. (2007). Numerical modeling of hydraulic fracture initiation and development, J. Min. Sci., 2007, vol. 43, no. 1, pp. 40-56. DOI: 10.1007/s10913-007-0006-6 [8] Timbrell, C., Maligno, A., Stevens, D. (2013). Simulation of complex 3D non-planar crack propagation using robust adaptive re-meshing and radial basis functions. Available at: http://www.zentech.co.uk/download/zentech-blos-nwc13.pdf. [9] Kramberger, J., Flasker, J. (2006). Numerical simulation of 3-D crack growth in thin-rim gears. Available at: http://www.gruppofrattura.it/ocs/index.php/esis/CP2006/paper/viewFile/9516/6139. T O. Y. Smetannikov et alii, Frattura ed Integrità Strutturale, 49 (2019) 140-155; DOI: 10.3221/IGF-ESIS.49.16 155 [10] Korolev, I.K., Petinov, S.V., Freidin, A.B. (2009). Numerical simulation of damage accumulation and fatigue crack growth in elastic materials, Vycisl. meh. splos. sred – Computational Continuum Mechanics, 2(3), pp. 34-43. DOI: 10.7242/1999-6691/2009.2.3.21 [11] Perkins, T.K., Kern, L.R. (1961). Widths of hydraulic fractures, J. Petrol. Technol., 13(9), pp. 937-949. DOI: 10.2118/89-PA. [12] Nordgren, R.P. (1972). Propagation of a vertical hydraulic fracture, Soc. Petrol. Eng. J., 12(4), pp. 306-314. DOI: 10.2118/3009-PA. [13] Economides, M. J., Oligney, R., Valko, P. (2002). Unified fracture design, Orsa Press, Alvin, TX. [14] Hoek, E. (1965). Rock fracture under static stress conditions, PhD Thesis in Philosophy and Engineering, The Faculty of Engineering of the University of Cape Town. [15] Murakami, Yu. (1990). Handbook of the stress intensity factors [Spravochnik po koehfficientam intensivnosti napryazhenij], 1, Mir, Moscow. << /ASCII85EncodePages false /AllowTransparency false /AutoPositionEPSFiles true /AutoRotatePages /None /Binding /Left /CalGrayProfile (Dot Gain 20%) /CalRGBProfile (sRGB IEC61966-2.1) /CalCMYKProfile (U.S. Web Coated \050SWOP\051 v2) /sRGBProfile (sRGB IEC61966-2.1) /CannotEmbedFontPolicy /Error /CompatibilityLevel 1.4 /CompressObjects /Tags /CompressPages true /ConvertImagesToIndexed true /PassThroughJPEGImages true /CreateJobTicket false /DefaultRenderingIntent /Default /DetectBlends true /DetectCurves 0.0000 /ColorConversionStrategy /CMYK /DoThumbnails false /EmbedAllFonts true /EmbedOpenType false /ParseICCProfilesInComments true /EmbedJobOptions true /DSCReportingLevel 0 /EmitDSCWarnings false /EndPage -1 /ImageMemory 1048576 /LockDistillerParams false /MaxSubsetPct 100 /Optimize true /OPM 1 /ParseDSCComments true /ParseDSCCommentsForDocInfo true /PreserveCopyPage true /PreserveDICMYKValues true /PreserveEPSInfo true /PreserveFlatness true /PreserveHalftoneInfo false /PreserveOPIComments true /PreserveOverprintSettings true /StartPage 1 /SubsetFonts true /TransferFunctionInfo /Apply /UCRandBGInfo /Preserve /UsePrologue false /ColorSettingsFile () /AlwaysEmbed [ true ] /NeverEmbed [ true ] /AntiAliasColorImages false /CropColorImages true /ColorImageMinResolution 300 /ColorImageMinResolutionPolicy /OK /DownsampleColorImages true /ColorImageDownsampleType /Bicubic /ColorImageResolution 300 /ColorImageDepth -1 /ColorImageMinDownsampleDepth 1 /ColorImageDownsampleThreshold 1.50000 /EncodeColorImages true /ColorImageFilter /DCTEncode /AutoFilterColorImages true /ColorImageAutoFilterStrategy /JPEG /ColorACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /ColorImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000ColorACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000ColorImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasGrayImages false /CropGrayImages true /GrayImageMinResolution 300 /GrayImageMinResolutionPolicy /OK /DownsampleGrayImages true /GrayImageDownsampleType /Bicubic /GrayImageResolution 300 /GrayImageDepth -1 /GrayImageMinDownsampleDepth 2 /GrayImageDownsampleThreshold 1.50000 /EncodeGrayImages true /GrayImageFilter /DCTEncode /AutoFilterGrayImages true /GrayImageAutoFilterStrategy /JPEG /GrayACSImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /GrayImageDict << /QFactor 0.15 /HSamples [1 1 1 1] /VSamples [1 1 1 1] >> /JPEG2000GrayACSImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /JPEG2000GrayImageDict << /TileWidth 256 /TileHeight 256 /Quality 30 >> /AntiAliasMonoImages false /CropMonoImages true /MonoImageMinResolution 1200 /MonoImageMinResolutionPolicy /OK /DownsampleMonoImages true /MonoImageDownsampleType /Bicubic /MonoImageResolution 1200 /MonoImageDepth -1 /MonoImageDownsampleThreshold 1.50000 /EncodeMonoImages true /MonoImageFilter /CCITTFaxEncode /MonoImageDict << /K -1 >> /AllowPSXObjects false /CheckCompliance [ /None ] /PDFX1aCheck false /PDFX3Check false /PDFXCompliantPDFOnly false /PDFXNoTrimBoxError true /PDFXTrimBoxToMediaBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXSetBleedBoxToMediaBox true /PDFXBleedBoxToTrimBoxOffset [ 0.00000 0.00000 0.00000 0.00000 ] /PDFXOutputIntentProfile () /PDFXOutputConditionIdentifier () /PDFXOutputCondition () /PDFXRegistryName () /PDFXTrapped /False /CreateJDFFile false /Description << /ARA /BGR /CHS /CHT /CZE /DAN /DEU /ESP /ETI /FRA /GRE /HEB /HRV (Za stvaranje Adobe PDF dokumenata najpogodnijih za visokokvalitetni ispis prije tiskanja koristite ove postavke. Stvoreni PDF dokumenti mogu se otvoriti Acrobat i Adobe Reader 5.0 i kasnijim verzijama.) /HUN /ITA /JPN /KOR /LTH /LVI /NLD (Gebruik deze instellingen om Adobe PDF-documenten te maken die zijn geoptimaliseerd voor prepress-afdrukken van hoge kwaliteit. De gemaakte PDF-documenten kunnen worden geopend met Acrobat en Adobe Reader 5.0 en hoger.) /NOR /POL /PTB /RUM /RUS /SKY /SLV /SUO /SVE /TUR /UKR /ENU (Use these settings to create Adobe PDF documents best suited for high-quality prepress printing. Created PDF documents can be opened with Acrobat and Adobe Reader 5.0 and later.) >> /Namespace [ (Adobe) (Common) (1.0) ] /OtherNamespaces [ << /AsReaderSpreads false /CropImagesToFrames true /ErrorControl /WarnAndContinue /FlattenerIgnoreSpreadOverrides false /IncludeGuidesGrids false /IncludeNonPrinting false /IncludeSlug false /Namespace [ (Adobe) (InDesign) (4.0) ] /OmitPlacedBitmaps false /OmitPlacedEPS false /OmitPlacedPDF false /SimulateOverprint /Legacy >> << /AddBleedMarks false /AddColorBars false /AddCropMarks false /AddPageInfo false /AddRegMarks false /ConvertColors /ConvertToCMYK /DestinationProfileName () /DestinationProfileSelector /DocumentCMYK /Downsample16BitImages true /FlattenerPreset << /PresetSelector /MediumResolution >> /FormElements false /GenerateStructure false /IncludeBookmarks false /IncludeHyperlinks false /IncludeInteractive false /IncludeLayers false /IncludeProfiles false /MultimediaHandling /UseObjectSettings /Namespace [ (Adobe) (CreativeSuite) (2.0) ] /PDFXOutputIntentProfileSelector /DocumentCMYK /PreserveEditing true /UntaggedCMYKHandling /LeaveUntagged /UntaggedRGBHandling /UseDocumentProfile /UseDocumentBleed false >> ] >> setdistillerparams << /HWResolution [2400 2400] /PageSize [612.000 792.000] >> setpagedevice