WO2011155862A1 - Способ определения напряженно-деформированного состояния слоистой среды - Google Patents

Способ определения напряженно-деформированного состояния слоистой среды Download PDF

Info

Publication number
WO2011155862A1
WO2011155862A1 PCT/RU2010/000299 RU2010000299W WO2011155862A1 WO 2011155862 A1 WO2011155862 A1 WO 2011155862A1 RU 2010000299 W RU2010000299 W RU 2010000299W WO 2011155862 A1 WO2011155862 A1 WO 2011155862A1
Authority
WO
WIPO (PCT)
Prior art keywords
layers
equations
layer
stress
parameters
Prior art date
Application number
PCT/RU2010/000299
Other languages
English (en)
French (fr)
Inventor
Максим Андреевич ЧЕРТОВ
Original Assignee
Щлюмберже Холдингс Лимитед
Шлюмберже Текнолоджи Б.В.
Шлюмберже Канада Лимитед
Сервисес Петролиерс Шлюмберже
Прад Рисеч Энд Девелопмент Лимитед
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Щлюмберже Холдингс Лимитед, Шлюмберже Текнолоджи Б.В., Шлюмберже Канада Лимитед, Сервисес Петролиерс Шлюмберже, Прад Рисеч Энд Девелопмент Лимитед filed Critical Щлюмберже Холдингс Лимитед
Priority to PCT/RU2010/000299 priority Critical patent/WO2011155862A1/ru
Priority to US13/701,447 priority patent/US20130151160A1/en
Publication of WO2011155862A1 publication Critical patent/WO2011155862A1/ru

Links

Classifications

    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V9/00—Prospecting or detecting by methods not provided for in groups G01V1/00 - G01V8/00
    • G—PHYSICS
    • G01—MEASURING; TESTING
    • G01M—TESTING STATIC OR DYNAMIC BALANCE OF MACHINES OR STRUCTURES; TESTING OF STRUCTURES OR APPARATUS, NOT OTHERWISE PROVIDED FOR
    • G01M5/00—Investigating the elasticity of structures, e.g. deflection of bridges or air-craft wings
    • G01M5/0066—Investigating the elasticity of structures, e.g. deflection of bridges or air-craft wings by exciting or detecting vibration or acceleration
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]

Definitions

  • the present invention relates to methods for determining the stress-strain state of layered media, with at least one porous-elastic layer filled with a liquid or gas, when changing the pressure of this liquid or gas.
  • the invention is applicable to a wide variety of such environments, which, in particular, include composite materials, geological formations, and bone tissues.
  • Poreelastic are layers in which their component medium is porous (the presence of pores, discontinuities) and elastic.
  • filling a poroelastic layer involves filling with a liquid or gas any such discontinuity.
  • the regions filled with liquid can be pores, cracks in the layer, voids of a cavity, etc.
  • the liquid can be, in particular, water, oil, etc. Natural gas, air, or any other gaseous substance can be considered as a gas.
  • An application of the method may be its implementation in the oil and gas industry, in particular, when determining the stress-strain state of a layered geological environment, characterized by the distribution of stresses and displacements when the pore pressure changes during the development of oil and gas fields (injection, production wells).
  • the present invention was the development of a method for determining the elastic stress-strain state of a layered medium, providing high accuracy of the results at a high speed of their receipt.
  • the method developed in the framework of the present invention can be used for, for example, preliminary monitoring of development projects to obtain preliminary estimates of the risk of large deformations at an early stage of development, as well as to decide on the need for full-scale 3D modeling.
  • the method can also be applied to conduct quick studies on sensitivity to changes in parameters and studies on adaptation by parameters, or when solving the inverse problem to determine the state of the inner zones of the collector based on available measurements of displacements and stresses.
  • a quick assessment of possible changes in stresses and displacements in the volume of the field allows planning necessary investments at an early stage of field development.
  • the well-known semi-analytical model used in the oil industry is the method of point sources of deformation, see D. Girtsma, Sagging of soil over compressible oil and gas reservoirs, Journal of Oil Technology, Volume 25, p. 7340744, 1973. (J. Geertsma Land subsidence over compacting oil and gas reservoirs. J. Pet. Tech. Vol. 25, pp. 734-744, 1973).
  • the depleted three-dimensional region the collector is divided into many small-sized sources of deformation of a cubic shape.
  • the contribution to the subsidence of the surface from each source is calculated based on the analytical expression for a single point source of deformation immersed in a semi-infinite elastic half-space.
  • Each such elementary solution is applicable at a large distance compared to the source size and, therefore, the method is not intended to calculate the deformation inside the depleted region itself. In addition, the method cannot be applied when the reservoir and surrounding rocks have different elastic properties.
  • the above method is also not applicable for calculating the displacement inside the depleted region, and, at the same time, requires additional calibration and adjustment in order to select the number, coordinates and magnitude of imaginary sources, which allow obtaining an acceptable residual displacement at the boundaries of the layers.
  • the known solutions have limited technical capabilities and, either, do not allow to determine the stress-strain state of the geological environment with sufficient completeness and accuracy, or, do not provide results with the required speed.
  • the technical result achieved by the claimed method can be attributed to the speed of the method when obtaining high accuracy of the desired values.
  • the claimed technical result is achieved in that: a) determine the initial parameters of the layered medium, b) determine the space-distributed load of the developed layers of the layered environment, c) compose a system of equations, including: equations of elastic equilibrium for each layer, equations that define the boundary conditions at the boundaries of the layers: continuity of stress between the layers, zero stress at the free boundary, zero displacement at infinity, d) obtain analytical solutions of the equations of elastic equilibrium in each layer, provided that the layers are plane-parallel, homogeneous and poroelastic, e) coordinate the analytical solutions obtained in different layers, with each other, taking into account the corresponding the existing boundary conditions, e) determine the values of stress and displacement at selected points of the layered medium in accordance with the required approximation accuracy based on agreed analytical solutions.
  • Initial parameters include geological parameters of the layers.
  • the geological parameters of the layers include: Lame elastic coefficients ( ⁇ , ⁇ ), porous elastic Biot parameter, layer thickness.
  • Geological parameters are measured using measuring instruments and / or are established empirically and / or are chosen arbitrarily, based on the requirements of the problem being solved.
  • the volume-distributed load is determined through the pore pressure gradient of the developed layers.
  • Pore pressure is calculated by solving the piezoelectric conductivity equation by the semi-analytical hydrodynamic method.
  • Steps d) and e) are implemented in Cartesian coordinates or in cylindrical coordinates.
  • the initial parameters of the layered medium are determined, which include, in particular, geological parameters characterizing the properties and condition of the geological environment of the field.
  • a layered medium is a vertical sequence of horizontal plane-parallel layers, where each layer is homogeneous and poroelastic.
  • SGM Layered Geomechanical Model
  • Plant parallel are considered layers bounded on two sides by parallel horizontal planes.
  • Porous are layers in which the medium, its component is porous (for example, contains interconnected pores distributed uniformly) and elastic (the medium is considered to be continuous, namely, continuously distributed in space and having elastic properties).
  • a medium is characterized by: elastic moduli ⁇ , ⁇ (Lame coefficients), porous elastic Bio parameter a and thickness of each layer (d).
  • Such parameters within the framework of the invention form geological parameters.
  • FIG. 2 A schematic sectional view of a layered three-dimensional field is shown in FIG. 2.
  • geological parameters is carried out by direct measurements using known means and methods used for such measurements in the oil and gas industry.
  • Such tools include, for example, sensors distributed over the respective sections of the section, seismic and acoustic measuring instruments, tiltmeters, and other downhole measuring equipment. So, for example, the thicknesses of the layers d are defined as the distances between the lithological boundaries, visible as jumps in the values on the log curves of various nature; elastic moduli ⁇ , ⁇ can be measured based on the analysis of acoustic logging; the bioelastic parameter Bio a can be estimated using laboratory mechanical tests of the extracted core.
  • geological parameters can also be established on the basis of existing knowledge, research obtained earlier, i.e. to choose arbitrarily, based on the requirements of the problem being solved, and also to determine empirically or using expert estimates.
  • the next step in the implementation of the method is to build SGM deposits, namely, a complete solution for a three-dimensional layered medium.
  • SGM deposits namely, a complete solution for a three-dimensional layered medium.
  • Such a model is described by a system of equations characterizing the elastic stress-strain state of equilibrium for each layer, i.e. a set of internal stresses and strains arising under the action of external loads, as well as equations defining the boundary conditions at the boundaries of the layers (between layers, as well as at the upper and lower boundaries of the layers).
  • the upper surface of the uppermost layer is free from stress.
  • the last, lowest layer is semi-infinite in depth with zero displacements at infinity.
  • the boundary conditions in the lateral directions are also zero displacements at infinity.
  • the boundary conditions between the layers are the continuity of the stress components acting on the elementary site of the interface (interface), and the continuity of the vertical displacements.
  • the tangential displacements between the layers are proportional to the normal component of stresses acting in the same direction, multiplied by the compliance coefficients specified during the construction of the GMS. By default, these compliance factors are equal to zero, which corresponds to rigidly interlocked layers.
  • Each individual poroelastic layer can have an independently specified profile of changes in the volume-distributed load, which is determined by the pressure of the fluid inside it, p (x, y, z).
  • This pressure profile is represented as the product of the horizontal pressure component, which varies as an arbitrary function of x and y, and the vertical component, which is a third-order polynomial in z.
  • the displacements and stresses inside the layered half-space are calculated using the two-dimensional Fourier transform by solving the corresponding system of algebraic equations. Ultimately, the complete solution, if necessary, is converted back into direct space in order to obtain the distribution of displacements and stresses. It is also possible to visualize the obtained values.
  • the accuracy of the approximation in our case is determined by how many discrete points are used to approximate the continuous functions that we obtain as input and which we calculate using the proposed method.
  • the information on the distribution of pore pressure necessary for predicting deformation can be obtained from external sources, i.e. calculated by any of the methods known for such a purpose, or established empirically.
  • the hydrodynamic parameters of the layers include: the initial pore pressure of the layer, the permeability of the layer (k), the porosity of the layer, the initial gas or oil saturation, the viscosity of the oil, the viscosity of the displacement agent and the displaced fluid, the initial and final saturation of the displacement agent.
  • the geological and technological parameters of the wells include: the coordinates of the well intersection for each formation, the depth of the roof and the bottom of each formation in the context of the well, the flow rate of the well, the dynamics of the current and cumulative oil or gas production, the dynamics of the current and accumulated injection of the displacement agent.
  • the task parameters simulated in GREAT can be set using either input data files in a text format using keywords, or using the application programming interface (API, Application Programming Interface) (GREAT, designed as a dynamically connected library (dll).
  • API Application Programming Interface
  • GREAT a dynamically connected library
  • Data exchange via the API is significantly faster, data exchange through text files is more transparent for debugging and requires less modifications when switching to newer versions of GREAT.
  • the calculation process in GREAT is based on summing the contributions from elementary solutions for point sources and flat boundaries and does not imply approximation on a discrete grid covering the entire computational domain in order to reduce the effective dimension of the problem and achieve a high calculation speed.
  • Input data for the associated SGM + GREAT calculation is organized as a text file consisting of a set of sections.
  • An example input file is provided in the description section of the “Preferred Embodiment”.
  • the FILES section defines the full path to the GREAT executables and the name of the output file.
  • the DIMENSIONS section describes the geometry and dimensions of the numerical problem in horizontal section.
  • the section defines the grid dimension of pressure sensor wells, the distance between them and the grid size in the SGM.
  • the porous collector modeled by GREAT must be immersed inside a larger elastic region, since the boundary conditions for the mechanical problem must be set at a sufficiently large distance from the sources of applied load in the model.
  • Distance from collector to the outer lateral boundary is specified in terms of factors Bx, By defining the relationship between the distance from the boundary of the collector to the outer boundary and the size of the collector, as shown in FIG. 4.
  • the boundary conditions at the external boundaries in the SGM are zero displacements.
  • the grid spacing between the approximation nodes in the SGM is chosen equal to the distance between the pressure sensors in GREAT. In total, the SGM calculation has (2Bx + l) - (2By + l) -MN nodes, where MN is the dimension of the grid of sensor wells.
  • the LAYERS section describes the sequence of layers, starting from the free surface in the direction of increasing depth, layer thickness and mechanical properties, as shown in FIG. 2.
  • the layers in which development is underway, the pressure changes and where it is necessary to read the pressure field calculated outside the SGM, are marked with the flag "PRESSURE_LOAD l".
  • two options for loading external data are possible: 1) a text file containing a table of MxN pressure values; and 2) a link to the input file of the GREAT task for this layer, which needs to be calculated.
  • the input file is expanded by adding to its end an automatically generated grid of MxN vertical closed pressure sensor wells, perforated sections of such wells are placed in the middle of the layers being developed.
  • the "OUTPUT" section defines a list of z coordinates of the layers in which you want to calculate the output.
  • the main steps of the algorithm are presented in FIG. 5.
  • Analysis of the results obtained after performing one SGM + GREAT simulation includes verification of the results, which can be a comparison with others reference solutions, either in the whole or in some part of the domain of definition of the model, or it can be limited only to verifying that the solution falls within a certain acceptable range of solutions. If the solution is not satisfactory, then the whole process returns to the initial step in FIG. 1 to correct the input parameters defining the model. Then, the simulation process shown in FIG. 1, repeated. A similar repetition of several runs of modeling with correction of input data is also carried out when setting tasks for adaptation by parameters. In this case, the resulting solution is compared with a set of target parameters selected by the operator, and the variation of the input parameters continues until the correspondence to the target values is obtained within a certain accuracy. Conducting a large number of SGM + GREAT simulations is also necessary when conducting parametric studies to determine how varying input parameters affects the final solution.
  • V is the differential operator, and is the displacement vector of the geological medium at the selected points of the layer
  • ⁇ , ⁇ are the Lame elastic coefficients
  • f is the volume-distributed load.
  • Equation (1) is written in invariant form. In what follows, (1) will be solved in Cartesian coordinates, in which the z axis is oriented perpendicular to the plane of the layers. The solution is built with taking into account the following assumptions: all externally applied load arises only due to the pore pressure gradient, where a
  • Equation (1) can be converted to the form:
  • the superscript for the variables is ⁇ and denotes the layer number; the numbering of the layers i— ⁇ ... n goes from top to bottom (see Fig. 2).
  • the internal boundaries are ideal elastic contacts between the layers and ideally impermeable barriers to fluid flow. Therefore, on the interfaces (surfaces) there is no jump in the following voltage components:
  • C is the coefficient of compliance of the interface layer in the tangent direction, isotropic in the plane (x, y).
  • Equations (27-32) express the boundary conditions at the internal interfaces between the layers as i— ⁇ ... n- ⁇ .
  • the Fourier transform is the most expensive part of the algorithm, because it requires (X V-log MAO operations for each layer with a non-zero pressure distribution, for each calculated component ( ⁇ #) and each unique ⁇ coordinate in which the output data are calculated.
  • Inverse matrix of boundary conditions requires only 0 (L MN) operations, where L is the number of layers.
  • FFTW library FFTW library
  • U.S. CUFFT Library NVIDIA Corp., 2008, http://www.nvidia.com/object/cuda_develop.html. It is also possible to use any other existing or developed FFT libraries that best suit the needs of the user.
  • FIG. 2 Schematic representation of a section of a layered three-dimensional field.
  • FIG. 3 The main stages of the calculation algorithm for the SGM.
  • FIG. 4 Schematic representation of a rectangular reservoir with a grid of pressure transducer wells immersed in surrounding elastic rocks.
  • FIG. 5 - The main stages of the algorithm associated SGM + GREAT calculation.
  • FIG. 6 Illustration of a simulated problem with a five-point development scheme in a two-layer reservoir.
  • the rectangular collector has a size of 10,000 by 10,000 feet, it is surrounded on all sides by impermeable elastic layers, which are necessary for the correct definition of boundary conditions at a great distance from the loaded central part of the problem. In this case, zero stresses are set on the upper free surface, and zero normal displacements and slip relative to tangential displacements are set on the lateral and lower boundaries.
  • Both the layers under development are crossed by four producing wells, designated P-1 - P-4, and one injection, indicated by 'INJ'.
  • the geological parameters of the layers, as well as their permeability, and the Biot parameter are shown in Table 1.
  • the determination of the external load applied to the layered geological environment of the field as a result of its development using a five-point scheme in the example under consideration is determined by calculating the pressure with the GREAT software package.
  • a constant fluid flow is set equal to 800 barrels / day at the producing wells and 1,600 barrels / day at the injection well.
  • Development modeling is carried out for 1440 days, the change in displacements and stresses is calculated at the end of this period.
  • the pressure distribution in the two-layer manifold is calculated using two independent GREAT launches individually for each layer.
  • flows through the intervals of wells crossing the formations are distributed in proportion to the permeabilities and thicknesses of the layers in the proportion of 1/4 for the upper layer and 3/4 for the lower.
  • the pressure distribution is calculated on a uniform grid of 50 to 50 fictitious wells with zero flow in the central part of the problem. Then, the calculated pressure field is converted into a format that can be loaded into the SGM and used as a mechanical load applied to the layers.
  • the full set of input parameters is duplicated in the example input file for GREAT (see below).
  • the given initial parameters and the load distribution are substituted into the compiled system of equations from the equations of elastic equilibrium and boundary conditions.
  • Analytical solutions for each layer were obtained in the Fourier space and written as an integral transform, therefore, a fast Fourier transform is performed to obtain solutions for individual layers at given values of their parameters. After this, decisions are made in each layer in order to achieve the fulfillment of boundary conditions between the layers. For this, the matrix of boundary condition equations is filled in, which is then inverted using the exact non-iterative matrix inversion algorithm. As a result of this treatment, the free coefficients of the analytical solutions recorded for each layer are determined.
  • FIG. 7 and FIG. Figure 8 shows examples of calculated vertical displacements and the normal components of the stress tensor in the direction of one of the horizontal coordinate axes, which were obtained with the above values of the input parameters.
  • the indicated displacement and stress components reach their maximum values near operating wells. So, the vertical displacements and 2 are about -9.5 cm near the producing wells and about 8 cm at the injection well.
  • the stress tensor component is of the order of -4.2 MPa for production wells and about 7.8 MPa for an injection well.
  • the input file for the lower development layer “5spot_GR_btm.txt” is similar.

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Aviation & Aerospace Engineering (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Software Systems (AREA)
  • Mathematical Physics (AREA)
  • Databases & Information Systems (AREA)
  • Data Mining & Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Life Sciences & Earth Sciences (AREA)
  • Geophysics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

Настоящее изобретение относится к способам определения напряженно-деформированного состояния слоистых сред, с, по крайней мере, одним пороупругим слоем, заполненным жидкостью или газом, при изменении давления этой жидкости или газа и может быть применимо, в частности, в нефтегазовой отрасли. При реализации способа осуществляют: а) определение исходных параметров слоистой среды, б) определение объемно-распределенной нагрузки разрабатываемых слоев слоистой среды, в) составление системы уравнений, включающей: уравнения упругого равновесия для каждого слоя, уравнения, определяющие граничные условия на границах слоев: непрерывность напряжений между слоями, нулевые напряжения на свободной границе, нулевые смещения на бесконечности, г) получение аналитических решений уравнений упругого равновесия в каждом слое, при условии, что слои являются плоскопараллельными, однородными и пороупругими, д) согласование между собой аналитических решений, полученных в разных слоях, с учетом соответствующих граничных условий, е) определение значений напряжения и смещения в выбранных точках среды в соответствии с требуемой точностью аппроксимации, на основе согласованных аналитических решений. Техническим результатом является быстродействие способа при получении высокой точности искомых значений.

Description

Способ определения напряженно-деформированного состояния слоистой среды
Область техники
Настоящее изобретение относится к способам определения напряженно-деформированного состояния слоистых сред, с, по крайней мере, одним пороупругим слоем, заполненным жидкостью или газом, при изменении давления этой жидкости или газа.
Изобретение применимо для широкого разнообразия таких сред, к которым, в частности, могут быть отнесены, композитные материалы, геологические формации, костные ткани.
«Пороупругими» являются слои, в которых среда их составляющая является пористой (наличие пор, несплошностей) и упругой.
Так, заполнение пороупругого слоя подразумевает заполнение жидкостью или газом любой такой несплошности. В частности, заполненными жидкостью областями могут быть поры, трещины в слое, пустоты каверны и пр. Жидкостью может являться, в частности, вода, нефть и др. В качестве газа может рассматриваться природный газ, воздух, или любое другое газообразное вещество.
Вариантом применения способа может быть его реализация в нефтегазовой отрасли, в частности, при определении напряженно- деформированного состояния слоистой геологической среды, характеризуемого распределением напряжений и смещений при изменении порового давления при разработке нефтяных и газовых месторождений (нагнетающие, добывающие скважины). Задачей настоящего изобретения явилась разработка способа по определению упругого напряженно-деформированного состояния слоистой среды, обеспечивающего высокую точность полученных результатов при высокой скорости их получения.
Предшествующий уровень техники
Так, в нефтегазовой отрасли неконтролируемые деформации и изменения напряжений могут приводить к многочисленным проблемам, включая повреждение и разрушение обсадных колонн, обрушение скважин, вынос песка, ухудшение проницаемости и пр. Таким образом, при разработке месторождений крайне важно осуществлять постоянный мониторинг деформации горных пород и изменений напряжений со временем.
Разработанный в рамках данного изобретения способ может быть использован в целях, например, предварительного мониторинга проектов разработки для получения предварительных оценок риска возникновения больших деформаций на ранней стадии разработки, а также для принятия решения о необходимости проведения полномасштабного 3D моделирования. Способ может быть также применен для проведения быстрых исследований на чувствительность к изменению параметров и исследований на адаптацию по параметрам, или при решении обратной задачи для определения состояния внутренних зон коллектора на основе имеющихся измерений смещений и напряжений. Наряду с другими преимуществами, проведение быстрой оценки возможных изменений напряжений и смещений в объеме месторождения позволяет осуществлять планирование необходимых инвестиции на ранней стадии разработки месторождений.
Из уровня техники известны методы по моделированию различных геомеханических процессов, связанных с изменением напряжений, деформаций, давлений и проницаемостей, полученных с использованием конечно-разностных и конечно-элементных механических и гидродинамических кодов, как коммерческого, так и исследовательского характера, в частности:
- П.Лонгмар, М.Мейнгай, П.Лемоньер, А.Онайси, Ч.Герард, Н. Кутсабелулис, Геомеханика при моделировании месторождений: обзор взаимосвязей и анализ проблем, Нефтегазовая наука и технология, Институт физики нефти, том 57, N°5, с.471-483, 2002 (P. Longuemare, М. Mainguy, P. Lemonnier, A. Onaisi, Ch. Gerard, N. Koutsabeloulis, Geomechanics in reservoir stimulation: overview of coupling methods and field case study, Oil & gas science and technology - Rev. IFP. Vol. 57 No. 5, pp. 471-483, 2002),
А.Сеттари, Ф.М.Моритц, Система для связанного гидродинамического и геомеханического моделирования месторождений, журнал SPE, с.219-226, 1998 (A. Settari, F.M. Mourits. А Coupled Reservoir and Geomechanical Simulation System, SPE Journal, pp. 219-226, 1998),
- Л.К.Томас, Л.Ю.Чин, Р.Г.Пирсон, Д.Е.Силт, L.K. Thomas, L.Y. Chin, R.G. Pierson, J.E. Sylte, Связанное геомеханическое и пластовое моделирование, журнал SPE, SPE-87339-PA, 2003 (Coupled geomechanics and reservoir simulation, SPE-87339-PA, 2003),
- П.Самьер, С.Дженнаро, Практическое итерационное связывание геомеханического и пластового и моделирования, журнал SPE, SPE- 106188, , 2007 (P. Samier, S. De Gennaro, Practical iterative coupling of geomechanics with reservoir simulation. SPE- 106188, 2007),
- патент US 7177764 B2, кл. GO 1 VI 1/00, опубл. 07.10.2004.
Данные методы основаны на аппроксимации моделируемой задачи и использовании дискретной сетки, что предоставляет значительную долю гибкости для моделирования различных конфигураций месторождения и последовательного улучшения качества аппроксимации.
Однако, такие методы характеризуются повышенной трудоемкостью и увеличенными временными затратами как на проведение вычислений, так и на подготовку и построение численной модели.
Известно так же использование полуаналитических моделей. Как правило, математические выкладки в этих моделях проведены для оценки только проседания геологической среды, вызванного откачкой жидкости из коллектора, то есть только вертикальной компоненты смещений, рассчитанной на свободной поверхности. Другие компоненты вектора смещений и тензора напряжений чаще всего не выводятся, хотя их получение в рамках большинства этих моделей вполне возможно.
Так, широко известной полуаналитической моделью, применяемой в нефтяной промышленности, является метод точечных источников деформации, см. Д.Гиртсма, Проседание почвы над сжимающимися нефтяными и газовыми пластами, Журнал Нефтяная Технология, том 25, с.7340744, 1973. (J. Geertsma, Land subsidence over compacting oil and gas reservoirs. J. Pet. Tech. Vol. 25, pp. 734-744, 1973). Согласно указанному методу истощенная трехмерная область коллектора разбивается на множество малых по размеру источников деформации кубической формы. Вклад в проседание поверхности от каждого источника вычисляется на основе аналитического выражения для единичного точечного источника деформации, погруженного в полубесконечное упругое полупространство. Каждое такое элементарное решение применимо на большом, по сравнению с размером источника, расстоянии и, следовательно, метод не предназначен для расчета деформации внутри самой истощенной области. Кроме этого, метод невозможно применить, когда коллектор и окружающие породы имеют различные упругие свойства.
Известен способ расчета проседания по, так-называемой, формуле Бруно, см. М.С.Бруно, Геомеханический анализ и анализ решений для предотвращения повреждения обсадных колонн, связанного со сжатием коллекторов, Журнал SPE, SPE 71695, 2001 (M.S. Bruno, Geomechanical analysis and decision analysis for mitigating compaction related casing damage. SPE 71695, 2001). Данная формула позволяет вычислить максимальное проседание над истощенной областью в форме диска.
Использование функции источника с кольцевым распределением нагрузки для расчета напряжений вокруг осесимметричных истощенных областей в однородном упругом полупространстве известно из см. П.Сегал, Напряжения вызванные откачкой жидкости из осесимметричных резервуаров, Теоретическая и Прикладная Геофизика, том 139 , 3/4, с.535-560, 1992 (Segall P., Induced stresses due to fluid extraction from axisymmetric reservoirs. Pure Appl. Geophys. 139, 3/4, pp. 535-560, 1992). Использование мнимых точечных источников деформации различных типов (точечные, дипольные, сферические источники), компенсирующих невязку смещений на границах между слоями с различными упругими модулями, известно из см. П.А.Фоккер, Б.Орлич, Полуаналитическое моделирование процесса проседания породы, Математическая еология, том 38, N°5, с.565-589, 2006 (P. A. Fokker, В. Orlie, Semi-analytic modeling of subsidence. Mathematical Geology, Vol. 38, No. 5, pp. 565-589, 2006).
Однако, вышеуказанный метод также не применим для вычисления смещения внутри истощенной области, и, вместе с тем, требует дополнительной калибровки и настройки, чтобы подобрать количество, координаты и величину мнимых источников, позволяющих получить приемлемую невязку смещений на границах слоев.
Таким образом, известные решения обладают ограниченными техническими возможностями и, либо, не позволяют определить напряженно-деформированное состояние геологической среды с достаточной полнотой и точностью, либо, не обеспечивают получение результатов с требуемым быстродействием.
Раскрытие изобретения
К техническому результату, достигаемому заявленным способом, можно отнести быстродействие способа при получении высокой точности искомых значений.
Заявленный технический результат достигается тем, что: а) определяют исходные параметры слоистой среды, б) определяют объемно-распределенную нагрузку разрабатываемых слоев слоистой среды, в) составляют систему уравнений, включающую: уравнения упругого равновесия для каждого слоя, уравнения, определяющие граничные условия на границах слоев: непрерывность напряжений между слоями, нулевые напряжения на свободной границе, нулевые смещения на бесконечности, г) получают аналитические решения уравнений упругого равновесия в каждом слое, при условии, что слои являются плоскопараллельными, однородными и пороупругими, д) осуществляют согласование аналитических решений, полученных в разных слоях, между собой с учетом соответствующих граничных условий, е) определяют значения напряжения и смещения в выбранных точках слоистой среды в соответствии с требуемой точностью аппроксимации на основе согласованных аналитических решений.
Возможны и другие варианты осуществления патентуемого способа, согласно которым:
Исходные параметры включают геологические параметры слоев.
Геологические параметры слоев включают: упругие коэффициенты Ламэ (λ, μ), пороупругий параметр Био, толщину слоя.
Геологические параметры измеряют с помощью средств измерения и/или устанавливают эмпирическим путем и/или выбирают произвольным образом, исходя из требований решаемой задачи.
Объемно-распределенная нагрузка определяется через градиент порового давления разрабатываемых слоев.
Поровое давление рассчитывают путем решения уравнения пьезопроводности полуаналитическим гидродинамическим методом.
Этапы г) и д) реализовывают в декартовых координатах или в цилиндрических координатах.
Дополнительно осуществляют визуализацию полученных значений напряжения и смещения.
Указанные признаки являются существенными и взаимосвязанными между собой причинно-следственной связью с образованием совокупности существенных признаков, достаточных для достижения указанного технического результата.
Для реализации способа определяют исходные параметры слоистой среды, к которым относятся, в частности, геологические параметры, характеризующие свойства и состояние геологической среды месторождения.
Следует отметить, что в рамках заявленного полуаналитического способа введены следующие важные предположения: слоистая среда представляет собой вертикальную последовательность горизонтальных плоскопараллельных слоев, где каждый слой является однородным и пороупругим. Такая слоистая модель названа здесь как Слоистая Геомеханическая Модель (СГМ).
«Плоскопараллельными» считаются слои, ограниченные с двух своих сторон параллельными горизонтальными плоскостями.
«Пороупругими» являются слои, в которых среда, их составляющая является пористой (например, содержит взаимосоединяющиеся поры, распределенные равномерно) и упругой (среда рассматривается как сплошная, а именно, непрерывно распределенная в пространстве и обладающая упругими свойствами). Такая среда характеризуется: упругими модулями λ, μ (коэффициенты Ламэ), пороупругим параметром Био а и толщиной каждого слоя (d). Такие параметры в рамках изобретения образуют геологические параметры.
«Однородными» считаются слои, в которых свойства, присущие слою, остаются неизменными (величина=соп81:) на всем протяжении слоя.
Схематичное изображение сечения слоистого трехмерного месторождения представлено на Фиг. 2.
Под объемно-распределенной нагрузкой понимается, в частности, распределение величины порового давления при добыче или нагнетании жидкости через скважины (добывающие* нагнетательные) в разрабатываемые слои месторождения.
Определение геологических параметров осуществляется путем их непосредственных измерений с помощью известных средств и методов, применяемых для подобных измерений в нефтегазовой отрасли. К таким средствам относятся, например, датчики, распределенные по соответствующим зонам разреза, сейсмические и акустические измерительные приборы, наклонометры, и иное скважинное измерительное оборудование. Так, например, толщины слоев d определяются как расстояния между литологическими границами, видимыми как скачки величин на каротажных кривых различной природы; упругие модули λ, μ можно измерить на основе анализа акустического каротажа; пороупругий параметр Био а можно оценить с помощью лабораторных механических испытаний извлеченного керна.
Вместе с тем, геологические параметры также можно установить на основании имеющихся знаний, исследований, полученных ранее, т.е. выбрать произвольным образом, исходя из требований решаемой задачи, а также определить эмпирическим путем или с использованием экспертных оценок.
Следующим шагом реализации способа является построение СГМ месторождения, а именно, полного решения для трехмерной слоистой среды. Такая модель описывается системой уравнений, характеризующих упругое напряженно-деформированное состояние равновесия для каждого слоя, т.е. совокупность внутренних напряжений и деформаций, возникающих при действии внешних нагрузок, а также уравнениями, определяющими граничные условия на границах слоев (между слоями, а также на верхней и нижней границах слоев).
Верхняя поверхность самого верхнего слоя свободна от нагрузок. Последний, самый нижний слой является полубесконечным в глубину с нулевыми смещениями на бесконечности. Граничные условия в боковых направлениях - также нулевые смещения на бесконечности. Граничные условия между слоями - непрерывность компонент напряжений, действующих на элементарную площадку интерфейсной границы (границы раздела), и непрерывность вертикальных смещений. Касательные смещения между слоями пропорциональны нормальной компоненте напряжений, действующему в том же направлении, умноженной на коэффициенты податливости, задаваемые при построении СГМ. По умолчанию, эти коэффициенты податливости равны нулю, чему соответствуют жестко сцепленные слои.
Каждый индивидуальный пороупругий слой может иметь независимо заданный профиль изменения объемно-распределенной нагрузки, который определяется давлением жидкости внутри него, р(х, у, z). Этот профиль давления представляется в виде произведения горизонтальной компоненты давления, которая изменяется как произвольная функция от х и у, и вертикальной компоненты, которая представляет собой полином 3-го порядка от z. Вышеперечисленные ограничения на геометрию и параметры модели необходимы для того, чтобы можно было построить быстрый полуаналитический метод. В объеме изобретения предположено, что слои имеют нумерацию сверху вниз ί = 1...п и ζ координата направлена вниз (Фиг. 2).
Далее, осуществляют получение аналитических решений указанных уравнений с последующим согласованием аналитических решений между собой с учетом граничных условий и получением искомых значений напряжения и смещения. Следует отметить, что в случае, когда уравнения, входящие в систему уравнений, представлены в инвариантной форме, а получение их аналитических решений реализуется, например, в декартовых или цилиндрических координатах, необходимо осуществить соответствую подстановку и заменить дифференциальные операторы в инвариантной форме их координатным представлением.
Вычисление смещений и напряжений внутри слоистого полупространства проводится с использованием 2-ух мерного преобразования Фурье путем решения соответствующей системы алгебраических уравнений. В конечном итоге, полное решение, при необходимости, преобразуется обратно в прямое пространство, чтобы получить распределения смещений и напряжений. Также возможно осуществить визуализацию полученных значений.
В качестве устройств, реализующих этапы по составлению массива уравнений, получения аналитических решений, их согласования и определения искомых значений, используются известные устройства для соответствующей математической обработки с техническими характеристиками, обеспечивающими получение заявленного технического результата. В качестве таких устройств применимы, например, компьютерные устройства, поддерживающие математические вычисления с вещественными числами с плавающей точкой, к которым, в частности, относятся Intel-совместимые компьютерные устройства архитектуры х86. Возможно использование компьютерных устройств, поддерживающих вычисления с плавающей точкой, других архитектур.
В качестве средств визуализации, позволяющих преобразовывать, полученные искомые величины в зрительные образы требуемого формата применимы различные программные продукты, реализующие графическую визуализацию математических данных, установленные и исполняемые на компьютерном оборудовании, снабженном средствами графического отображения информации, например, компьютерным дисплеем. В качестве нескольких примеров таких программных продуктов можно упомянуть Gnuplot, Techplot, Microsoft Excel, Matlab.
Точность аппроксимации в нашем случае определятся тем, насколько много дискретных точек используется для аппроксимации непрерывных функций, которые мы получаем как входные данные, и которые мы рассчитываем с помощью предложенного способа.
Следует отметить, что информация о распределении порового давления, необходимая для предсказания деформации, может быть получена из внешних источников, т.е. рассчитана по любому из известных для такой цели способов, или установлена эмпирическим путем. В то же время, возможно также использование различных конечно-разностных и конечно-элементных вычислительных методов для расчета давлений, или импорт распределений давления в ручном режиме. Так, в рамках данного изобретения, в качестве наиболее предпочтительного варианта осуществления способа, предложено рассчитывать изменения давления в соответствие с полуаналитическим способом оценки газовых залежей пласта (Gas Reservoir Evaluation and Assessment Tool), далее - GREAT (см. US7069148 B2, кл. G 06 F 19/00, опубл. 27.06.2006; Дж. Буссвелл, P. Банерджи, P.K.M. Тамбьянагам, Дж. Спаф. Обобщенное аналитическое решение для моделирования месторождений с несколькими скважинами и различными граничными условиями, СПЕ-99288, 2006. G. Busswell, R. Banerjee, R. .M. Thambynayagam, J. Spath, Generalized analytical solution for reservoir problems with multiple wells and boundary conditions, SPE- 99288, 2006). Этот метод основан на использовании оригинальных выражений для аналитических решений уравнения пьезопроводности для однофазной системы с однородными коэффициентами диффузии и позволяет получать решение такой задачи на порядки быстрее, чем конечно-разностные методы.
Связывание этих двух методов (СГМ и GREAT) с помощью специального программного модуля, рассматриваемого как часть данного изобретения, позволяет построить быстрый программный инструмент для связанного моделирования геомеханики резервуара. Функциональность, реализованная в разработанном программном модуле для связывания, может быть представлена в виде 3-х основных этапов следующих друг за другом:
1) экспорт информации о свойствах слоистого месторождения, необходимой при расчете полей давлений из СГМ в формат, поддерживаемый GREAT;
2) запуск GREAT с использованием экспортированных входных параметров;
3) импорт выходных данных GREAT, определяющих приложенную нагрузку со стороны давления и запуск СГМ для расчета смещений и напряжений.
Так, для определения распределения порового давления рассматриваемого месторождения с помощью GREAT необходимо предварительно определить гидродинамические параметры слоев, а также геолого-технологические параметры скважин месторождения.
Гидродинамические параметры слоев включают: начальное поровое давление слоя, проницаемость слоя (к), пористость слоя, начальную газо- или нефтенасыщенность, вязкость нефти, вязкость агента вытеснения и вытесняемой жидкости, начальную и конечную насыщенности агентом вытеснения.
Геолого-технологические параметры скважин включают: координаты пластопересечения скважины для каждого пласта, глубина кровли и подошвы каждого пласта в разрезе скважины, дебит скважины, динамика текущей и накопленной добычи нефти или газа, динамика текущей и накопленной закачки агента вытеснения.
Задание параметров задачи, моделируемой в GREAT, может быть выполнено с использованием либо входных файлов данных в текстовом формате с использованием ключевых слов, либо с использованием программного интерфейса приложения (API, Application Programming Interface) (GREAT, оформленного как динамически подключаемая библиотека (dll). Обмен данных через API существенно быстрее. Обмен данными через текстовые файлы более прозрачен для отладки и требует меньше модификаций при переходе на более новые версии GREAT. Процесс расчета в GREAT основан на суммировании вкладов от элементарных решений для точечных источников и плоских границ и не предполагает аппроксимацию на дискретной сетке, покрывающей всю расчетную область, чтобы уменьшить эффективную размерность задачи и достичь высокой скорости расчета. В связи с этим, GREAT рассчитывает значения давлений в точках, где расположены источники течения, тогда как для расчета деформаций в механический численный метод требуется ввести пространственное распределение давлений. Чтобы построить распределение давления, действующее как приложенная нагрузка в механическом численном методе, предложено ввести регулярную сетку заглушённых скважин с нулевыми потоками, и использовать их в качестве датчиков давления, как показано на Фиг. 4.
Входные данные для связанного СГМ+GREAT расчета организованы в виде текстового файла, состоящего из набора разделов. Пример входного файла приведен в разделе описания «Предпочтительное осуществление изобретения».
Раздел "FILES" определяет полный путь к исполняемым файлам GREAT и имя файла вывода.
Раздел "DIMENSIONS" описывает геометрию и размеры численной задачи в горизонтальном сечении. В частности, в разделе определяется размерность сетки скважин- датчиков давления, расстояние между ними и размер сетки в СГМ. Пористый коллектор, моделируемый GREAT, должен быть погружен внутрь упругой области большего размера, поскольку граничные условия для механической задачи должны быть заданы на достаточно большом расстоянии от источников приложенной нагрузки в модели. Расстояние от коллектора до внешней боковой границы задается в терминах множителей Вх, By определяющих отношение между расстоянием от границы коллектора до внешней границы и размером коллектора, как показано на Фиг. 4. Граничные условия на внешних границах в СГМ - нулевые смещения. Шаг сетки между узлами аппроксимации в СГМ выбран равным расстоянию между скважинами-датчиками давления в GREAT. Всего в СГМ расчете имеется (2Bx+l)-(2By+l)-M-N узлов, где M-N - размерность сетки скважин-датчиков.
Раздел "LAYERS" описывает последовательность слоев, начиная от свободной поверхности в направлении увеличения глубины, толщины слоев и механические свойства, как показано на Фиг. 2. Слои, в которых ведется разработка, меняется давление и где требуется считать поле давления, рассчитанное вне СГМ, помечаются флагом "PRESSURE_LOAD=l". Для таких слоев возможно два варианта загрузки внешних данных: 1) текстовый файл, содержащий таблицу из MxN значений давления; и 2) ссылка на входной файл задачи GREAT для данного слоя, которую нужно просчитать. Перед выполнением GREAT во втором случае, входной файл расширяется путем добавления в его конец автоматически сгенерированной сетки из MxN вертикальных закрытых скважин- датчиков давления, перфорированные секции таких скважин помещаются в середину разрабатываемых слоев.
Раздел "OUTPUT" определяет список z координат слоев, в которых требуется рассчитать выходные данные. Основные этапы алгоритма представлены на Фиг. 5.
Анализ результатов, полученных после выполнения одного СГМ+GREAT моделирования (Фиг. 1) включает в себя верификацию результатов, которая может представлять собой сравнение с другими эталонными решениями, либо во всей, либо в какой-то части области определения модели, или же может ограничиваться только проверкой того, что решение попадает внутрь некого допустимого интервала решений. Если решение не является удовлетворительным, то весь процесс возвращается к первоначальному шагу на Фиг. 1 , чтобы провести коррекцию входных параметров, определяющих модель. Затем процесс моделирования, представленный на Фиг. 1 , повторяется. Подобное повторение нескольких запусков моделирования с коррекцией входных данных проводится и при постановке задач на адаптацию по параметрам. В этом случае получаемое решение сравнивается с набором целевых параметров, выбираемых оператором, и варьирование входных параметров продолжается до тех пор, пока не будет получено соответствие целевым значениям в пределах некоторой точности. Проведение большого количества СГМ+GREAT моделирований также необходимо при проведении параметрических исследований для определения того, как варьирование входных параметров влияет на конечное решение.
Уравнение для отдельного пороупругого слоя, характеризующее упругое напряженно-деформированное равновесие слоя, представлено в виде:
Figure imgf000019_0001
где V - дифференциальный оператор, и - вектор смещения геологической среды в выбранных точках слоя, λ, μ - упругие коэффициенты Ламэ, f - объемно-распределенная нагрузка.
Уравнение (1) записано в инвариантной форме. В дальнейшем, (1) будет решаться в декартовых координатах, в которых ось z ориентирована перпендикулярно плоскости слоев. Решение строится с учетом следующих предположений: вся внешне приложенная нагрузка возникает только за счет градиента порового давления , где а
Figure imgf000020_0003
- пороупругая константа; давление может быть разбито на произведение z функции и (х, у) функции; изменение давления в z направлении описывается полиномом 3-ей степени. Уравнение (1) может быть преобразовано к виду:
Figure imgf000020_0001
где введены новые обозначения для безразмерного коэффициента β, для вариации давления в ζ направления ρ(ζ), и (х, у) вариации давления h(x, у). Для решения системы из 3 уравнений второго порядка (2) в 3D, используется метод потенциальной функции с последующим интегральным преобразованием как, например, см. В.М. Ентов, Т.А. Малахова, Об изменении напряженно-деформированного состояния горных пород при изменении давления в насыщенном жидкостью пласте, Известия РАН, Механика твердого тела, N° 6, стр. 53-65, 1974; а также И.Н.Снеддон, Преобразования Фурье, Делавэр, Нью-Йорк 1995 (I.N. Sneddon, Fourier Transforms, Dover, New York, 1995), где подобный подход был использован для построения осесимметричного решения.
Анзатц решения с использованием потенциальной функции Ф и бездивергентного векторного поля Ψ:
Figure imgf000020_0002
Figure imgf000021_0001
где
Figure imgf000021_0002
После подстановки (3-5) в (2), можно получить эквивалентный набор уравнений:
Figure imgf000021_0005
Введение бездивергентного векторного поля Ψ в выбранной форме (3-4) позволяет разделить исходное уравнение (2) в систему независимых уравнений (8-9).
Решение (8-9) ищется с использованием двумерного интегрального преобразования Фурье в плоскости (х, у):
Figure imgf000021_0003
с базисными функциями
Figure imgf000021_0006
После по становки (10-11) в (8-9) можно получить
Figure imgf000021_0004
где волновой вектор к образ функции источника h {m,n) определены в виде
Figure imgf000022_0002
Решение обыкновенных дифференциальных уравнений под интегралами (13-14) позволяет получить общий вид кернфункций:
Figure imgf000022_0001
Чтобы получить действительные значения смещений в прямом пространстве, необходимо ввести 6 граничных условий, которые бы зафиксировали свободные параметры АГА6, выполнить дифференцирование в (3-5), и применить прямое преобразование Фурье в (10-11). Напряжения можно найти путем дифференцирования смещений по пространству и использовать закон Гука.
Чтобы построить полное решение для трехмерной слоистой среды, необходимо согласовать между собой индивидуальные аналитические решения в каждом слое путем надлежащего выбора и решения контактных граничных условий на интерфейсных плоскостях между слоями. В 3D требуется иметь по 6 условий между каждой парой соседних слоев и по 3 условия на верхней и нижней внешних границах. В нижнем направлении все компоненты смещений равны нулю при
Figure imgf000022_0003
На верхней свободной поверхности приложены нулевые усилия:
Figure imgf000023_0003
Предполагается, что локальная система координат внутри каждого пороупругого слоя центрирована, т.е. плоскость ζ - 0 расположена в середине слоя; точки (х - 0; у = 0) в разных слоях лежат вдоль одной вертикальной оси. Верхний индекс у переменных σ и и обозначает номер слоя; нумерация слоев i— \ ...п идет сверху вниз (см. Фиг. 2). Внутренние границы представляют собой идеальные упругие контакты между слоями и идеально-непроницаемые барьеры для перетоков жидкости. Следовательно, на интерфейсах (поверхностях) отсутствует скачок следующих компонент напряжений:
Figure imgf000023_0001
Для компонент смещений возможен только скачок в касательном направлении, который полностью определяется значением сдвигового напряжения на интерфейсной границе:
Figure imgf000023_0002
где С, - коэффициент податливости интерфейсного слоя в касательном направлении, изотропный в плоскости (х, у).
Граничные условия (20, 21) должны выполняться в любой точке (х, у). Можно показать, что это эквивалентно требованию того, что (20, 21) должны также выполняться для Фурье образов физических компонент в каждой точке (т, п). После подстановки (10, 11, 18) в (3- 5), и использования для выражения компонент
Figure imgf000024_0003
напряжений:
Figure imgf000024_0001
μ μ
Подставляя выражение для кернфункции G (17) и группируя члены при коэффициентах А ГА6:
Figure imgf000024_0002
Figure imgf000025_0001
где введены новые обозначения для функций при источниковых членах й(т,п) , учитывающих полиномиальное изменение распределения давления в ζ-направлении:
Figure imgf000025_0002
В конечном выражении для граничных условий (20, 21), используемых для расчета коэффициентов АгАб, использованы эквивалентные линейные преобразования вида:
Figure imgf000025_0003
Используя эти выражения, контактные граничные условия (20- 21) преобразуются в :
Figure imgf000026_0001
Следует отметить, что для краткости изложения, dj в (27-32) обозначает половину толщины слоя / в отличие от толщин слоев, изображенных на Фиг. 2. Уравнения (27-32) выражают граничные условия на внутренних интерфейсах между слоями при i—\ ...n-\ . Чтобы задать свободные граничные условия (19) на поверхности следует использовать выражения (28, 29, 31) при i = и затем подставить А°...А^ = 0 , чтобы оставить только члены из первого верхнего слоя. Нижний слой ί = η имеет бесконечную толщину = оо. В этом слое удобно выбрать внутреннюю локальную координату ζ - 0 на верхней границе этого слоя, а не в его центре, как в слоях конечной толщины. В этом случае, система (27-32) по-прежнему справедлива для последней межслойной границы при / = п, если подставить dn— 0 и А" = А1 = А" = 0, что вытекает из нулевых граничных условий на бесконечности.
После соответствующего учета особенностей граничных условий на первой и последней межслойной границе можно построить ленточную 4-диагональную матрицу размерности 6*(и-1)+3, которая является полной и позволяет вычислить коэффициенты Aj-A6 внутри всех слоев для каждого Фурье-индекса (т, п). Подсистема уравнений содержащих As~A6 полностью независима от системы А \-А4. Поскольку в первой подсистеме A s-A6 источниковый член в правой части системы уравнений, определяющий внешнее нагружение, является нулевым, решение этой подсистемы тривиально: А5 = А6 = 0. Таким образом, матрица граничных условий сокращается до 4*(и-1)+2 уравнений (27- 30), которые можно разрешить любым алгоритмом по обращению матриц.
Схематично (Фиг. 3) расчет компонент смещений и напряжений в заданных точках (х, у, z)j требует выполнения следующих шагов: 1) выполнение прямого интегрального преобразования распределений давлений в слоях с приложенной внешней нагрузкой; 2) обращение матрицы граничных условий для построения решения в преобразованном пространстве; 3) выполнение обратного интегрального преобразования построенного решения в заданных точках вывода (х, у, ζ
При фактическом проведение численных расчетов более удобно заменить континуальное интегрирование преобразований (10, 11, 16) на быстрое преобразование Фурье (Fast Fourier Transform, FFT), которое эквивалентно численному интегрированию по методу Эйлера в прямоугольной области. Это подразумевает существование регулярной сетки из MxN аппроксимационных узлов и прямоугольной расчетной области размерности Ах-М на Ay-N, в которой пространственные шаги сетки выбираются пользователем исходя из требуемой точности аппроксимации.
Преобразование Фурье является наиболее затратной частью алгоритма, поскольку оно требует (X V-log МАО операций для каждого слоя с ненулевым распределением давления, для каждой рассчитываемой компоненты ( σ#) и каждой уникальной ζ координаты, в которой рассчитываются выходные данные. Обращение матрицы граничных условий требует только 0(L MN) операций, где L - число слоев.
Очень важно использовать максимально эффективный FFT алгоритм. Одним из вариантов является использование хорошо зарекомендовавших себя и используемых в настоящее время мульти- платформенных исполняемых библиотек, как например, FFTW Library (FFTW-библиотека), www.fftw.org; США CUFFT Library, NVIDIA Corp., 2008, http://www.nvidia.com/object/cuda_develop.html. Возможно также использование любых других имеющихся или разрабатываемых FFT библиотек, которые наилучшим образом соответствуют потребностям пользователя.
Краткое описание чертежей
Фиг. 1 - Этапы способа.
Фиг. 2 - Схематичное изображение сечения слоистого трехмерного месторождения.
Фиг. 3 - Основные этапы расчетного алгоритма для СГМ.
Фиг. 4 - Схематичное изображение прямоугольного коллектора с сеткой скважин-датчиков давления, погруженного в окружающие упругие породы.
Фиг. 5 - Основные этапы алгоритма связанного СГМ+GREAT расчета.
Фиг. 6 - Иллюстрация моделируемой задачи с пятиточечной схемой разработки в двухслойном коллекторе.
Фиг. 7 - Распределение значений вертикальных смещений Uz на верхней границе верхнего разрабатываемого слоя на глубине z=6800 футов.
Фиг. 8 - Распределение значений нормальных горизонтальных напряжений σ~ на глубине ζ=6850 футов.
Предпочтительное осуществления изобретения
В качестве примера численной реализации предложенного способа быстрого определения напряженно-деформированного состояния слоистой среды, рассмотрим определение смещений и напряжений, вызванных разработкой месторождения с заводнением по классической пятитоточечной схеме. В рамках примера рассмотрим единичный пятиточечный шаблон, состоящий из четырех добывающих скважин и одной нагнетающей.
Определим исходные параметры, рассматриваемой в примере 3-х мерной задачи по геомеханическому моделированию разработки по пятиточечной схеме. Значения параметров выбраны произвольно в рамках возможных значений свойств горных пород и типичных характеристик нефтегазовых месторождений. Геометрия задачи проиллюстрирована на Фиг. 6. Графическая визуализация на рисунке выполнена с использованием программных продуктов, принадлежащих компании Schlumberger, которые не имеют непосредственного отношения к данной заявке. Различные типы штриховки соответствуют различным значениям проницаемости горных пород (к), которая изменяется от 0 до 1 мД. Прямоугольный коллектор с пористостью 20% состоит из 2 разрабатываемых слоев проницаемостью 0.5 мД (верхний) и 1 мД (нижний) соответственно. Между разрабатываемыми проницаемыми слоями находится непроницаемый барьер. Прямоугольный коллектор имеет размер 10000 на 10000 футов, он окружен со всех сторон непроницаемыми упругими слоями, которые необходимы для корректного задания граничных условий на большом расстоянии от нагружаемой центральной части задачи. В данном случае на верхней свободной поверхности заданы нулевые напряжения, на боковых и на нижней границе заданы нулевые нормальные смещения и скольжение относительно тангенциальных смещений. Оба разрабатываемых слоя пересечены четырьмя добывающими скважинами, обозначенными Р-1 - Р-4, и одной нагнетающей, обозначенной 'INJ'. Геологические параметры слоев, а также их проницаемости, и параметр Био приведены в Таблице 1. Отметим, что упругие свойства слоев заданы в терминах модуля Юнга Е и коэффициента Пуассона v, которые связаны взаимооднозначным алгебраическим соотношением с параметрами Ламэ λ, μ. Полный набор входных параметров продублирован в примере входного файла для СГМ-расчета и в примере входного файла для GREAT, которые представлены ниже.
Figure imgf000031_0001
Определение внешней нагрузки, приложенной к слоистой геологической среде месторождения в результате его разработки по пятиточечной схеме в рассматриваемом примере, определяется с помощью расчета давления программным комплексом GREAT. На скважинах задан постоянный поток жидкости, равный 800 баррелей/день на добывающих и 1600 баррелей/день на нагнетающей скважине. Моделирование разработки ведется на протяжении 1440 дней, изменение смещений и напряжений рассчитывается по завершении этого периода. Расчет распределения давления в двухслойном коллекторе проводится с помощью двух независимых запусков GREAT по отдельности для каждого слоя. При этом, потоки через интервалы скважин, пересекающих пласты, распределяются пропорционально проницаемостям и мощностям слоев в пропорции 1/4 для верхнего слоя и 3/4 для нижнего. Распределение давления рассчитывается на равномерной сетке из 50 на 50 фиктивных скважин с нулевым потоком в центральной части задачи. Затем рассчитанное поле давления конвертируется в формат, который может быть загружен в СГМ и используется как механическая нагрузка, приложенная к слоям. Полный набор входных параметров продублирован в примере входного файла для GREAT (см. ниже).
Далее, заданные исходные параметры и распределение нагрузки подставляются в составленную систему уравнений из уравнений упругого равновесия и граничных условий. Аналитические решения для каждого слоя получены в Фурье-пространстве и записаны в виде интегрального преобразования, поэтому выполняется быстрое преобразование Фурье, чтобы получить решения для отдельных слоев при заданных значениях их параметров. После этого проводится согласование решений в каждом слое для того, чтобы добиться выполнения граничных условий между слоями. Для этого заполняется матрица уравнений граничных условий, которая затем обращается с использованием точного неитерационного алгоритма обращения матриц. В результате этого обращения определяются свободные коэффициенты аналитических решений, записанных для каждого слоя. В дальнейшем, полученное полное решение, справедливое одноврмененно для всех слоев преобразуется с помощью обратного быстрого преоборазования Фурье, чтобы получить значения искомых смещений и напряжений в декартовом пространстве. На Фиг. 7 и Фиг. 8 приведены примеры рассчитанных вертикальных смещений и нормальны компонент тензора напряжений в направлении одной из горизонтальных координатных осей, которые были получены при вышеупомянутых значениях входных параметров. На Фиг. 7 и Фиг. 8 видно, что максимальных значений указанные компоненты смещений и напряжений достигают вблизи работающих скважин. Так, вертикальные смещения и2 составляют около -9.5 см вблизи добывающих скважин и около 8 см у нагнетающей скважины. Компонента тензора напряжений составляет порядка -4.2 МПа у добывающих скважин и около 7.8 МПа у нагнетающей скважины.
Пример входного файла СГМ:
Figure imgf000033_0001
Figure imgf000034_0001
Figure imgf000035_0001
Пример входного файла "5spot_GR_top.txt" определяющего расчет давления с помощью GREAT в верхнем пористом разрабатываемом слое. Входной файл для нижнего разрабатываемого слоя "5spot_GR_btm.txt" аналогичен.
Figure imgf000036_0001
Figure imgf000037_0001
Промышленная применимость
Промышленная применимость изобретения подтверждается приведенным выше примером осуществления способа, а также возможностью его реализации при использовании широко известных в промышленности и других областях технологического оборудования и материалов.

Claims

Формула
1. Способ определения напряженно-деформированного состояния слоистой среды, включающий этапы:
а) определение исходных параметров слоистой среды,
б) определение объемно-распределенной нагрузки разрабатываемых слоев слоистой среды,
в) составление системы уравнений, включающей:
-уравнения упругого равновесия для каждого слоя,
-уравнения, определяющие граничные условия на границах слоев: непрерывность напряжений между слоями, нулевые напряжения на свободной границе, нулевые смещения на бесконечности,
г) получение аналитических решений уравнений упругого равновесия в каждом слое, при условии, что слои являются плоскопараллельными, однородными и пороупругими,
д) согласование между собой аналитических решений, полученных в разных слоях, с учетом соответствующих граничных условий,
е) определение значений напряжения и смещения в выбранных точках среды в соответствии с требуемой точностью аппроксимации на основе согласованных аналитических решений.
2. Способ по п.1, отличающийся тем, что исходные параметры включают геологические параметры слоев.
3. Способ по п.2, отличающийся тем, что геологические параметры слоев включают: упругие коэффициенты Ламэ (Я, μ), пороупругий параметр Био, толщину слоя.
4. Способ по п.З, отличающийся тем, что параметры измеряют с помощью средств измерения и/или устанавливают эмпирическим путем и/или выбирают произвольным образом, исходя из требований решаемой задачи.
5. Способ по п.1, отличающийся тем, что объемно-распределенная нагрузка определяется через градиент порового давления разрабатываемых слоев.
6. Способ по п.5, отличающийся тем, что поровое давление рассчитывают путем решения уравнения пьезопроводности полуаналитическим гидродинамическим методом.
7. Способ по п.1, отличающийся тем, что этапы г) и д) реализованы в декартовых координатах.
8. Способ по п.1, отличающийся тем, что этапы г) и д) реализованы в цилиндрических координатах.
9. Способ по п.1 , отличающийся тем, что дополнительно осуществляют визуализацию полученных значений напряжения и смещения.
PCT/RU2010/000299 2010-06-08 2010-06-08 Способ определения напряженно-деформированного состояния слоистой среды WO2011155862A1 (ru)

Priority Applications (2)

Application Number Priority Date Filing Date Title
PCT/RU2010/000299 WO2011155862A1 (ru) 2010-06-08 2010-06-08 Способ определения напряженно-деформированного состояния слоистой среды
US13/701,447 US20130151160A1 (en) 2010-06-08 2010-06-08 Method for determining the stress-strain state of a stratified medium

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/RU2010/000299 WO2011155862A1 (ru) 2010-06-08 2010-06-08 Способ определения напряженно-деформированного состояния слоистой среды

Publications (1)

Publication Number Publication Date
WO2011155862A1 true WO2011155862A1 (ru) 2011-12-15

Family

ID=45098280

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/RU2010/000299 WO2011155862A1 (ru) 2010-06-08 2010-06-08 Способ определения напряженно-деформированного состояния слоистой среды

Country Status (2)

Country Link
US (1) US20130151160A1 (ru)
WO (1) WO2011155862A1 (ru)

Cited By (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113505484A (zh) * 2021-07-13 2021-10-15 西南科技大学 深部采场底板岩层的最大破坏深度确定方法
CN116305451A (zh) * 2023-03-02 2023-06-23 中国地质大学（北京） 连续-非连续地质模型建立方法及装置

Families Citing this family (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US8374836B2 (en) * 2008-11-12 2013-02-12 Geoscape Analytics, Inc. Methods and systems for constructing and using a subterranean geomechanics model spanning local to zonal scale in complex geological environments
CN110017135B (zh) * 2019-02-15 2022-05-20 西南石油大学 一种裂缝性地层井壁裂缝扩展压力预测方法
CN111324960A (zh) * 2020-02-26 2020-06-23 海南大学 一种含有浅埋隐伏球状空洞地基承载力的量化方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SU953506A1 (ru) * 1981-01-16 1982-08-23 Ярославский политехнический институт Способ определени напр жени и деформации при расслаивании слоистых резинокордных материалов
SU1442904A1 (ru) * 1987-05-12 1988-12-07 Предприятие П/Я А-7840 Способ определени физико-механических характеристик слоистых анизотропных материалов
RU1810783C (ru) * 1991-03-04 1993-04-23 Луганский Машиностроительный Институт Способ испытани плоских образцов слоистых материалов
CN1511247A (zh) * 2001-05-25 2004-07-07 加州理工学院 确定包括体积力作用的层状和渐变结构的大形变和应力

Family Cites Families (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US7177764B2 (en) * 2000-07-14 2007-02-13 Schlumberger Technology Corp. Simulation method and apparatus for determining subsidence in a reservoir
GB0017227D0 (en) * 2000-07-14 2000-08-30 Schlumberger Ind Ltd Fully coupled geomechanics in a commerical reservoir simulator
GB2372567B (en) * 2001-02-22 2003-04-09 Schlumberger Holdings Estimating subsurface subsidence and compaction

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
SU953506A1 (ru) * 1981-01-16 1982-08-23 Ярославский политехнический институт Способ определени напр жени и деформации при расслаивании слоистых резинокордных материалов
SU1442904A1 (ru) * 1987-05-12 1988-12-07 Предприятие П/Я А-7840 Способ определени физико-механических характеристик слоистых анизотропных материалов
RU1810783C (ru) * 1991-03-04 1993-04-23 Луганский Машиностроительный Институт Способ испытани плоских образцов слоистых материалов
CN1511247A (zh) * 2001-05-25 2004-07-07 加州理工学院 确定包括体积力作用的层状和渐变结构的大形变和应力

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113505484A (zh) * 2021-07-13 2021-10-15 西南科技大学 深部采场底板岩层的最大破坏深度确定方法
CN116305451A (zh) * 2023-03-02 2023-06-23 中国地质大学（北京） 连续-非连续地质模型建立方法及装置
CN116305451B (zh) * 2023-03-02 2024-01-09 中国地质大学（北京） 连续-非连续地质模型建立方法及装置

Also Published As

Publication number Publication date
US20130151160A1 (en) 2013-06-13

Similar Documents

Publication Publication Date Title
Cacace et al. Modelling of fractured carbonate reservoirs: outline of a novel technique via a case study from the Molasse Basin, southern Bavaria, Germany
US10846447B2 (en) Method and system for stacking fracture prediction
US6662109B2 (en) Method of constraining by dynamic production data a fine model representative of the distribution in the reservoir of a physical quantity characteristic of the subsoil structure
US8768672B2 (en) Method for predicting time-lapse seismic timeshifts by computer simulation
Zambrano et al. Analysis of fracture roughness control on permeability using SfM and fluid flow simulations: implications for carbonate reservoir characterization
BRPI0714028A2 (pt) métodos para refinar uma propriedade fìsica e para produzir hidrocarbonetos a partir de uma região de subsolo
Shahin et al. Viscoplastic interpretation of localized compaction creep in porous rock
WO2011155862A1 (ru) Способ определения напряженно-деформированного состояния слоистой среды
Tsang et al. Modeling shear rigidity of stratified bedrock in site response analysis
Al-Zubaidi et al. 3D mechanical earth model for Zubair oilfield in southern Iraq
Abdallah et al. Compaction banding in high‐porosity carbonate rocks: 2. A gradient‐dependent plasticity model
US8718992B2 (en) Method for history matching of a geological model comprising a sub-seismic fault network
Das et al. Compressibility predictions using digital thin-section images of rocks
Krylov et al. Numerical modeling of nonlinear response of seafloor porous saturated soil deposits to SH-wave propagation
Aji et al. 3D hybrid model of foundation‐soil‐foundation dynamic interaction
Cier et al. Automatically adaptive stabilized finite elements and continuation analysis for compaction banding in geomaterials
Benetatos et al. Fully integrated hydrocarbon reservoir studies: myth or reality?
CN107609265A (zh) 一种基于蚂蚁追踪的地层应力场有限元模拟方法及系统
Congro et al. Determination of fault damage zones in sandstone rocks using numerical models and statistical analyses
Yang Calibrated Fragility Functions for Seismic Loading of Sacramento-San Joaquin River Delta Levees
Duda et al. Anisotropic poroelastic modelling of depletion-induced pore pressure changes in Valhall overburden
Panza et al. Earthquakes site effects modeling by hybrid MS-BIEM: the case study of Sofia, Bulgaria
Vasco et al. On the use of adjoints in the inversion of observed quasi-static deformation
Sekhavatian et al. Reliability analysis of deep excavations by RS and MCS methods: case study
Latychev et al. On the compliance method and the assessment of three-dimensional seafloor gas hydrate deposits

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 10852973

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

WWE Wipo information: entry into national phase

Ref document number: 13701447

Country of ref document: US

122 Ep: pct application non-entry in european phase

Ref document number: 10852973

Country of ref document: EP

Kind code of ref document: A1