Modeling the processes of deformation and destruction of the rock sample during its extraction from great depths

Article investigates the change in the geophysical properties of rocks in the process of extracting the rock sample from great depths. Evaluation of changes in effective elastic properties, porosity and permeability of rock samples during extraction was carried out by means of finite element modeling. Assessment of the critical dimensions and orientation of internal defects, leading to the destruction of the rock samples during extraction from great depths, has been made based on the methods of linear destruction mechanics. Approach that makes it possible to calculate the change in the mechanical properties, porosity and fracturing of reservoir rocks in the process of extracting the rock sample from depths to the surface is proposed. Use of refined data on the mechanical properties of recoverable rock samples makes it possible to increase the accuracy of digital geological models required for geological exploration, determination of reservoir properties and oil and gas saturation of a field, and development of oil and gas deposits. Application of such models is especially relevant at all stages of the fields development with hard-to-recover reserves.

to clarify the physical and mechanical properties of rocks in reservoir conditions, which can be used to calibrate geophysical data.
Statement of the problem. Use of deformation models for rock sample, which directly take into account the presence of pores, allows direct assessment of the change in their volume when the rock sample is extracted to the surface. Digital rock sample model, taking into account its microstructure (Fig.1), is considered as an object of research. Construction of a three-dimensional digital rock sample model based on computed tomography (CT) data involves frame-by-frame digitization of CT data. It is done to transform a set of pixel-by-pixel 2D sections of the studied rock sample (Fig.1, b) into a digital 3D model with subsequent comparison of individual voxels in the 3D model with the developed finite element (FE) of rock sample model (Fig.1, c). This approach allows constructing a FE model reliably, taking into account the presence of micropores, to determine the effective mechanical characteristics of the studied oil and gas source rocks by the FE homogenization method (Fig.1).
The paper considers a rock sample simulation model, which assumes the presence of stochastically located isolated and interconnected pores (Fig.1, d). The resolution of the digital model is 750 m. It was assumed that connected pores are located along randomly oriented lines. Volume share of isolated and interconnected pores is 2 % of the total rock sample volume. Rock sample with the following geometric characteristics was considered ( Fig.1, a): radius R is 15 mm, height H is 80 mm.
Increasing requirements for the accuracy of modeling geomechanical processes lead to the need to use verified models for deformation and destruction of rocks. A wide range of approaches to modeling is now used: finite element method (FEM) [3], method of movable cellular automata [5,20]; approaches based on molecular dynamics [6,24] and percolation methods [29]. It is necessary to highlight the approaches that consider geomaterials as a homogeneous environment; directly, as well as indirectly, take into account the presence of pores (e.g., using the models of poroelasticity and poroplasticity [21]). FEM in the form of the displacement method is used in this study to simulate the processes of deformation and destruction of the rock sample during its extraction.
At deep bedding, the reservoir rock is in a non-uniform viscoelastic-plastic state. Reason for the occurrence of plastic (residual) deformations is the significant weight effect of the overlying rocks. When the rock sample is extracted, the load is removed. Thus, for the considered loading conditions, the solution of the boundary problem within the framework of the non-isothermal infinitesimal linear elastic formulation, which was used in the study, is admissible. State of the rock sample in the formation conditions was considered as a reference configuration. To determine the stress-strain state during extraction, the vertical and lateral tensile pressures were considered, as well as the change in temperature in the volume, which corresponds to the change in thermobaric conditions during extraction from great depths.
Process of unloading the rock sample during extraction has a two-stage nature. First, when drilling-over the rock sample, rock pressure is replaced by the hydrostatic pressure of the drilling fluid. This process ends at a distance of 1.0-1.5 rock sample diameter from the borehole bottom. Further, during extraction, the rock sample is in conditions of decreasing uniform all-round compression by a hydrostatic column of the drilling fluid. After that, the hydrostatic pressure is replaced by atmospheric pressure. Use of the superposition principle within the framework of the linear formulation of the problem makes it possible to disregard the sequence of loads application. Possibility of applying this approach is also confirmed by the experimental data obtained in [4]. Dependence of the change in open porosity on the magnitude of the rock pressure has a pronounced linear character. Methodology. Removal of compressive stresses during extraction was modeled by setting tensile vertical and lateral pressures; the temperature change was set negative (cooling). Initial vertical pressure was taken equal to 100 MPa, and lateral pressure was 30 MPa. Due to the linear formulation of the problem, the results for other initial thermobaric conditions can be obtained by linear scaling of the acquired results. For example, for Bazhenov-type field Bolshoy Salym [9], initial axial pressure is 50 MPa; radial pressure is 15 MPa.
During the calculations, it was assumed that both the skeletal material of the reservoir rock and the effective pore material were isotropic. In FE calculations, following parameters of materials were used [7]: skeleton/pores -Young's modulus 20/2 GPa; Poisson's ratio 0.22/0.49; thermal expansion coefficient is 5•10 -6 /2•10 -4 1/K.
Relative change in the pore volume р during extraction of the porous geomaterial to the surface is determined on the assumption of small deformations by the first invariant of the strain tensor : first invariant of the strain tensor averaged over the volume of all pores; first invariant of a tensor equal to the sum of its diagonal elements. Using the decomposition of the strain tensor into elastic  е and temperature  th components, the dependence is obtained: Data on the change in pore volume make it possible to judge the change in the permeability of rocks based on the Kozeny -Karman equation [13,14]: where kpermeability; porosity; fshape factor of a circular capillary cross-section; Sfspecific surface area of filter channels; Th -hydraulic tortuosity of channels. Relationship between changes in pore volume р and porosity  can be expressed as follows: Assessment of the parameters for internal defects (namely, dimensions and orientation) leading to the destruction of the rock sample during extraction is a relevant task that requires the use of computer modeling of the destruction process. Analysis of the conditions for the propagation of cracks was carried out based on methods of linear fracture mechanics by calculating the stress intensity factors (SIF) using the results of finite element modeling for the stress state of rock samples with cracks of various initial configurations.
Case of a rock sample with an idealized central disc-shaped crack was considered (Fig.2). Influence of the crack radius l and the angle between the plane of the crack and the horizontal plane , and the magnitude of the vertical rock pressure were investigated. Lateral pressure was set proportional to the vertical one with a coefficient of 0.3 in all cases of loading. Determination of the stress-strain state of a rock sample with a crack Fig was carried out within the boundaries of a linear-elastic formulation. In the considered model of a macrohomogeneous medium the temporal variation of a spatially homogeneous temperature field does not affect the growth of cracks. In this regard, the study did not take into account such a change.
Crack edges were considered free. SIF was determined based on the analysis for the distribution of displacement fields in the vicinity of the crack front.
First, the influence of the spatial orientation of the central inner disc-shaped crack in the case of its radius invariable was investigated. Change in the angle between the plane of the crack and the horizontal plane  in the range from 0 to 45 was considered. Crack radius l was taken equal to 3 mm, rock sample radius R -15 mm (l/R = 0.2).
All SIF (KI, KII and KIII) are non-zero. KII and KIII allow variation along the crack front. In the case of the simultaneous presence of several modes of destruction (a combination of normal detach, transverse and longitudinal shears), the effective SIF Keff was determined by the formula that corresponds to the criterion of the release rate of elastic deformation energy [28]: In the case of a horizontal crack () KII and KIII are equal to zero. Growth of the crack and, accordingly, destruction of the rock sample will begin under the condition Keff  KIС, where KIС is the critical SIF characterizing the crack resilience of the material.
Discussion. Determination of the change in pore volume and porosity was studied for two structural levels: (i) a representative volume with a separate isolated spherical pore in an idealized periodic porous structure (Fig.3) and (ii) rock sample directly, taking into account the presence of isolated and interconnected pores (see Fig.1, c and d). Calculations were performed using the FE software ANSYS. Twenty-node three-dimensional isoparametric FE with quadratic approximation inside the FE in the fields of nodal displacements were used in the calculations.
Initially, changes in soil porosity during extraction were estimated based on the analysis of the relative change in the volume of an isolated spherical pore in an idealized periodic structure (Fig.3).
Problem was solved in a three-dimensional statement. Due to symmetry, 1/8 of the elementary representative volume (ERV) of a porous material in the form of a cube with a central spherical pore was considered with the setting of symmetry conditions on the corresponding edges. Tensile pressure was applied to the upper edge corresponding to the removed vertical rock pressure р  . Pressure, corresponding to the removed lateral rock pressure р  , was applied to the side edges. Ratio of the ERV sizes and spherical pores provides a porosity of 4 %, which corresponds to the porosity of reservoir rocks in the Bazhenov formation [15,17,18].
Contribution to the change in pores by individual factors was determined by means of FE modeling of the rise: changes in pressure and temperature. Thus, the removal of pressure leads to an increase in pore volume by 0.85 % in relation to the initial volume, a decrease in temperature by 100 С leads to a decrease in pore volume by 0.2 %, which does not contradict [18]. Figure 4 shows the fields characterizing the stressed and deformed states of an elementary representative volume of a porous rock sample material containing a single pore during a rise characterized by a change in where σij, εij, i, j = 1.3 ̅̅̅̅stress and strain tensor components; lateral deformation coefficient. Direct FE modeling of rock sample deformation processes, taking into account its microstructure, was performed to assess the change in the volume of oil-saturated pores of reservoir rocks during their extraction to the surface. Problem was solved with the following boundary conditions: a tensile pressure was applied to the upper plane of the rock sample, corresponding to the removed vertical rock pressure; pressure, corresponding to the removed side rock pressure, was applied to the side edge (of the cylindrical surface). Zero equality of axial displacements was set on the lower edge. Center of the lower section of the rock sample was fixed in two directions in the plane and fixing one more degree of freedom in the direction along the generatrix for a point on the outer radius to exclude solid (translational and rotational) movements. Model contains 328 704 FE. Investigations of ERV deformation showed that temperature change does not significantly affect the change in pore volume and in this formulation it was not taken into account.  Removal of rock pressure during extraсtion leads to an increase in pore volume by 1 %. Fields characterizing the stress-strain state of the rock sample and separately for the pores are shown in Fig.5.
Use of refined data on the mechanical properties of recovered rock samples improves the accuracy of digital geological models. One of the promising methods for determining effective mechanical characteristics is the FE homogenization method, which ensures the equivalence of the energies of heterogeneous and homogeneous media, carried out using data on the microporous structure of rocks obtained using CT.
Effective elastic properties of the homogenized rock sample material are determined based on the equality (generalized Hooke's law): where C 4 effective elastic compliance tensor of the 4th order. Tensors of deformations and stresses, averaged over the volume of the representative element, are determined by the following equalities: The line above the introduced tensors means the correspondence to the homogenized (averaged) material.
Elastic characteristics of the rock sample correspond to the anisotropic material. Model of a transversely isotropic material with an anisotropy axis along the direction of vertical weight effect can be considered as an adequate approximation. As a first approximation, on the assumption that there is no influence of the weight effect on the formation of pores and the isotropy of the extrapore skeleton properties, rock sample material can be considered as isotropic. In this case, as a consequence of equation (7) under uniaxial tension for the axial components of stresses and strains, the simplest relation is valid: 11 11 Effective Young's modulus of the material was calculated for the investigated rock sample (Fig.1): For the considered porosity value, the Young's modulus value determined using FE homogenization practically coincides (the difference is less than 0.5 %) with the value determined based on the rule of mixtures: where E p and E k -elastic moduli for effective material of pores and skeleton. However, application of this rule is only possible with porosity values less than 10 %. At large values, an increase in the range of possible changes in the elastic properties of the heterogeneous rock sample material (Hill's fork [26]) leads to significant errors in determining the constants.
Proposed approach can also be applied when lowering the resolution of the digital rock sample model. In this case, it is required to determine the effective properties of enlarged voxels based on a model with a high resolution, which can be performed using FE homogenization.
Results obtained indicate a significant degree of inhomogeneity for the stress and strain states of the rock sample caused by changes in pressure and temperature. Temperature has no significant effect on the change in pore volume. Comparison of the results showed that under equivalent thermobaric conditions, the increase in pore volume determined on the basis of a full-size rock sample with a microporous structure is 15 % greater than that determined on the basis of ERV rock sample with a single pore. This is due to the fact that the full-size model takes into account the stochastic (nonperiodic) arrangement of pores, the presence of related porosity, and the mutual influence of pores on each other. Based on the results obtained, the elastic properties of the rock sample material can be estimated and the effective elastic properties of individual components of the model can be recalculated when its dimension is reduced.
Rock sample material was assumed to be macrohomogeneous and isotropic during solving the problems regarding a crack. Determination of the effective elastic properties of the rock sample was carried out by the finite element homogenization method based on the developed digital rock sample model.
Dependencies of KI and KII on the angle between the plane of the crack and the horizontal plane  are shown in Fig.6. At the top and bottom points of the crack front (point of the maximum Keff realization), KIII value is zero. In this case, with an increase in the inclination angle of the crack, KI decreases monotonically, and KII monotonically increases.  Dependence of the effective SIF Keff on the deviation angle of the crack plane from the horizontal is shown in Fig.7, a. KIС coefficients for reservoir rocks vary with in 0.5-1.7 MPa·m 1/2 [11]. In this study, the SIF critical value KIС was taken equal to 1.46 MPa·m 1/2 , which corresponds to the value of the crack resilience of the rock sample material from the Akkum gas-condensate field (Republic of Kazakhstan). Comparing the obtained dependences Keff with KIС, it can be concluded that for a crack with a radius of 3 mm at an axial pressure of 50 MPa, destruction occurs at any orientation of the crack within the considered range. At a pressure of 30 MPa, destruction occurs at any inclination angles of the crack less than 35°. At a pressure of less than 20 MPa, destruction does not occur for cracks with any orientation. Thus, for each orientation of the crack and its radius, the upper limit of the axial pressure (depth of bedding) can be indicated, below which the destruction of rock sample failure does not occur.
At the next stage, the influence of the relative crack size l/R on SIF was investigated in the case of a horizontally oriented crack (  . Change in the crack radius in the range from 3 to 12 mm was considered. For small values of l/R, under the assumption of the lateral pressure absence, the results showed good agreement with the known analytical solution for an inclined crack in an unlimited space [28]. Difference does not exceed 0.5 %. Dependence of KI on the relative crack size is shown in Fig.7, b. Its analysis allows concluding that at axial pressure of 50 MPa, destruction occurs when the ratio l/R is more than 0.05 (l > 0.75 mm), at a pressure of 30 MPamore than 0.11, at 20 MPamore than 0.29.
Detection of the destruction fact of the rock sample during extraction from great depths carries information about the excess of the defect size for the critical values, the assessment of which can be made on the basis of nomograms similar to that shown in Fig.7 Obtained results make it possible to evaluate qualitatively and quantitatively the influence of various factors (size and orientation of internal defects, mechanical properties of the material, rocks depth) on the crack resilience of rock samples during extraction from great depths to the surface. In order to refine the results, the methods of non-linear destruction mechanics [2,22] and more complex material models should be used [21,23,30].
Conclusion. An approach that allows assessing the change in porosity and fracturing of reservoir rocks during rock sample extraction to the surface has been proposed. Change in the volume of oilsaturated rock samples pores during its extraction from a depth of 2900 m, corresponding to Bazhenov-type fields, was estimated by means of direct finite element modeling of core deformation processes. In order to determine the critical dimensions and orientation of internal defects leading to the destruction of the rock sample during extraction, the results of modeling the destruction process of rock samples with cracks of various initial configurations were obtained.
Dependence of the change in the volume of oil-saturated pores in the rock sample on the change in temperature and pressure conditions during its extraction from great depths was obtained based on the results of the finite element analysis of the change in the stress-strain state of the rock sample.
Calculations took into account in detail the microstructural heterogeneity of the rock sample properties based on the use of the developed digital rock sample model. Method of finite element homogenization was used to determine the effective properties.
Results of the influence of the size and spatial orientation for the internal circular crack in the rock sample on its destruction during extraction, found based on the results of finite element calculations using the criteria of linear destruction mechanics, are obtained and analyzed. Results of the influence of axial and lateral pressures on the critical values of crack size and deviation angle of the crack plane normal from the rock sample axis are presented.