Rock crack propagation simulation method under pulse fracturing effect

Abstract

The invention relates to the technical field of rock damage simulation, in particular to a rock crack propagation simulation method under the action of pulse fracturing, which comprises the following steps: introducing a phase field variable, and constructing a phase field method model; establishing a friction damage physical mechanism based on a rate state dependent friction law and a sliding law; introducing a friction damage physical mechanism into the phase field method model, establishing a coupling relationship between a friction state variable and the phase field method model, and constructing a phase field-friction damage coupling control equation set; constructing a rock loading model under the action of pulse fracturing; solving a phase field-friction damage coupling control equation set by adopting an adaptive alternating iteration algorithm introducing a dynamic local grid; a friction damage physical mechanism is introduced into a phase field method model, a rock crack network under the action of pulse fracturing can be simulated more truly and accurately, a dynamic local grid adaptive method is adopted, and waste of computing resources is reduced while rock crack propagation simulation precision is guaranteed.

Description

Technical Field

[0001] This invention relates to the field of rock damage simulation technology, and specifically to a method for simulating rock crack propagation under pulsed fracturing. Background Technology

[0002] Pulse fracturing is a novel fracturing technology that uses a cyclic loading process of "pressurization-stabilization-decompression" to modify reservoirs. Its core mechanism is to utilize the fatigue effect of rocks to induce the initiation, propagation and connection of microcracks under cyclic stresses that are far below the static strength of the rocks, thereby forming a more complex fracture network and increasing the reservoir modification volume.

[0003] When designing reservoir stimulation projects, the key lies in simulating rock crack propagation and predicting the morphology of the rock fracture network formed by fracturing in order to optimize fracturing parameters. Existing rock crack propagation simulation methods mainly target the fracturing process under continuous loading. They usually assume that cracks only propagate when the load increases or reaches its peak, ignoring the unloading stage. The relative slippage and friction between pre-existing microfracture surfaces inside the rock also consume energy, cause damage, and may drive the physical mechanism of new crack initiation. This makes it impossible for existing rock crack propagation simulation methods to capture the cumulative effect of damage and the gradual propagation of cracks under pulsed loading, and they cannot simulate the fracture network under pulsed fracturing.

[0004] Furthermore, due to the large scale and random distribution of rock fracture networks under pulsed fracturing, extremely high spatial resolution is required in simulations to capture the initiation and fine propagation paths of microcracks. Most existing rock crack propagation simulation methods employ either a globally uniform fine mesh or a pre-defined local fine mesh. A globally uniform fine mesh results in a huge number of degrees of freedom and high computational costs, making it difficult to apply to three-dimensional reservoir-scale simulations in actual reservoir stimulation engineering design. On the other hand, the method of pre-defined local fine mesh cannot predict the specific location and propagation direction of a large number of randomly initiating and dynamically evolving microcracks, thus failing to achieve adaptive mesh refinement, leading to either wasted computational resources or insufficient resolution in critical areas.

[0005] Therefore, it is necessary to provide a simulation method for rock crack propagation under pulsed fracturing that can match the physical mechanism of rock crack propagation under pulsed fracturing, and simulate the rock fracture network under pulsed fracturing with high efficiency and high accuracy at an affordable computational cost, so as to provide strong support for the design of reservoir stimulation projects using pulsed fracturing technology. Summary of the Invention

[0006] To address the aforementioned problems, this invention provides a method and application for simulating rock crack propagation under pulsed fracturing, thereby resolving the technical issues raised in the background section.

[0007] The technical solution adopted by this invention to solve its technical problem is as follows: A method for simulating rock crack propagation under pulsed fracturing includes the following steps: S1: Introduce phase field variables to describe the fracture state of rock materials, combine them with the displacement field of the rock, establish a total energy functional including elastic strain energy density and fracture energy, and construct a phase field method model. S2: Introducing frictional state variables, based on the rate-state dependent frictional law and the sliding law, the evolution law of frictional damage under pulse fracturing is determined, and the physical mechanism of frictional damage regarding the total damage variable is established. S3: Introduce the physical mechanism of friction damage into the phase-field method model, and establish the coupling relationship between friction state variables and elastic strain energy density and fracture energy based on the stiffness and strength properties of rock materials, respectively, and construct the phase-field-friction damage coupled control equation set. S4: Establish a basic rock geological model based on the geological conditions of the reservoir to be modified, set relevant parameters such as rock material parameters, pulse fracturing load conditions, boundary conditions, and friction damage physical mechanism, and construct a rock loading model under pulse fracturing. S5: For the rock loading model, an alternating iterative algorithm with dynamic local mesh adaptation is adopted to decouple the displacement field, friction state variables, total damage variables, and phase field variables, solve the phase field-friction damage coupled control equations, and output the simulation results of rock crack propagation.

[0008] Furthermore, the total energy functional Π established in S1 is: ; In the formula, u is the displacement field. For phase field variables, Indicates that the rock material is intact. This indicates that the rock material has completely fractured. This indicates that the rock material is in the process of fracture. For elastic strain energy density, Let γ be the fracture energy, Ω be the phase field surface density function, Ω be the rock solid domain, and ε be the strain field; a degradation function is introduced into the elastic strain energy density. To reduce the elastic strain energy of the damage, ,in, It represents tensile strain energy.

[0009] Furthermore, the specific method of S2 is as follows: Based on the rate-state dependent friction law, the relationship between the friction coefficient μ, the equivalent sliding speed V, and the friction state variable θ is constructed as follows: ; where \(a\) is the direct rate effect parameter, \(a > 0\), \(b\) is the state evolution effect parameter, \(b > a\), and and are the reference friction coefficient, reference sliding rate and reference friction state variable respectively; Based on the sliding law, establish the evolution law relation of the friction state variable \(\theta\): ; where \(L\) is the characteristic sliding distance. When \(V\approx0\), it means the current position is in the "healing" or "aging" state. When \(V > 0\), it means the current position is being updated; It is set that the total damage variable \(D\) under pulsed fracturing is obtained by gradually accumulating the friction damage increment in time. The friction damage increment follows the rate-state dependent friction law, that is, ; where \(\lambda\) is the damage accumulation coefficient, \(0 < \lambda < 1\).

[0010] Furthermore, the specific method for establishing the coupling relationship between the friction state variable and the elastic strain energy density in step S3 is as follows: Modify the degradation function to a bivariate function related to the phase field variable and the friction state variable \(\theta\), ; where \(q\) is a regularization constant with a value of , \(h(\theta)\) is the friction damage factor, , \(0 < h(\theta)\leq1\), \(\alpha\) is the coupling coefficient, and \(( )\) + means taking the positive value; establish the stiffness degradation coupling relation: , \(C\) is the elastic stiffness tensor, is the initial stiffness tensor in the intact state of the rock material; establish the coupling relationship between the friction state variable and the elastic strain energy density, that is .

[0011] Furthermore, the specific method for establishing the coupling relationship between the friction state variable and the fracture energy in step S3 is as follows: Modify the fracture energy to a function related to the friction state variable \(\theta\), ; where is the initial fracture energy, \(\beta\) is the wear coefficient, \(0\leq\beta<1\).

[0012] Furthermore, the specific steps of step S5 are as follows: S51: Spatial discretization of the rock loading model is performed, dividing it into a uniform coarse grid. Temporal discretization of the phase field-friction damage coupled control equation set is performed, dividing it into multiple time steps. S52: At the current time step, using the current mesh of the rock loading model, the alternating iterative solution algorithm is used to iteratively solve and update the displacement field, friction state variables, total damage variables, and phase field variables until convergence. S53: Determine whether the current time step meets the grid update conditions. If yes, adaptively refine and coarse the current grid and update the current grid. If not, go to S55. S54: Map the displacement field, friction state variables, total damage variables, and phase field variables from the previous mesh to the current mesh using linear interpolation; S55: Determine whether all time steps have been solved. If yes, output the simulation results of rock crack propagation. If not, return to S52.

[0013] Furthermore, the specific method for iteratively solving and updating the displacement field, friction state variables, total damage variables, and phase field variables in S52 is as follows: Freeze the phase field variables, friction state variables, and total damage variables obtained in the previous iteration, construct the reduced stiffness matrix using the stiffness degradation coupling relationship, substitute it into the discretized mechanical equilibrium equations, and solve for the displacement field in this iteration. The mechanical equilibrium equations are as follows: Where σ is the stress tensor, , For volume forces, satisfying the following conditions at the boundary of the rock solid domain: , , Given the displacement boundary function; The equivalent sliding rate is updated using the displacement field obtained from this iteration. The friction state variables for this iteration are solved based on the velocity-state-dependent friction law. The total damage variable at the current time step is then calculated and updated. Based on establishing the coupling relationship between friction state variables and fracture energy, the reduced fracture energy is updated using the friction state variables and total damage variables obtained in this iteration. The energy is then substituted into the discretized phase field evolution equation to solve for the phase field variables in this iteration. The phase field evolution equation is obtained through total energy functional variation.

[0014] Furthermore, the specific method for adaptively refining and coarsening the current mesh in S53 is as follows: Calculate the phase field gradient magnitude of each cell in the current mesh. ; Set the first threshold With the second threshold , The grid cells satisfy Marked as crack zone elements, the mesh elements satisfy... Marked as process area element, the mesh element satisfies Marked as a complete region unit; The elements marked as crack zones and their adjacent elements are refined; the elements marked as intact zones are examined for their relationship with adjacent elements, and these elements are merged to form coarse mesh elements while satisfying geometric constraints and numerical accuracy; the elements marked as process zones are set to meshes with gradually changing sizes.

[0015] Furthermore, the specific method for establishing the basic rock geological model in S4 is as follows: based on the geological conditions of the reservoir to be modified, the three-dimensional spatial structure of the reservoir to be modified is extracted, and a corresponding basic rock geological model is established. A main fracture of a certain length is preset in the center of the basic rock geological model to simulate the initial fracture formed by the initiation of the perforation cluster. Several natural weak surfaces are randomly distributed around the main fracture, with a length of 2.5%-15% of the main fracture. The dip angle of the randomly distributed natural weak surfaces is reduced by a reduction factor, which is in the range of 0 to 1.

[0016] Furthermore, the specific method for setting the pulse fracturing load conditions and boundary conditions in S4 is as follows: a pulsed fluid pressure load p(t) is applied to the main fracture surface. ,in, Based on pressure, denoted as pressure amplitude and f as pulse frequency; a constant minimum horizontal principal stress is applied around the rock geological foundation model, with vertical displacement constrained at the bottom and horizontal displacement constrained on both sides.

[0017] Compared with the prior art, the beneficial effects of the present invention are: 1. The present invention provides a method for simulating rock crack propagation under pulsed fracturing. It uses the physical mechanism of frictional damage to describe the physical process of rock crack propagation under pulsed fracturing and introduces this physical mechanism into the phase field method model. Compared with existing crack propagation simulation methods, it can more realistically and accurately simulate the rock fracture network under different pulsed fracturing load conditions, which is beneficial for optimizing pulsed fracturing parameters when designing reservoir stimulation projects.

[0018] 2. The present invention provides a method for simulating rock crack propagation under pulsed fracturing. In the simulation solution, a dynamic local mesh adaptive method is adopted, which can automatically and accurately track the crack propagation front. The mesh is refined in the region of crack propagation front to ensure the simulation accuracy of the rock crack network, and a coarse mesh is used in the region far from the crack propagation front to reduce the computational cost of the simulation. While ensuring the simulation accuracy of rock crack propagation, the waste of computing resources is reduced. Attached Figure Description

[0019] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, for those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0020] Figure 1 This is a schematic diagram of the process for simulating rock crack propagation under pulsed fracturing as described in this embodiment; Figure 2 for Figure 1 A flowchart of the S5 process; Figure 3 This is a schematic diagram of the simulation results of rock crack propagation described in this embodiment; Figure 4 This is a schematic diagram of the damage evolution results described in this embodiment; Figure 5 This is a schematic diagram of the pressure-crack volume curve described in this embodiment. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] This invention provides a method for simulating rock crack propagation under pulsed fracturing, such as... Figure 1 As shown, the method includes the following steps: S1: Introduce phase field variables to describe the fracture state of rock materials, combine them with the displacement field of the rock, establish a total energy functional including elastic strain energy density and fracture energy, and construct a phase field method model. Furthermore, the established total energy functional Π is: ; In the formula, the phase field variable As a continuous order parameter field in the phase-field method, it is used to describe the fracture state of rock materials. Indicates that the rock material is intact. This indicates that the rock material has completely fractured. This indicates that the rock material is in a fracture process, where x and t are the spatial coordinate parameters and the time parameter, respectively; u is the displacement field. For elastic strain energy density, The fracture energy is the critical energy release rate of the rock material, γ is the phase field surface density function, Ω is the rock solid domain, and ε is the strain field. A degradation function is introduced into the elastic strain energy density. To reduce the elastic strain energy of the damage, that is ,in, It represents tensile strain energy.

[0023] S2: Based on the rate-state dependent friction law and the sliding law, define the friction state variables, determine the evolution law of friction damage under pulse fracturing, and establish the physical mechanism of friction damage regarding the total damage variable; Furthermore, this invention attributes the fatigue damage evolution of rocks under pulsed fracturing to a frictional damage physical mechanism, which is a frictional physical process controlled by the sliding rate and frictional state that occurs during the sliding process at the crack surface or in the diffuse crack process zone. The specific method of S2 is as follows: Based on the rate-state dependent friction law, a relationship is constructed between the friction coefficient μ, the equivalent sliding speed V, and the friction state variable θ: ; In the formula, a is the direct rate effect parameter, a>0, and b is the state evolution effect parameter, b>a, used to represent the velocity weakening during unstable sliding, and the equivalent sliding rate V is related to the local shear strain rate. Related, , l C It is a phase field length scale parameter, related to the size of the crack surface or the region of diffuse cracking. The frictional state variable has a time dimension and is used to characterize the maturity at various locations of the contact interface where the frictional slip process occurs. , , These represent the reference friction coefficient, reference sliding speed, and reference friction state variable, respectively, with values ​​ranging from 0.6 to 0.85. , ; Based on the law of sliding, the evolution law relationship of the friction state variable θ is established: ; In the formula, L is the characteristic sliding distance. When the current position of the contact interface is stationary, that is, V≈0, dθ / dt≈1, the friction state variable θ increases linearly with time, indicating that the current position is in a "healing" or "aging" state, and the strength of the current position is restored. When the current position of the contact interface slides, that is, V>0, the friction state variable θ evolves towards the steady-state value L / V, indicating that the current position is being updated. The total damage variable D under pulse fracturing is defined as the frictional damage increment in each load cycle. It is obtained by gradually accumulating over time, and then the damage evolution law under pulsed fracturing is determined; according to the rate-state dependent friction law, the friction damage increment is related to the equivalent sliding rate V and the friction state variable θ, ; In the formula, λ is the damage accumulation coefficient, 0 < λ < 1; the total loss variable is used to reduce the effective elastic modulus of the rock material, that is, , is the reduced effective elastic modulus of the rock material, and E is the complete elastic modulus of the rock material when it is not damaged.

[0024] S3: Introduce the physical mechanism of friction damage into the phase field method model, and establish the coupling relationship between the friction state variable and the elastic strain energy density and fracture energy respectively based on the stiffness property and strength property of the rock material, and construct the phase field-friction damage coupling control equations; Furthermore, the specific method for establishing the coupling relationship between the friction state variable and the elastic strain energy density is: modify the degradation function into a bivariate function related to the phase field variable and the friction state variable θ, ; In the formula, q is a regularization constant, and its value is , which is used to ensure the stability of numerical solution, h(θ) is the friction damage factor, , 0 < h(θ) ≤ 1, α is the coupling coefficient, ( ) + means taking the positive value. The coupling relationship between the friction state variable and the elastic strain energy density shows that when the contact interface is updated due to sliding, that is, when the friction state variable decreases, the friction damage factor decreases accordingly, resulting in additional and irreversible degradation of the material stiffness, thus reflecting the cumulative damage caused by frictional slip; establish the stiffness degradation coupling relationship: , C is the elastic stiffness tensor, is the initial stiffness tensor of the rock material in the complete state; according to the relationship between the degradation function and the elastic strain energy density, establish the coupling relationship between the friction state variable and the elastic strain energy density, that is .

[0025] Furthermore, the specific method for establishing the coupling relationship between the friction state variable and the fracture energy is: modify the fracture energy into a function related to the friction state variable θ, ; In the formula, Let β be the initial fracture energy and β be the wear coefficient, 0≤β<1. The coupling relationship between friction state variables and fracture energy indicates that as frictional sliding proceeds under cyclic loading, the friction damage factor decreases, the crack surface is worn and deteriorated, and the fracture energy required to generate a new surface decreases accordingly, enabling the crack to propagate even under lower energy drive.

[0026] Furthermore, by combining the coupling relationship between frictional state variables and elastic strain energy density, and the coupling relationship between frictional state variables and fracture energy, with the total energy functional established in the phase-field method model, we obtain the phase-field-friction damage coupled control equation set.

[0027] S4: Establish a basic rock geological model based on the geological conditions of the reservoir to be modified, set relevant parameters such as rock material parameters, pulse fracturing load conditions, boundary conditions, and friction damage physical mechanism, and construct a rock loading model under pulse fracturing. Furthermore, the specific method for constructing the rock loading model is as follows: Geological conditions of the reservoir to be modified are obtained, the three-dimensional spatial structure of the reservoir to be modified is extracted, and a basic rock geological model is constructed by scaling the three-dimensional spatial structure. The basic rock geological model is used as the computational domain. In this embodiment, the scaling ratio of the basic rock geological model to the three-dimensional spatial structure of the reservoir to be modified is 1:2.5. A main fracture of a certain length is preset at the center of the basic rock geological model to simulate the initial fracture formed by the initiation of the perforation cluster. Several natural weak surfaces are randomly distributed around the main fracture, with a length of 2.5%-15% of the main fracture. The dip angle of the natural weak surfaces is randomly distributed. The normal stiffness and tangential stiffness of the natural weak surfaces are weakened by a reduction factor, with the reduction factor ranging from 0 to 1. For example, if the normal reduction factor is 0.1, the normal stiffness is 0.1 times the stiffness of the intact rock. In this embodiment, the basic rock geological model is a scaled-down shale gas field with an overall size of 100m. The main fracture is set in the horizontal direction with a length of 20m, and there are 50 natural weak surfaces with lengths of 0.5-3m. The rock material parameters are set, including the elastic modulus, Poisson's ratio, and density of the rock matrix, as well as the initial fracture energy of the rock. In this embodiment, the laboratory core mechanical test report of the selected shale gas field shows that the elastic modulus of the rock material in this shale gas field ranges from 25 to 35 GPa, and the Poisson's ratio ranges from 0.2 to 0.3. Therefore, this embodiment sets the elastic modulus to 30 GPa, the Poisson's ratio to 0.25, and the density to 2500 kg / m³. 3 The initial fracture energy of the rock was set to 200 J / m. 2 ; The pulse fracturing load condition is set as applying a pulsed fluid pressure load p(t) on the main fracture surface. ,in, Based on pressure, f is the pressure amplitude, and f is the pulse frequency. In this embodiment, the monitoring results of the selected shale gas field pulse fracturing test show that the pressure fluctuation frequency is mainly distributed in the range of 0.5Hz to 2Hz, and the pressure amplitude is about 15% to 25% of the base pressure. Therefore, in this embodiment, the pulse frequency is set to f=1Hz, the base pressure is 18MPa, the pressure amplitude is 14MPa, and the total application time of the pulsed fluid pressure load is 10 pulse cycles, that is, 10s. The boundary condition is set as applying a constant minimum horizontal principal stress around the rock geological foundation model to simulate the confining pressure of the formation. Vertical displacement is constrained at the bottom of the rock geological foundation model, and horizontal displacement is constrained on both sides. In this embodiment, the minimum horizontal principal stress gradient of the selected shale gas field is about 0.018 MPa / m, and the stress corresponding to a burial depth of 2000m is about 36 MPa. Therefore, based on the scaling ratio of the rock geological foundation model relative to the shale gas field, this embodiment sets the minimum horizontal principal stress to 15 MPa.

[0028] Parameters related to the physical mechanism of frictional damage are defined, including static damage parameters and cyclic load damage parameters. Static damage parameters describe the initial damage state of the rock before it is subjected to load. In this embodiment, based on the geological data of the selected shale gas field, it is assumed that the initial damage state of the rock follows a Weibull distribution, and the static damage parameters are determined. These static damage parameters include static scale parameters. And static shape parameter m, set m=5; Cyclic load damage parameters are used to clarify the physical quantities of rock under cyclic loading and friction damage physical mechanisms. These parameters are determined based on the number of pulse fracturing cycles and pressure amplitude, and include the reference friction coefficient, direct rate effect parameters, state evolution effect parameters, characteristic slip distance, damage accumulation coefficient, phase field length scale parameters, etc. In this embodiment, m is set as follows: a=0.01, b=0.015 With λ=0.1, based on the initial fracture energy and the existing AT2 type phase field fracture model, the phase field length scale parameter is determined to be 0.5m; S5: For the rock loading model, an alternating iterative algorithm with dynamic local mesh adaptation is adopted to decouple the displacement field, friction state variables, total damage variables, and phase field variables, solve the phase field-friction damage coupled control equations, and output the simulation results of rock crack propagation.

[0029] Furthermore, the specific steps for solving the phase-field-friction damage coupling control equations are as follows: S51: Spatial discretization of the rock loading model is performed, dividing it into a uniform coarse grid. Temporal discretization of the phase field-friction damage coupled control equation set is performed, dividing it into multiple time steps. Furthermore, if the basic rock geological model is a two-dimensional planar model, the mesh elements are three-node triangular elements or four-node quadrilateral elements; if the basic rock geological model is a three-dimensional planar model, the mesh elements are four-node tetrahedral elements or eight-node hexahedral elements. The same linear Lagrangian function is used to interpolate and distribute the displacement field, friction state variables, total damage variables, and phase field variables in each mesh element to simplify coupled calculations. The implicit backward Euler method is used for time discretization to ensure solution stability.

[0030] S52: At the current time step, using the current mesh of the rock loading model, the alternating iterative solution algorithm is used to iteratively solve and update the displacement field, friction state variables, total damage variables, and phase field variables until convergence. Furthermore, the specific method for iteratively and alternately solving and updating the displacement field, friction state variables, total damage variables, and phase field variables is as follows: At the current time step Next, freeze the phase field variables solved in the previous iteration. Friction state variables Total damage variable n = 1, 2, 3, ..., N, where N is the total number of time steps, and k = 1, 2, 3, ..., K, where K is the upper limit of the number of iterations; based on the stiffness degradation coupling relationship, the reduced stiffness matrix is ​​constructed using the phase field variables and friction state variables, and substituted into the discretized mechanical equilibrium equations to solve for the displacement field of this iteration. The mechanical equilibrium equation is: ; Where σ is the stress tensor , For volume forces, the solution domain of the mechanical equilibrium equations is Ω. At the boundary of the solution domain, the following conditions are met: , , Given the displacement boundary function; The displacement field obtained by this iteration Calculate the displacement field increment, and update the equivalent sliding rate based on the displacement field increment. The friction state variables for this iteration are solved based on the velocity-state-dependent friction law. Then calculate and update the friction damage increment at the current time step. The total damage variable is updated to , This represents the total damage variable obtained at the previous time step; Based on establishing the coupling relationship between friction state variables and fracture energy, the reduced fracture energy is updated using the friction state variables and total damage variables obtained in this iteration. , ; In the formula, δ is the coupling coefficient, and in this embodiment, δ=0.3; based on the historical maximum tensile strain energy The fracture driving force is calculated based on the current stress state; the reduced fracture energy and fracture driving force are then substituted into the phase field evolution equation to solve for the phase field variables in this iteration. The phase field evolution equation is obtained through total energy functional variational analysis, and is as follows: ; in, This is the internal length scale.

[0031] S53: Determine whether the current time step meets the grid update conditions. If yes, adaptively refine and coarse the current grid and update the current grid. If not, go to S55. Furthermore, the grid update condition is that the current time step and the time step of the last grid update meet a set time step interval. In this embodiment, the set time step interval is 5 time steps. The specific method for adaptive densification and coarsening of the current grid is as follows: Calculate the phase field gradient magnitude of each cell in the current mesh. It is used to identify the location of cracks and the damage front region. The phase field gradient mode reaches its maximum value near the crack surface and decays rapidly towards the intact material region. Set the first threshold With the second threshold , The grid cells satisfy Marked as crack zone elements, the mesh elements satisfy... Marked as process area element, the mesh element satisfies Marked as a complete region unit; The elements marked as crack zones and their adjacent elements are refined. For example, a four-node tetrahedral element is divided into four sub-elements, and an eight-node hexahedral element is subdivided into an octree. For elements marked as complete zones, their relationship with adjacent elements is checked, and these elements are merged to form coarse mesh elements while satisfying geometric constraints and numerical accuracy. For elements marked as process zones, a mesh with gradually changing size is set to ensure numerical stability.

[0032] S54: Map the displacement field, friction state variables, total damage variables, and phase field variables from the previous mesh to the current mesh using linear interpolation; S55: Determine whether all time steps have been solved. If yes, output the simulation results of rock crack propagation. If not, return to S52.

[0033] Furthermore, the pulsed fracturing load conditions set in this embodiment are replaced with conventional fracturing load conditions, that is, the pressure amplitude is set to 0, and rock crack propagation simulation is performed again. The simulation results of rock crack propagation under conventional fracturing load conditions and pulsed fracturing load conditions are compared. Figure 3 As shown: Under conventional fracturing load conditions, only one main fracture extending along the direction of the maximum principal stress is formed on the rock loading model, with a length of about 35m, and the natural weak surfaces are hardly activated; while under pulse fracturing load conditions, the main fracture bifurcates multiple times during the propagation process and successfully activates 12 natural weak surfaces, forming a complex fracture network with the main fracture as the trunk and multiple secondary fractures as branches. The width of the modified area is about 25m, which is about 62.5m on the field scale according to the scaling ratio. Actual microseismic monitoring data from this shale gas field shows that during a conventional fracturing operation in the same block, the width of the fracture zone obtained by microseismic inversion was approximately 40m to 60m. In a pulse fracturing test, microseismic data showed that the width of the fracture zone increased to 60m to 80m, and the fracture complexity increased by about 30%, verifying the accuracy of the rock crack propagation simulation method used in this embodiment. The solution process outputs the damage evolution results, such as Figure 4 The damage evolution results shown in the figure are as follows: during the pressure rise phase, the damage is mainly concentrated at the tip of the main crack and at the ends of some naturally weak surfaces with favorable orientations; while during the pressure unloading phase, under the action of the friction damage physical mechanism, significant new damage still grows in the middle of some naturally weak surfaces that are obliquely related to the main crack, thus explaining the physical mechanism of new cracks induced during the unloading phase. This phenomenon cannot be captured by existing rock crack propagation simulation methods, indicating that using the friction damage physical mechanism to describe the physical process of rock crack propagation under pulse fracturing and introducing this physical mechanism into the phase field method model can more realistically simulate the rock crack network compared with existing crack propagation simulation methods. The pressure-crack volume curve output during the solution process, such as Figure 5 As shown, under the action of pulsed fracturing load, the crack volume curve shows a sawtooth upward trend, and the total crack volume growth rate is significantly higher than that under conventional fracturing load. This indicates that pulsed fracturing load can more effectively open the fracture space. This invention can accurately simulate the rock fracture network under different pulsed fracturing load conditions, which is beneficial for optimizing the pulsed fracturing parameters when designing reservoir stimulation projects. The final mesh distribution after the solution is completed shows that the mesh is highly refined near the crack path, with the smallest element size being approximately 0.125m, while a coarse mesh of 2m is maintained in the elastic region far from the crack. The total number of degrees of freedom of the rock loading model increases from the initial approximately 5,000 to approximately 18,000. Compared with a globally uniform fine mesh with all elements of 0.125m and more than 320,000 degrees of freedom, this embodiment improves the computational efficiency by approximately 18 times while ensuring high resolution in the crack region. This invention reduces the waste of computational resources while ensuring the accuracy of rock crack propagation simulation.

[0034] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A method for simulating rock crack propagation under pulsed fracturing, characterized in that, Includes the following steps: S1: Introduce phase field variables to describe the fracture state of rock materials, combine them with the displacement field of the rock, establish a total energy functional including elastic strain energy density and fracture energy, and construct a phase field method model. S2: Introducing frictional state variables, based on the rate-state dependent frictional law and the sliding law, the evolution law of frictional damage under pulse fracturing is determined, and the physical mechanism of frictional damage regarding the total damage variable is established. S3: Introduce the physical mechanism of friction damage into the phase-field method model, and establish the coupling relationship between friction state variables and elastic strain energy density and fracture energy based on the stiffness and strength properties of rock materials, respectively, and construct the phase-field-friction damage coupled control equation set. S4: Establish a basic rock geological model based on the geological conditions of the reservoir to be modified, set relevant parameters such as rock material parameters, pulse fracturing load conditions, boundary conditions, and friction damage physical mechanism, and construct a rock loading model under pulse fracturing. S5: For the rock loading model, an alternating iterative algorithm with dynamic local mesh adaptation is adopted to decouple the displacement field, friction state variables, total damage variables, and phase field variables, solve the phase field-friction damage coupled control equations, and output the simulation results of rock crack propagation.

2. The method for simulating rock crack propagation under pulsed fracturing according to claim 1, characterized in that, The total energy functional Π established in S1 is: ； In the formula, u is the displacement field. For phase field variables, Indicates that the rock material is intact. This indicates that the rock material has completely fractured. This indicates that the rock material is in the process of fracture. For elastic strain energy density, Let γ be the fracture energy, Ω be the phase field surface density function, Ω be the rock solid domain, and ε be the strain field; a degradation function is introduced into the elastic strain energy density. To reduce the elastic strain energy of the damage, ,in, It represents tensile strain energy.

3. The method for simulating rock crack propagation under pulsed fracturing according to claim 2, characterized in that, The specific method of S2 is as follows: Based on the rate-state dependent friction law, the relationship between the friction coefficient μ, the equivalent sliding speed V, and the friction state variable θ is constructed as follows: ； In the formula, a is the direct rate effect parameter, a>0, and b is the state evolution effect parameter, b>a. , , These are the reference friction coefficient, reference sliding speed, and reference friction state variable, respectively. Based on the law of sliding, the evolution law relationship of the friction state variable θ is established: ； In the formula, L is the characteristic sliding distance. When V≈0, it means that the current position is in a "healing" or "aging" state. When V>0, it means that the current position is being updated. The total damage variable D under pulse fracturing is defined as the frictional damage increment in each load cycle. The increase in frictional damage is obtained by gradually accumulating it over time, and it follows the rate-state-dependent friction law, that is, ； In the formula, λ is the damage accumulation coefficient, 0 < λ < 1.

4. The method for simulating rock crack propagation under pulsed fracturing according to claim 3, characterized in that, The specific method for establishing the coupling relationship between frictional state variables and elastic strain energy density in S3 is as follows: The degenerate function... Corrected to phase field variables Bivariate function related to friction state variable θ , ； where \(q\) is a regularization constant with a value of , \(h(\theta)\) is the friction damage factor, , \(0 < h(\theta)\leq1\), \(\alpha\) is the coupling coefficient, and \(()\) + denotes taking the positive value; establish the stiffness degradation coupling relation: , \(C\) is the elastic stiffness tensor, is the initial stiffness tensor in the intact state of the rock material; establish the coupling relation between the friction state variable and the elastic strain energy density, that is .

5. The method for simulating rock crack propagation under pulsed fracturing according to claim 4, characterized in that, The specific method for establishing the coupling relationship between frictional state variables and fracture energy in S3 is as follows: The fracture energy... Modified to a function related to the friction state variable θ , ； In the formula, Let β be the initial fracture energy, and β be the wear coefficient, where 0 ≤ β < 1.

6. The method for simulating rock crack propagation under pulsed fracturing according to claim 5, characterized in that, The specific steps of S5 are as follows: S51: Spatial discretization of the rock loading model is performed, dividing it into a uniform coarse grid. Temporal discretization of the phase field-friction damage coupled control equation set is performed, dividing it into multiple time steps. S52: At the current time step, using the current mesh of the rock loading model, the alternating iterative solution algorithm is used to iteratively solve and update the displacement field, friction state variables, total damage variables, and phase field variables until convergence. S53: Determine whether the current time step meets the grid update conditions. If yes, adaptively refine and coarse the current grid and update the current grid. If not, go to S55. S54: Map the displacement field, friction state variables, total damage variables, and phase field variables from the previous mesh to the current mesh using linear interpolation; S55: Determine whether all time steps have been solved. If yes, output the simulation results of rock crack propagation. If not, return to S52.

7. The method for simulating rock crack propagation under pulsed fracturing according to claim 6, characterized in that, The specific method for iteratively solving and updating the displacement field, friction state variables, total damage variables, and phase field variables in S52 is as follows: Freeze the phase field variables, friction state variables, and total damage variables obtained in the previous iteration, construct the reduced stiffness matrix using the stiffness degradation coupling relationship, substitute it into the discretized mechanical equilibrium equations, and solve for the displacement field in this iteration. The mechanical equilibrium equations are as follows: Where σ is the stress tensor, , For volume forces, satisfying the following conditions at the boundary of the rock solid domain: , , Given the displacement boundary function; The equivalent sliding rate is updated using the displacement field obtained from this iteration. The friction state variables for this iteration are solved based on the velocity-state-dependent friction law. The total damage variable at the current time step is then calculated and updated. Based on establishing the coupling relationship between friction state variables and fracture energy, the reduced fracture energy is updated using the friction state variables and total damage variables obtained in this iteration. The energy is then substituted into the discretized phase field evolution equation to solve for the phase field variables in this iteration. The phase field evolution equation is obtained through total energy functional variation.

8. The method for simulating rock crack propagation under pulsed fracturing according to claim 7, characterized in that, The specific method for adaptive refinement and coarsening of the current mesh in S53 is as follows: Calculate the phase field gradient magnitude of each cell in the current mesh. ; Set the first threshold With the second threshold , The grid cells satisfy Marked as crack zone elements, the mesh elements satisfy... Marked as process area element, the mesh element satisfies Marked as a complete region unit; The elements marked as crack zones and their adjacent elements are refined; the elements marked as intact zones are examined for their relationship with adjacent elements, and these elements are merged to form coarse mesh elements while satisfying geometric constraints and numerical accuracy; the elements marked as process zones are set to meshes with gradually changing sizes.

9. The method for simulating rock crack propagation under pulsed fracturing according to claim 1, characterized in that, The specific method for establishing the basic rock geological model in S4 is as follows: Based on the geological conditions of the reservoir to be modified, the three-dimensional spatial structure of the reservoir to be modified is extracted, and a corresponding basic rock geological model is established. A main fracture of a certain length is preset in the center of the basic rock geological model to simulate the initial fracture formed by the initiation of the perforation cluster. Several natural weak surfaces are randomly distributed around the main fracture, with a length of 2.5%-15% of the main fracture. The dip angle of the randomly distributed natural weak surfaces is reduced by a reduction factor, which is in the range of 0 to 1.

10. The method for simulating rock crack propagation under pulsed fracturing according to claim 9, characterized in that, The specific method for setting the pulse fracturing load conditions and boundary conditions in S4 is as follows: a pulsed fluid pressure load p(t) is applied to the main fracture surface. ,in, Based on pressure, denoted as pressure amplitude and f as pulse frequency; a constant minimum horizontal principal stress is applied around the rock geological foundation model, with vertical displacement constrained at the bottom and horizontal displacement constrained on both sides.

Images