A method for simulating fracture propagation based on an adaptive multi-scale phase field model

Abstract

The application discloses a fracturing fracture extension simulation method based on an adaptive multi-scale phase field model and relates to the field of oil and gas reservoir reconstruction. The technical scheme is as follows: parameters of reservoir physical properties and rock mechanics are acquired, and a reservoir geological model is established according to the parameters; the reservoir geometric model is divided into a plurality of finite element grids, and the grid attributes of the finite element grids are valued by using the parameters, so that preset hydraulic fractures in the finite element grids and initial variable values on nodes of the finite element grids are obtained; a multi-scale adaptive grid for dynamic expansion of the hydraulic fractures is constructed according to the position relationship between the preset hydraulic fractures and the finite element grids; a coupling mathematical model of stress field, seepage field and phase field of hydraulic fracturing is established, a numerical calculation expression of the coupling mathematical model is constructed according to the multi-scale adaptive grid; and the coupling mathematical model is solved based on the numerical calculation expression, with the parameters of reservoir physical properties and rock mechanics and the boundary conditions of the multi-scale adaptive grid as constraints, so that the hydraulic fracturing fracture extension simulation results under different reservoirs and construction conditions are obtained.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of oil and gas reservoir reconstruction, more particularly, it relates to a fracturing fracture extension simulation method based on an adaptive multi-scale phase field model. BACKGROUND

[0002] Hydraulic fracturing technology is a key means to realize efficient development of oil and gas resources. By pumping high-pressure fluid into the reservoir to build artificial fractures, the percolation state of the reservoir is changed, and the oil and gas production of a single well is improved. The hydraulic fracture propagation in the reservoir is affected by various complex factors such as in-situ stress characteristics, rock mechanical properties, and fluid flow properties. Revealing the hydraulic fracture propagation law under the complex action of geology and engineering has important theoretical significance and engineering value for optimizing hydraulic fracturing design, improving reservoir reconstruction effect, and reducing environmental risk. With the development of unconventional oil and gas resources such as shale oil and gas and tight oil and gas, the demand for hydraulic fracture propagation simulation under complex geological conditions is also increasing. The display crack model represented by the finite element method and the boundary element method has difficulty in effectively dealing with complex expansion behaviors such as crack deflection, branching, and intersection. The phase field method describes the geometric shape of the crack by introducing a continuous phase field variable, avoids the complexity of explicitly tracking the crack interface in traditional methods, and can naturally handle complex expansion situations such as crack deflection, branching, and intersection, providing a more flexible framework for hydraulic fracture propagation simulation under complex geological conditions.

[0003] The phase field model describes the evolution process of the system through the minimization of the free energy functional. When using the phase field model to simulate hydraulic fracture propagation, the phase field variable has a high gradient near the hydraulic fracture. The quality and resolution of the grid will directly affect the accuracy and efficiency of the numerical simulation. To balance the numerical calculation accuracy and efficiency, adaptive grid technology is widely used in the solution of the phase field model of hydraulic fracture propagation. Local grid refinement schemes are used in the crack area, while grid coarsening schemes are used in other areas. When simulating the dynamic propagation of hydraulic fractures, adaptive grid technology needs to dynamically adjust the density and distribution of the grid to capture the changes in the phase field value, increasing the complexity of algorithms such as grid generation, interpolation, and data transfer. In addition, adaptive grid technology often uses unstructured grid schemes, making grid generation and adjustment more difficult. Frequent changes in the grid can cause problems such as grid quality degradation, reduced computational efficiency, and numerical calculation oscillation. Although adaptive grid technology has significant advantages in improving the computational efficiency of the phase field model, it still has many shortcomings in terms of technical implementation, computational efficiency, and accuracy control. Therefore, it is necessary to develop simpler, more efficient, and more stable adaptive grid algorithms. SUMMARY

[0004] The application aims to provide a fracturing fracture extension simulation method based on an adaptive multi-scale phase field model, solve the problems of great difficulty in realizing an adaptive grid algorithm, difficulty in dynamic grid division, and insufficient numerical solution stability in the prior art, realize efficient and stable solution of a multi-physical field coupling problem by using a multi-scale mixed finite element grid under a finite element method, accurately capture the complex expansion behavior of a hydraulic fracture by a phase field method to optimize design of a hydraulic fracturing construction parameter, and achieve the goal of improving fracturing reconstruction effect of an oil and gas reservoir and reducing environmental risk.

[0005] The application provides a fracturing fracture extension simulation method based on an adaptive multi-scale phase field model, and the method comprises the following steps:

[0006] Reservoir physical properties and rock mechanics parameters are acquired, and a reservoir geological model is established according to the parameters;

[0007] A reservoir geometry model is created based on the reservoir geological model and a preset hydraulic fracturing position attribute, the reservoir geometry model is divided into a plurality of finite element grids, and the grid attribute of the finite element grid is valued by using the parameters, so that an initial variable value of the preset hydraulic fracture in the finite element grid and the finite element grid node is obtained;

[0008] A multi-scale adaptive grid of dynamic expansion of the hydraulic fracture is constructed according to the position relationship between the preset hydraulic fracture and the finite element grid;

[0009] A coupling mathematical model of a stress field, a seepage field and a phase field of hydraulic fracturing is established, and a numerical calculation expression of the coupling mathematical model is constructed according to the multi-scale adaptive grid;

[0010] Based on the numerical calculation expression, the reservoir physical properties and the rock mechanics parameters and the boundary grid node variable value of the multi-scale adaptive grid are taken as constraints to solve the coupling mathematical model, so that a hydraulic fracturing fracture extension simulation result under different reservoir and construction conditions is obtained.

[0011] In an implementation scheme, the parameters include a horizontal maximum principal stress, a horizontal minimum principal stress, a reservoir pore pressure, a reservoir thickness, a matrix porosity, a matrix permeability, a Young's modulus, a Poisson's ratio, a matrix bulk modulus, a skeleton bulk modulus and a rock tensile strength.

[0012] In an implementation scheme, the reservoir geological model is established according to the parameters, and specifically, the reservoir geological model is established according to a mapping relationship between the spatial position and the physical attribute of the parameters.

[0013] In an implementation scheme, the reservoir geometry model is created based on the reservoir geological model and the preset hydraulic fracturing position attribute, and the reservoir geometry model is divided into a plurality of finite element grids, and specifically:

[0014] selecting a rectangular or cuboid region from the reservoir geological model according to the hydraulic fracture location attribute to create a reservoir geometric model;

[0015] The reservoir geometric model is meshed by using finite elements to obtain a calculation domain composed of finite element meshes.

[0016] In an implementation scheme, the grid attributes of the finite element meshes are valued by using parameters to obtain preset hydraulic fractures in the finite element meshes and initial variable values on the nodes of the finite element meshes, including:

[0017] The hydraulic fractures are preset in the finite element meshes; wherein the hydraulic fractures include one-dimensional line segments or two-dimensional planes, and the hydraulic fractures do not coincide with the interfaces of the finite element meshes;

[0018] Initial variable values are preset on the nodes of the finite element meshes, and the initial variable values include initial displacement values, initial pressure values and initial phase field values.

[0019] In an implementation scheme, the initial displacement values are calculated by applying stress boundary conditions, the initial pressure values are equal to reservoir pore pressure values, and the initial phase field values are calculated by using hydraulic fracture trajectories and historical state variables.

[0020] In an implementation scheme, a multi-scale adaptive mesh for dynamic expansion of the hydraulic fractures is constructed according to the position relationship between the preset hydraulic fractures and the finite element meshes, including:

[0021] The crack tip elements and the non-crack tip elements are determined by using the relative position relationship between the preset hydraulic fractures and the finite element meshes;

[0022] The center of the finite element mesh in which the crack main body region is located is taken as a circle point, a pre-configured critical radius value is taken as a radius to perform region division and set intersection, the finite element meshes in the region are identified as refined elements, and the finite element meshes outside the region are identified as non-refined elements;

[0023] The non-refined elements adjacent to the refined elements are searched as transition elements, and the remaining non-refined elements that are not adjacent are taken as coarsened elements; wherein the number of nodes of the refined elements and the coarsened elements is the same, there is a multiple relationship between the sizes of the refined elements and the coarsened elements, and the size of the refined elements matches the length scale parameter of the phase field.

[0024] In an implementation scheme, a coupling mathematical model of the stress field, the seepage field and the phase field of the hydraulic fracturing is established, and a numerical calculation expression of the coupling mathematical model is constructed according to the multi-scale adaptive mesh, and specifically:

[0025] The numerical base functions, the local stiffness matrix and the local equivalent load matrix of the transition elements are constructed;

[0026] The local stiffness matrix and the local equivalent load matrix of the stress field, the seepage field and the phase field of the calculation refinement unit and the calculation coarsening unit are calculated, the local stiffness matrix and the local equivalent load matrix of the transition unit are combined, the global stiffness matrix and the equivalent load matrix of the multi-scale adaptive grid are assembled, and the control equations of the stress field, the seepage field and the phase field are composed according to the global stiffness matrix and the equivalent load matrix.

[0027] The node value conversion between the local fine grid and the coarse grid of the transition unit is realized through the numerical basis function, the Newton-Raphson iteration method is used to iteratively solve the control equations of the stress field, the seepage field and the phase field, and the trial solution of each iteration step is obtained, and the numerical solution of the coupled mathematical model is realized.

[0028] If the trial solution of the adjacent iteration step at the current time step meets the convergence condition, the solution result of the numerical solution converges, and the displacement value, the pressure value and the phase field value of the node are output.

[0029] In an implementation scheme, the numerical basis function, the local stiffness matrix and the local equivalent load matrix of the transition unit are constructed, including:

[0030] Taking a single transition unit as a construction object, the fine grid and the coarse grid are locally divided, and the solid phase stiffness matrix of the stress field, the seepage matrix of the seepage field and the evolution matrix of the phase field are taken as the left end items to construct the local balance equation.

[0031] The boundary number and the node number of the transition unit are determined according to the connection condition of the transition unit and the adjacent coarsening unit and refinement unit, and the boundary of the transition unit sharing the same node is taken as a group to apply the linear boundary condition.

[0032] The numerical basis function of the transition unit in displacement, pressure and phase field is sequentially solved by combining the local balance equation of the transition unit and the linear boundary condition.

[0033] The finite element stiffness matrix, the equivalent load matrix and the node value conversion relationship between the coarse grid and the fine grid are established through the numerical basis function.

[0034] In an implementation scheme, based on the numerical calculation expression, the parameters of reservoir properties and rock mechanics and the boundary conditions of the multi-scale adaptive grid are taken as constraints to solve the coupled mathematical model, and the simulation results of the hydraulic fracturing crack extension under different reservoirs and construction conditions are obtained, specifically:

[0035] The initial variable value preset on the finite element grid node is taken as input, the numerical calculation expression is iteratively solved in combination with the reservoir properties and rock mechanics parameters and the boundary grid node variable value, and the trial solution of the displacement, the pressure and the phase field value of the grid node at the current time step is obtained.

[0036] After obtaining the progressive solution of the phase field value of the grid node at the current time step, the crack tip element and the crack non-tip element are divided according to the phase field value of the grid element, different critical radius values are used for element refinement operation and element coarsening operation, the global grid and node number are adaptively adjusted, and the adaptive division of the multi-scale grid in the dynamic expansion process of the hydraulic fracture is realized.

[0037] If the trial solution of the adjacent iteration step at the current time step satisfies the convergence condition, the solution converges, the displacement value, the pressure value and the phase field value of the node are output, and the calculation of the next time step is started.

[0038] The calculation is repeated until the simulation termination condition is met, the displacement, pressure and phase field calculation results of the multi-scale adaptive grid node are output and data visualization processing is performed, and the dynamic extension simulation results of the hydraulic fracture under different reservoirs and construction conditions are obtained.

[0039] Compared with the prior art, the present application has the following beneficial effects:

[0040] (1) The introduction of the multi-scale adaptive grid effectively balances the calculation efficiency and result accuracy of the phase field model in solving the crack expansion problem, and improves the universality of the phase field method in engineering application.

[0041] (2) The algorithm complexity of the adaptive grid dynamic adjustment is reduced under the structured grid framework, and the flexibility of the matching between the coarse grid element and the fine grid element is enhanced by using the multi-node transition element.

[0042] (3) The adaptive multi-scale grid construction method is highly similar to the traditional finite element method, which facilitates the construction and solution of the finite element model of the multi-physical field coupling problem. BRIEF DESCRIPTION OF DRAWINGS

[0043] The drawings described herein are used to provide further understanding of the embodiments of the present application, constitute a part of the present application, and do not constitute a limitation of the embodiments of the present application. In the drawings:

[0044] Figure 1 A flowchart of a fracturing crack extension simulation method based on an adaptive multi-scale phase field model is provided for the embodiments of the present application;

[0045] Figure 2 A grid division schematic diagram of a reservoir geometric model is provided for the embodiments of the present application;

[0046] Figure 3 A schematic diagram of a preset initial hydraulic fracture is provided for the embodiments of the present application;

[0047] Figure 4 A schematic diagram of refined and non-refined elements near a local crack is provided for the embodiments of the present application;

[0048] Figure 5 A local crack vicinity refinement unit, a transition unit and a coarsening unit provided for an embodiment of the present application are shown in the schematic diagram.

[0049] Figure 6 A two-dimensional four-node quadrilateral mesh isoparametric transformation provided for an embodiment of the present application is shown in the schematic diagram.

[0050] Figure 7 A multi-scale adaptive mesh size and node distribution provided for an embodiment of the present application is shown in the schematic diagram.

[0051] Figure 8 A mixed finite element mesh at an injection time of 1s provided for an embodiment of the present application is shown in the diagram.

[0052] Figure 9 A displacement distribution along the y direction at an injection time of 1s provided for an embodiment of the present application is shown in the diagram.

[0053] Figure 10 A phase field value distribution at an injection time of 1s provided for an embodiment of the present application is shown in the diagram.

[0054] Figure 11 A mixed finite element mesh at an injection time of 10s provided for an embodiment of the present application is shown in the diagram.

[0055] Figure 12 A displacement distribution along the y direction at an injection time of 10s provided for an embodiment of the present application is shown in the diagram.

[0056] Figure 13 A mixed finite element mesh at an injection time of 10s provided for an embodiment of the present application is shown in the diagram. DETAILED DESCRIPTION

[0057] In order to make the objectives, technical solutions and advantages of the present application clearer, further detailed description will be made to the present application in combination with embodiments and drawings, and the schematic embodiments of the present application and the description thereof are only used to explain the present application, and not as a limitation to the present application.

[0058] It should be noted that the term "include" or "may include" used in various embodiments of the present application indicates the existence of the claimed function, operation or element, and does not limit the addition of one or more functions, operations or elements. In addition, as used in various embodiments of the present application, the terms "include", "have" and their homonyms are only intended to indicate a specific feature, number, step, operation, element, component or combination of the foregoing, and should not be understood as first excluding the existence or addition of one or more other features, numbers, steps, operations, elements, components or combinations of the foregoing.

[0059] It should be understood that terms such as "first", "second" are used only for descriptive purposes and cannot be understood as indicating or implying relative importance or implicitly indicating the number of the technical features indicated. Therefore, the features defined as "first", "second" can be explicitly or implicitly included one or more features. In the description of the present application, the meaning of "a plurality of" is two or more, unless otherwise explicitly specified.

[0060] Reference is made to Figure 1 , Figure 1 A flowchart of a fracturing fracture extension simulation method based on an adaptive multi-scale phase field model is provided for an embodiment of the present application, as shown in Figure 1 The method comprises the following steps:

[0061] S101, obtain the parameters of reservoir properties and rock mechanics, and establish a reservoir geological model according to the parameters.

[0062] Specifically, (1) the parameters of reservoir properties and rock mechanics include the horizontal maximum principal stress, the horizontal minimum principal stress, the pore pressure, the reservoir thickness, the matrix porosity, the matrix permeability, the Young's modulus, the Poisson's ratio, the matrix bulk modulus, the skeleton bulk modulus, and the rock tensile strength. In this embodiment, the simulation basic parameters of reservoir properties and rock mechanics parameters are shown in Table 1.

[0063] (2) According to the collected parameters of reservoir properties and rock mechanics, a hydraulic fracturing reservoir geological model is established through the mapping relationship between spatial position and physical properties. It should be understood that this embodiment adopts the assumption of homogeneous reservoir, that is, the physical property parameters and rock mechanics parameters are the same in the reservoir.

[0064] Table 1 Simulation basic parameters

[0065]

[0066] S102, create a reservoir geometric model based on the reservoir geological model and the preset hydraulic fracturing position attribute, divide the reservoir geometric model into a plurality of finite element grids, and assign values to the grid attributes of the finite element grids using the parameters, to obtain the initial variable values of the preset hydraulic fractures in the finite element grids and the finite element grid nodes.

[0067] Based on the reservoir geological model and the preset hydraulic fracturing position attribute, a reservoir geometric model is created, and the reservoir geometric model is divided into a plurality of finite element grids, specifically:

[0068] According to the hydraulic fracturing position attribute, a rectangular or cuboid region is selected from the reservoir geological model to create a reservoir geometric model;

[0069] The reservoir geometric model is meshed by using finite elements to obtain a calculation domain composed of finite element grids.

[0070] The grid properties of the finite element grid are valued by using parameters, and the preset hydraulic fracture in the finite element grid and the initial variable value on the node of the finite element grid are obtained, specifically:

[0071] The hydraulic fracture is preset in the finite element grid; wherein the hydraulic fracture includes one-dimensional line segment or two-dimensional plane, and the hydraulic fracture does not coincide with the interface of the finite element grid;

[0072] The initial variable value is preset on the grid node of the finite element grid, and the initial variable value includes initial displacement value, initial pressure value and initial phase field value.

[0073] Specifically, for a two-dimensional model, a rectangle is selected as a geometric model, referring to Figure 2 , a coarse square grid is used to divide the geometric model. Wherein the horizontal maximum principal stress direction and the horizontal minimum principal stress direction are along the global coordinate x axis and y axis direction respectively. In the two-dimensional model, referring to Figure 3 , a one-dimensional line segment is used to preset the initial fracture hydraulic fracture, to ensure that the fracture does not coincide with the interface of the grid element. Wherein the initial hydraulic fracture length is 0.15m, along the horizontal maximum principal stress direction, and the injection point is located at the origin (0, 0). The initial variable value of the finite element grid node is preset, wherein the initial displacement value is calculated by the stress boundary condition, the initial pressure value is equal to the reservoir pore pressure value, and the initial phase field value is calculated by the hydraulic fracture trajectory and the historical state variable.

[0074] S103, according to the position relationship between the preset hydraulic fracture and the finite element grid, a multi-scale adaptive grid for dynamic expansion of hydraulic fracture is constructed;

[0075] Specifically, the multi-scale adaptive grid includes the following steps:

[0076] (1) determine the crack tip element and the non-crack tip element by the relative position relationship between the preset hydraulic fracture and the finite element grid;

[0077] (2) take the center of the finite element grid in the main body area of the crack as a circle point, and divide the area with a given critical radius value as the radius and take the union, mark the finite element grid in the area as refinement element, and mark the finite element grid outside the area as non-refinement element;

[0078] (3) search for the non-refinement element adjacent to the refinement element as transition element, and take the remaining non-adjacent non-refinement element as coarsening element; wherein the number of nodes of the refinement element and the coarsening element is the same, the size of the refinement element and the coarsening element has a multiple relationship, and the size of the refinement element matches the length scale parameter of the phase field.

[0079] The refining unit is used to finely depict the complex expansion behavior of the hydraulic fracture, and the coarsening unit exists in the area far from the hydraulic fracture, and a larger coarse grid can be used.

[0080] In the embodiment, the crack tip unit and the crack non-tip unit are determined according to the position of the initial hydraulic fracture, the crack tip unit selects 2 times the size of the coarse grid unit as the critical radius value, the crack non-tip unit selects 1 times the size of the coarse grid unit as the critical radius value, then a circle is drawn with the grid center as the center point, and then the refining unit and the non-refining unit are divided, and the schematic diagram is referred to Figure 4 . According to the division rule, the refining unit, the transition unit and the coarsening unit are divided, and the schematic diagram is referred to Figure 5 . Among them, the refining unit and the coarsening unit adopt four-node square grid units, and the transition unit adopts multi-node square grid units.

[0081] S104, a coupling mathematical model of stress field, seepage field and phase field of hydraulic fracturing is established, and a numerical calculation expression of the coupling mathematical model is constructed according to the multi-scale adaptive grid.

[0082] In the embodiment, the mathematical models of the stress field, the seepage field and the phase field are described respectively, as follows:

[0083] (1) Mathematical model of stress field

[0084] The stress equilibrium equation under the influence of body force and inertial force is not considered, and the reservoir is regarded as an elastic porous medium: .

[0085] According to the effective stress principle of porous medium, the following formula can be obtained: , wherein σ represents the total stress; σ eff represents the effective stress; α represents the Biot coefficient; c represents the phase field value; I represents the unit tensor; p represents the reservoir pore pressure; represents the divergence operator;

[0086] (2) Mathematical model of seepage field

[0087] The fluid flow continuity equation in the saturated porous medium is: .

[0088] The fluid volume increment calculation formula is: .

[0089] The Darcy seepage equation of the porous medium is: .

[0090] The above formula can be obtained: , wherein, ξrepresents fluid volume increment; t represents time; v represents fluid flow rate; ε v represents volume strain; represents Biot modulus; k represents matrix permeability; μ f represents fluid viscosity; represents Biot coefficient; represents fluid viscosity;

[0091] (3) Mathematical model of phase field

[0092] According to the variational principle, the total free energy density equation in the hydraulic fracturing process of porous media is:

[0093]

[0094] The elastic energy density calculation formula is:

[0095] The pore fluid energy density calculation formula is:

[0096] The fracture energy density calculation formula is: , wherein, ψ elas represents elastic energy density; ψ fluid represents pore fluid energy density; ψ frac represents fracture energy density; ψ total represents total free energy density; ε represents solid elastic strain tensor; l 0 represents length scale parameter; G c represents critical energy release rate; represents Lame constant.

[0097] The micro-force balance equation of porous media under the condition of ignoring micro-inertia and external micro-force is:

[0098] The internal energy balance equation in the hydraulic fracture propagation process of porous media under elastic deformation condition is:

[0099] ; wherein; H m represents micro-traction force acting on the fracture, K m represents internal micro-force; ρ represents density; represents internal energy.

[0100] ​​​​The rate of change of the total free energy density is:

[0101] ;

[0102] According to the principle of thermodynamic consistency, we have:

[0103] The above equations can be combined to give: ;

[0104] ;

[0105] , where g (c) represents the decay function; represents the tensile stress; represents the compressive stress; H frac represents the history state variable.

[0106] (4) Evolution model of pore elastic parameters

[0107] The mathematical model of porosity evolution in porous media is: .

[0108] The mathematical model of Biot coefficient and Biot modulus evolution is:

[0109] , the expression of the decay function g (c) is: , where ϕ 0 represents the porosity of the matrix before damage occurs; K , K s , K f represent the bulk modulus of the porous medium, the bulk modulus of the solid skeleton, and the bulk modulus of the liquid, respectively; χ represents the minimum value, used to avoid numerical singularity;

[0110] (5) Boundary conditions of control equations

[0111] The boundary conditions of hydraulic fracture propagation in porous media include:

[0112] ;

[0113] ;

[0114] ; where ∂Ω u , ∂Ω u , Γ f represent the Dirichlet boundary of the stress field, seepage field, and phase field, respectively; ∂Ω t∂Ω q ∂Ω and ∂Ω represent the Neumann boundaries of the stress field, seepage field, and phase field, respectively; u represents the displacement vector; and The vectors represent displacement and pressure at the boundary; t and q represent the boundary stress vector and flow rate, respectively; n represents the normal vector.

[0115] (6) Numerical calculation expression

[0116] First, the equivalent weak integral form of the governing equations for the stress field, seepage field, and phase field under a multi-scale adaptive grid is constructed using the finite element method, as follows:

[0117] ;

[0118] ;

[0119] In the formula, w u , w p , w c These represent the weighting functions for displacement, pressure, and phase field, respectively, with the superscript T indicating transpose.

[0120] The time term is discretized using a backward Euler difference scheme, and then the equivalent weak integral form of the governing equations is linearized by matrix transformation to obtain the numerical expression, as follows:

[0121] ;

[0122] ;

[0123] ;

[0124] ;

[0125] ;

[0126] ;

[0127] In the above formula, the superscript n represents the value at the element node; N u N p N c B represents the shape functions representing the differences between displacement, pressure, and phase field, respectively; u , B represents the strain matrix and the volumetric strain matrix, respectively; p B c R represents the shape function gradient of the difference between pressure and phase field, respectively; u R p R crespectively represent the displacement, pressure and phase field; δu, δp, δc respectively represent the increment of displacement, pressure and phase field; subscript (t) represents the last time step; t -Δ t ) represents the last time step;

[0128] The numerical calculation expression is solved by numerical iteration, specifically: preset initial values of displacement, pressure and phase field; the initial values are taken as input, the Newton-Raphson iteration method is used to iteratively solve the control equations of the stress field, the seepage field and the phase field, and the trial solution of each iteration step is obtained; if the trial solutions of adjacent iteration steps at the current time step satisfy the convergence condition, the solution converges, and the displacement, pressure and phase field values of the node are output.

[0129] Specifically, the numerical base function, local stiffness matrix and local equivalent load matrix of the transition unit are constructed.

[0130] Specifically, taking a single transition unit as the construction object, the local fine mesh and coarse mesh are divided, and the solid phase stiffness matrix of the stress field, the seepage matrix of the seepage field, and the phase evolution matrix are taken as the left end items to construct the local balance equation;

[0131] The boundary number and node number of the transition unit are determined by the connection condition of the transition unit and the adjacent coarse unit and fine unit, and the boundary of the transition unit sharing the same node is taken as a group to apply linear boundary conditions;

[0132] Combining the local balance equation of the transition unit and the linear boundary condition, the numerical base function of the transition unit in displacement, pressure and phase field is sequentially solved;

[0133] The conversion relationship between the finite element stiffness matrix, the equivalent load matrix and the node value between the coarse mesh and the fine mesh is established through the numerical base function. It should be noted that the matrix form of the control equation of the transition unit is:

[0134]

[0135] Wherein,

[0136] ;

[0137] ;

[0138] ;

[0139] ;

[0140] ;

[0141] ;

[0142] ;

[0143] ;

[0144] The conversion relationship between the local fine mesh unit and the coarse mesh unit of the transition unit is:

[0145] ;

[0146] Wherein, ;

[0147] ;

[0148] ;

[0149] ;

[0150] ;

[0151] ; In the formula, the superscripts C and F represent the local coarse mesh unit and the fine mesh unit, and the subscripts u, p and c represent the displacement field, the pressure field and the phase field respectively; K and F represent the stiffness matrix and the equivalent load matrix respectively; , , represent the displacement numerical base function, the pressure numerical base function and the phase field numerical base function respectively; u f , p f , c f represent the displacement, the pressure and the phase field matrix of the local fine mesh node; u c , p c , c c represent the displacement, the pressure and the phase field matrix of the local coarse mesh node.

[0152] The local stiffness matrix and the local equivalent load matrix of the stress field, the seepage field and the phase field of the refined unit and the coarsened unit are calculated, the global stiffness matrix and the equivalent load matrix of the multi-scale adaptive grid are assembled by combining the local stiffness matrix and the local equivalent load matrix of the transition unit, and the global stiffness matrix and the equivalent load matrix constitute the control equation;

[0153] The numerical solution of the coupling model is realized by combining the numerical iteration mode of the mathematical model and the preset initial variable value, and the node value conversion between the local fine mesh and the coarse mesh of the transition unit is realized by the numerical base function.

[0154] Specifically, the control of the stress field, the seepage field and the phase field is solved by presetting the initial values of displacement, pressure and phase field, using Newton-Raphson iteration method, and using Picard iteration to enhance the stability of numerical solution, and the firsti The trial solution of the +1 iteration step is: ; wherein, represents the retardation coefficient of the Picard iteration; , , displacement matrix, pressure matrix and phase field matrix of the n-th iteration step, respectively.

[0155] When the displacement, pressure and phase field values of the nodes of the adjacent iteration step at the current time step satisfy the following convergence condition, the numerical iteration solution converges, and the displacement, pressure and phase field values of the nodes are output:

[0156] , wherein, , , represent the iteration convergence tolerances of the displacement, pressure and phase field, respectively.

[0157] In the present embodiment, a two-dimensional adaptive multi-scale phase field model is used to simulate the dynamic extension of hydraulic fractures, a structured quadrilateral grid is used as a calculation unit, numerical solution is carried out based on finite element isoparametric transformation, a two-dimensional 4-node quadrilateral element isoparametric transformation diagram is shown in Figure 6 , and the convergence tolerances of the displacement, pressure and phase field are 1.0x10 -7 , 1.0x10 -5 and 1.0x10 -5 , respectively.

[0158] S105, based on the numerical calculation expression, using the parameters of reservoir properties and rock mechanics and the boundary grid node variable values of the multi-scale adaptive grid as constraints, solving the coupled mathematical model to obtain the simulation results of the hydraulic fracture extension under different reservoir and construction conditions.

[0159] Specifically, the way to solve the coupled mathematical model is as follows:

[0160] S1051, taking the preset initial variable values on the finite element grid nodes as input, combining the reservoir properties and rock mechanics parameters and the boundary conditions, carrying out numerical iteration solution on the numerical calculation expression to obtain the trial solution of the displacement, pressure and phase field values of the grid nodes at the current time step;

[0161] S1052, after obtaining the asymptotic solution of the phase field values of the grid nodes at the current time step, dividing the crack tip elements and the crack non-tip elements according to the phase field values of the grid elements, carrying out element refinement operation and element coarsening operation using different critical radius values, adaptively adjusting the global grid and node number to realize adaptive division of the multi-scale grid in the process of dynamic extension of the hydraulic fracture;

[0162] S1053, if the trial solution of the adjacent iteration step at the current time step satisfies the convergence condition, the solution converges, and the displacement value, pressure value and phase field value of the node are output, and the calculation of the next time step is started;

[0163] S1054, the loop is calculated until the simulation termination condition is met, and the displacement, pressure and phase field calculation results of the multi-scale adaptive grid nodes are output and visualized, and the simulation results of the hydraulic fracture dynamic extension under different reservoir and construction conditions are obtained.

[0164] In this embodiment, the number of nodes of the coarse grid element and the fine grid element in the two-dimensional model is 4, and the length of the coarse grid element is 5 times the length of the fine grid element. The number of nodes of the transition element is different due to the different sizes of the associated coarse grid element and fine grid element. The schematic diagram is referred to Figure 7 When 3 of the 4 edges of the transition element are adjacent to the coarse grid element and 1 edge is adjacent to the fine grid element, the total number of nodes of the transition element is 8. Similarly, when 2 of the 4 edges of the transition element are adjacent to the coarse grid element, the total number of nodes of the transition element is 12. Through the change of the nodes of the transition element, the coarse grid element and the fine grid element can be adaptively combined to realize multi-scale adaptive grid partitioning of the calculation domain. According to the different number and distribution of nodes of the transition element, different displacement, pressure and field numerical basis functions are constructed using linear boundary conditions, and further, the global stiffness matrix and equivalent load matrix are constructed through the local stiffness matrix and equivalent load matrix of the fine grid element, the transition element and the coarse grid element. The displacement, pressure and phase field values of the multi-scale adaptive grid nodes at different times are obtained by solving the coupled mathematical model through numerical iteration. After obtaining the progressive solution of the phase field value at the current time step, the crack tip element and the non-crack tip element are further divided, and the global grid and node number are adaptively adjusted. In this embodiment, for the crack tip element, the critical radius value of element refinement is 5 times the length of the coarse grid element; and for the non-crack tip element, the critical radius value of element refinement is 2 times the length of the coarse grid element.

[0165] The hybrid finite element grid, the displacement value along the y direction and the phase field value at the injection time of 1s are respectively referred to Figure 8 , Figure 9 , Figure 10 . Among them, through the results of the hybrid finite element grid, we can observe the multi-scale adaptive grid partitioning method of the calculation domain, and the grid refinement operation is performed in the circular area with a radius of 5 times the length of the coarse grid element at the hydraulic fracture tip. Figure 9 is the displacement value along the y direction, which is obtained by stress field calculation. Figure 10 is the phase field value, and the area where the phase field value is not 0 represents the hydraulic fracture extension trajectory. The hybrid finite element grid, the displacement value along the y direction and the phase field value at the injection time of 10s are respectively referred to Figure 11 ,Figure 12 , Figure 13 Compared with the results of 1s, the adaptive division of multi-scale grid in the dynamic extension process of hydraulic fracture is realized, and the simulation of the dynamic extension of hydraulic fracture is completed. In the simulation process, the adaptive division of coarse grid cells and fine grid cells is realized, the calculation amount of the phase field model is effectively reduced, and the extension trajectory of the hydraulic fracture at different times is obtained under the given reservoir and construction parameters.

[0166] The above specific embodiments further illustrate the purpose, technical solutions and beneficial effects of the present application. It should be understood that the above description is only a specific embodiment of the present application and is not intended to limit the protection scope of the present application. Any modification, equivalent replacement, improvement, etc. within the spirit and principles of the present application should be included in the protection scope of the present application.

Claims

1. A method for simulation of fracture extension based on an adaptive multi-scale phase-field model, the method comprising: The method comprises the following steps: ​ obtaining parameters of reservoir properties and rock mechanics, and establishing a reservoir geological model according to the parameters; creating a reservoir geometric model based on the reservoir geological model and a preset hydraulic fracturing position attribute, dividing the reservoir geometric model into a plurality of finite element grids, and assigning values to grid attributes of the finite element grids by using the parameters to obtain preset hydraulic fractures in the finite element grids and initial variable values on nodes of the finite element grids; constructing a multi-scale self-adaptive grid for dynamic expansion of the hydraulic fractures according to the position relationship between the preset hydraulic fractures and the finite element grids; wherein the multi-scale self-adaptive grid for dynamic expansion of the hydraulic fractures is constructed according to the position relationship between the preset hydraulic fractures and the finite element grids, which comprises the following steps: determining crack tip elements and non-crack tip elements through the relative position relationship between the preset hydraulic fractures and the finite element grids; taking the center of the finite element grid in the main body region of the crack as a circle point, and performing regional division and set operation by taking a pre-configured critical radius value as a radius, and identifying the finite element grids in the region as refined elements, and identifying the finite element grids outside the region as non-refined elements; searching for non-refined elements adjacent to the refined elements as transition elements, and taking the remaining non-adjacent non-refined elements as coarsened elements; wherein the number of nodes of the refined elements and the coarsened elements is the same, there is a multiple relationship between the sizes of the refined elements and the coarsened elements, and the size of the refined elements matches the length scale parameter of the phase field; establishing a coupled mathematical model of stress field, seepage field and phase field of hydraulic fracturing, and constructing a numerical calculation expression of the coupled mathematical model according to the multi-scale self-adaptive grid; wherein the coupled mathematical model of stress field, seepage field and phase field of hydraulic fracturing is constructed, and the numerical calculation expression of the coupled mathematical model is constructed according to the multi-scale self-adaptive grid, which specifically comprises the following steps: constructing a numerical base function, a local stiffness matrix and a local equivalent load matrix of the transition elements; calculating the local stiffness matrix and the local equivalent load matrix of the stress field, the seepage field and the phase field of the refined elements and the coarsened elements, combining the local stiffness matrix and the local equivalent load matrix of the transition elements to assemble a global stiffness matrix and an equivalent load matrix of the multi-scale self-adaptive grid, and composing control equations of the stress field, the seepage field and the phase field according to the global stiffness matrix and the equivalent load matrix; realizing the conversion of node values between the local fine grid and the coarse grid of the transition elements by using the numerical base function, and iteratively solving the control equations of the stress field, the seepage field and the phase field by using the Newton-Raphson iteration method to obtain a progressive solution of each iteration step, and realizing the numerical solution of the coupled mathematical model; if the progressive solutions of adjacent iteration steps at the current time step satisfy the convergence condition, the solution result of the numerical solution converges, and the displacement value, the pressure value and the phase field value of the node are outputted; based on the numerical calculation expression, taking the parameters of reservoir properties and rock mechanics and the boundary grid node variable values of the multi-scale self-adaptive grid as constraints, solving the coupled mathematical model to obtain hydraulic fracturing crack extension simulation results under different reservoir and construction conditions.

2. The method of claim 1, wherein, The parameters include a horizontal maximum principal stress, a horizontal minimum principal stress, a reservoir pore pressure, a reservoir thickness, a matrix porosity, a matrix permeability, a Young's modulus, a Poisson's ratio, a matrix bulk modulus, a skeleton bulk modulus, and a rock tensile strength.

3. The method of claim 1, wherein, The reservoir geological model is established according to the parameters, specifically by establishing the reservoir geological model according to a mapping relationship between spatial positions and physical properties of the parameters.

4. The method of claim 1, wherein, The reservoir geometry model is created based on the reservoir geological model and a preset hydraulic fracturing position attribute, and the reservoir geometry model is divided into a plurality of finite element grids, specifically by: selecting a rectangular or cuboid region from the reservoir geological model according to the hydraulic fracturing position attribute to create the reservoir geometry model; and dividing the reservoir geometry model into grids by using finite elements to obtain a calculation domain composed of the finite element grids.

5. The method of claim 4, wherein, The grid attributes of the finite element grids are valued by using the parameters to obtain preset hydraulic fractures in the finite element grids and initial variable values on nodes of the finite element grids, including: presetting the hydraulic fractures in the finite element grids; wherein the hydraulic fractures include one-dimensional line segments or two-dimensional planes, and the hydraulic fractures do not coincide with interfaces of the finite element grids; presetting initial variable values on nodes of the finite element grids; wherein the initial variable values include initial displacement values, initial pressure values, and initial phase field values.

6. The method of claim 5, wherein, The initial displacement values are calculated by applying stress boundary conditions, the initial pressure values are equal to reservoir pore pressure values, and the initial phase field values are calculated by using hydraulic fracture trajectories and historical state variables.

7. The method of claim 1, wherein, The numerical base functions, local stiffness matrices, and local equivalent load matrices of the transition elements are constructed, including: taking a single transition element as a construction object, locally dividing fine grids and coarse grids, and constructing local balance equations by taking a solid phase stiffness matrix of a stress field, a permeation matrix of a seepage field, and a phase field evolution matrix as left end terms; determining the number of boundaries and nodes of the transition element according to a connection condition of the transition element with adjacent coarse elements and fine elements, and applying linear boundary conditions to boundaries of the transition element sharing a same node as a group; sequentially solving the numerical base functions of the transition element in displacement, pressure, and phase field by combining the local balance equations of the transition element with the linear boundary conditions; establishing conversion relationships of finite element stiffness matrices, equivalent load matrices, and node values between the coarse grids and the fine grids by using the numerical base functions.

8. The method of claim 1, wherein, The coupled mathematical model is solved based on the numerical calculation expression and by taking the parameters of reservoir properties and rock mechanics and the variable values of boundary grid nodes of the multi-scale adaptive grids as constraints to obtain simulation results of hydraulic fracturing fracture extension under different reservoir and construction conditions, specifically by: taking the preset initial variable values on the nodes of the finite element grids as inputs, combining the parameters of reservoir properties and rock mechanics and the variable values of the boundary grid nodes, and numerically iteratively solving the numerical calculation expression to obtain progressive solutions of displacement, pressure, and phase field values of the nodes at a current time step. After obtaining the asymptotic solution of the phase field value of the grid node at the current time step, the fracture tip element and the non-tip element of the fracture are divided according to the phase field value of the grid element, different critical radius values are used for element refinement and element coarsening operations, and the global grid and node number are adaptively adjusted to realize the adaptive division of the multi-scale grid in the dynamic expansion process of the hydraulic fracture; If the asymptotic solution of the adjacent iteration step at the current time step satisfies the convergence condition, the solution converges, the displacement value, the pressure value and the phase field value of the node are output, and the calculation of the next time step is started; The loop calculation is performed until the simulation termination condition is met, the calculation results of the displacement, pressure and phase field of the multi-scale adaptive grid node are output and data visualization processing is performed, and the dynamic extension simulation results of the hydraulic fracture under different reservoirs and construction conditions are obtained.

Images