Hydraulic fracture evolution simulation method and device and electronic equipment

Abstract

The invention provides a hydraulic fracture evolution simulation method and device and electronic equipment.The hydraulic fracture evolution simulation method comprises the steps that a damaged rock mass is divided into a plurality of different calculation areas, and a seepage control equation of each calculation area is determined; in combination with seepage anisotropy, determining a seepage tensor of porous medium seepage-stress coupling in each calculation area; associating the seepage control equations of all the calculation areas with phase fields to obtain a target seepage control equation; and determining the displacement, fluid pressure and phase field of the damaged rock mass based on the target seepage control equation. Through the deduced seepage-stress coupling phase field model, the strain energy is decomposed through spectral decomposition, it is guaranteed that crack propagation occurs under the energy generated by tensile strain, and the problem that an existing hydraulic crack research method is limited is solved.

Description

Technical Field

[0001] The present invention relates to the technical field of geological disaster prevention and control, and particularly to a method, device and electronic equipment for simulating the evolution of hydraulic fractures. Background Art

[0002] In geotechnical engineering, water is a key factor affecting engineering safety, and fractures are important seepage channels for water, which play a promoting role in the initiation and propagation of fractures. For example, in the case of a rock slope on the reservoir bank under the action of a dynamic water level, due to the seepage effect of the rock mass, a water level difference will be formed between the interior of the slope and the reservoir water surface, thereby generating a dynamic water pressure, which promotes the sliding of the slope rock mass and causes crack propagation, resulting in landslides. Since the operation of the Three Gorges Reservoir area, there have been more than 4,000 landslides, such as the Bazimen Landslide and Liangshuijing Landslide in the Three Gorges Reservoir area. In addition, in oil and gas exploitation, hydraulic fracturing technology is also used to promote the formation of a fracture network in the formation, increase the permeability of the rock formation, and thus increase oil and gas production. Therefore, in order to further study the influence of water seepage on the safety of slope engineering and the effect of hydraulic fracturing technology on oil and gas exploitation, the present invention will focus on the influence of seepage-stress coupling on crack propagation.

[0003] At present, indoor experimental research cannot meet the research needs, and numerical simulation technology is currently the mainstream means to solve geotechnical engineering problems. The important methods for numerical simulation of seepage-stress coupling problems are the discrete method and the continuous method. Among the continuous methods, the phase field method is one of the most popular research methods in recent years. Therefore, this patent adopts the phase field method to study the seepage-stress coupling model, aiming to consider the influence of water seepage on crack initiation and propagation.

[0004] Regarding the problem of limitations in existing research methods for hydraulic fractures, no good solution has been proposed yet. Summary of the Invention

[0005] The present invention provides a method, device and electronic equipment for simulating the evolution of hydraulic fractures to solve the defect of limitations in existing research methods for hydraulic fractures.

[0006] In a first aspect, the present invention provides a method for simulating the evolution of hydraulic fractures, including: Dividing the damaged rock mass into several different calculation regions and determining the seepage control equation for each of the calculation regions; Combining seepage anisotropy to determine the permeability tensor of porous media seepage-stress coupling for each of the calculation regions; Associating the seepage control equations of all the calculation regions with the phase field to obtain a target seepage control equation; Based on the target seepage control equation, determining the displacement, fluid pressure and phase field of the damaged rock mass.

[0007] A method for simulating the evolution of hydraulic fractures provided by the present invention divides the damaged rock mass into several different calculation regions, including: Dividing the damaged rock mass into a matrix region, a crack region, and a transition region; Wherein, the transition region is a calculation region located between the matrix region and the crack region.

[0008] A method for simulating the evolution of hydraulic fractures provided by the present invention determines the seepage control equation of the matrix region, including: Based on Darcy's law and Biot's theory, obtaining the storage coefficient and fluid velocity of the matrix region; Generating the seepage control equation of the matrix region according to the storage coefficient and fluid velocity of the matrix region.

[0009] A method for simulating the evolution of hydraulic fractures provided by the present invention determines the permeability tensor of the seepage-stress coupling of the porous medium in the matrix region, including: Obtaining the initial permeability coefficient of the matrix region; Based on the product of the initial permeability coefficient of the matrix region and the right Cauchy-Green strain tensor, determining the permeability tensor of the seepage-stress coupling of the porous medium in the matrix region.

[0010] A method for simulating the evolution of hydraulic fractures provided by the present invention determines the permeability tensor of the seepage-stress coupling of the porous medium in the crack region, including: Determining the first permeability tensor in the open state and the second permeability tensor in the closed state of the crack region; Combining the first permeability tensor and the second permeability tensor to obtain the permeability tensor of the seepage-stress coupling of the porous medium in the crack region.

[0011] A method for simulating the evolution of hydraulic fractures provided by the present invention, before correlating the seepage control equations of all the calculation regions with the phase field, includes: Interpolating the permeability tensors of the crack region and the matrix region to obtain the permeability tensor of the entire calculation region of the damaged rock mass.

[0012] A method for simulating the evolution of hydraulic fractures provided by the present invention correlates the seepage control equations of all the calculation regions with the phase field to obtain a target seepage control equation, including: Determining the total internal energy of the system where the damaged rock mass is located; Combining the total internal energy of the system, and according to the variational principle and the principle of minimum energy, obtaining the phase field control equation of the seepage-stress coupling; Introducing a historical energy field to determine the weak form of each physical field in the seepage-stress coupling model; the seepage-stress coupling model includes a stress field, a phase field, and a seepage field; Perform finite element discretization on the weak form of each of the physical fields to obtain the residuals of the corresponding physical fields; Generate the target seepage control equation by combining the residuals of the physical fields.

[0013] According to a hydraulic fracture evolution simulation method provided by the present invention, based on the target seepage control equation, determine the displacement, fluid pressure, and phase field of the damaged rock mass, including: Solve for the displacement and fluid pressure of each physical field in the seepage-stress coupling model by means of an alternating iteration method, and determine the displacement vector of the nodes by the fluid pressure of the fixed nodes until a preset convergence criterion is satisfied; Determine the field variables of the historical energy field by the displacement and fluid pressure of the nodes; Update the field variables of the historical energy field, and obtain the phase field of the damaged rock mass by the Newton-Raphson method.

[0014] In a second aspect, the present invention further provides a hydraulic fracture evolution simulation device, including: A division module, configured to divide the damaged rock mass into a plurality of different calculation regions, and determine the seepage control equation of each of the calculation regions; A processing module, configured to determine the permeability tensor of the porous medium seepage-stress coupling of each of the calculation regions by combining seepage anisotropy; An association module, configured to associate the seepage control equations of all the calculation regions with the phase field to obtain the target seepage control equation; A determination module, configured to determine the displacement, fluid pressure, and phase field of the damaged rock mass based on the target seepage control equation.

[0015] In a third aspect, the present invention further provides an electronic device, including a memory, a processor, and a computer program stored on the memory and executable on the processor, where when the processor executes the program, the hydraulic fracture evolution simulation method described in the first aspect above is implemented.

[0016] In a fourth aspect, the present invention further provides a non-transitory computer-readable storage medium, on which a computer program is stored, and when the computer program is executed by a processor, the hydraulic fracture evolution simulation method described in the first aspect above is implemented.

[0017] In a fifth aspect, the present invention further provides a computer program product, including a computer program, and when the computer program is executed by a processor, the hydraulic fracture evolution simulation method described in the first aspect above is implemented.

[0018] Compared with the prior art, the present invention has the following beneficial effects: The hydraulic fracture evolution simulation method provided by the present invention deduces the seepage-stress coupling phase field model, decomposes the strain energy by spectral decomposition, ensures that crack propagation occurs under the energy generated by tensile strain, considers the influence of water work on crack propagation, and the permeability coefficient changes with the crack aperture and rock mass deformation, solving the problem of limitations in existing research methods for hydraulic fractures. Description of the Drawings

[0019] In order to more clearly illustrate the technical solutions in the present invention or the prior art, the following will briefly introduce the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings in the following description are some embodiments of the present invention. For those of ordinary skill in the art, without creative efforts, other drawings can also be obtained based on these drawings.

[0020] Figure 1 is the flowchart of the hydraulic fracture evolution simulation method provided by the present invention; Figure 2 is a schematic diagram of the boundary and geometric shape of an impermeable rectangular plate in an embodiment of the present invention; Figure 3 is a schematic diagram of the comparison results of the theoretical and numerical solutions with different time lengths in an embodiment of the present invention; Figure 4 is a schematic diagram of the influence of different roughness coefficients on the water pressure distribution when T is taken as 0.1 in an embodiment of the present invention; Figure 5 is a schematic diagram of the geometric dimensions and boundary conditions of the Sneddon analytical solution model in an embodiment of the present invention; Figure 6 is a schematic diagram of the influence of element size on the numerical simulation results in the Sneddon analytical solution in an embodiment of the present invention; Figure 7 is a schematic diagram of the relationship curve between the fluid pressure and fluid volume at the injection point in an embodiment of the present invention; Figure 8 is a schematic diagram of the phase field, water pressure, and Y-direction displacement distributions of fluid injection-induced crack propagation in an embodiment of the present invention; Figure 9 is a schematic diagram of the change process of the crack aperture in an embodiment of the present invention; Figure 10 is a schematic diagram of the multi-layer rock mass hydraulic fracture propagation model in an embodiment of the present invention; Figure 11 is the phase field cloud map and horizontal displacement cloud map when the rock mass elastic modulus ratio is 2 at t = 200s in an embodiment of the present invention; Figure 12 is the water pressure curve graph when the rock mass elastic modulus ratio E1 / E2 = 2 in an embodiment of the present invention; Figure 13 Schematic diagram of the horizontal displacement along the Parh-1 path on the crack surface when the elastic modulus ratio E1 / E2 = 2; Figure 14 It is the phase field cloud map and the horizontal position cloud map with a rock stratum elastic modulus ratio of 4 at t = 200 s in the embodiment of the present invention; Figure 15 It is the water pressure curve graph when the rock stratum elastic modulus ratio E1 / E2 = 4 in the embodiment of the present invention; Figure 16 Schematic diagram of the horizontal displacement along the Parh-1 path on the crack surface when the elastic modulus ratio E1 / E2 = 2 in the embodiment of the present invention; Figure 17 It is the structural block diagram of the hydraulic fracture evolution simulation device provided by the present invention; Figure 18 It is the structural schematic diagram of the electronic device provided by the present invention. Detailed implementation manners

[0021] To make the objectives, technical solutions, and advantages of the present invention clearer, the technical solutions in the present invention will be clearly and completely described below with reference to the accompanying drawings in the present invention. Apparently, the described embodiments are some, but not all, of the embodiments of the present invention. All other embodiments obtained by those of ordinary skill in the art based on the embodiments in the present invention without making creative efforts shall fall within the protection scope of the present invention.

[0022] The present invention provides a method for simulating hydraulic fracture evolution, Figure 1 It is the flow chart of the method for simulating hydraulic fracture evolution provided by the present invention. As shown in the figure, the method includes the following steps: Step S101: Divide the damaged rock mass into several different calculation regions and determine the seepage control equation for each calculation region; Step S102: Combine seepage anisotropy to determine the permeability tensor of the porous medium seepage-stress coupling for each calculation region; Step S103: Correlate the seepage control equations of all calculation regions with the phase field to obtain the target seepage control equation; Step S104: Based on the target seepage control equation, determine the displacement, fluid pressure, and phase field of the damaged rock mass.

[0023] In this method, first, the damaged rock mass is divided into calculation regions according to its structural components, and the seepage control equation for each calculation region is determined. Then, considering seepage anisotropy, the permeability tensor of the seepage-stress coupling in porous media is derived. Next, the seepage control equations of all calculation regions are correlated with the phase field to generate a unified target seepage control equation. Finally, the target seepage control equation is derived and the algorithm is implemented to obtain the displacement, fluid pressure, and phase field of the damaged rock mass. In this process, a seepage-stress coupling phase field model is derived, the spectral decomposition is used to decompose the strain energy, ensuring that crack propagation occurs under the energy generated by tensile strain, considering the influence of water work on crack propagation, and the permeability coefficient changes with the crack aperture and rock mass deformation, solving the problem of limitations in existing research methods for hydraulic fractures.

[0024] In some of these embodiments, in step S101, the damaged rock mass is divided into several different calculation regions, including: dividing the damaged rock mass into a matrix region , a crack region , and a transition region . Among them, the transition region is the calculation region located between the matrix region and the crack region .

[0025] Based on this embodiment, the seepage control equation of the matrix region is determined, including: based on Darcy's law and Biot's theory, obtaining the storage coefficient and fluid flow velocity of the matrix region; according to the storage coefficient and fluid flow velocity of the matrix region, generating the seepage control equation of the matrix region.

[0026] Furthermore, the permeability tensor of the seepage-stress coupling in the porous media of the matrix region is determined, including: obtaining the initial permeability coefficient of the matrix region; based on the product of the initial permeability coefficient of the matrix region and the right Cauchy-Green strain tensor, determining the permeability tensor of the seepage-stress coupling in the porous media of the matrix region.

[0027] Exemplarily, for the matrix region , regarding it as a porous medium, using Darcy's law and Biot's theory to describe the flow of fluid in the porous medium, according to the mass conservation equation, the seepage control equation can be expressed as:

[0028] Among them, is the fluid density, represents the storage coefficient of the matrix region , represents the Biot coefficient of the matrix region , , Represents the strain energy density caused by fluid pressure, Represents the matrix region of the volume strain, Represents the fluid source or sink flux, p Represents the fluid pressure, t Represents time, Represents the fluid flow velocity, storage coefficient Can be calculated according to the following formula:

[0029] Wherein, Represents the porosity of the rock, Represents the compressibility coefficient of the fluid, Represents the bulk modulus of the fluid.

[0030] For the fluid flow velocity, it can be calculated by Darcy's law, and the formula is as follows:

[0031] Wherein, Represents the viscosity of the fluid, Represents the permeability tensor at this point, Represents the acceleration of gravity. Considering the influence of rock mass deformation on the permeability coefficient, the matrix region of the permeability tensor can be expressed as:

[0032] Wherein, Represents the matrix region Initial permeability coefficient, Represents the deformation gradient matrix Determinant value ( ), Represents the right Cauchy-Green strain tensor.

[0033] Determine the permeability tensor of the seepage-stress coupling in the crack zone porous medium, including determining the first permeability tensor in the open state and the second permeability tensor in the closed state of the crack zone; combining the first permeability tensor and the second permeability tensor to obtain the permeability tensor of the seepage-stress coupling in the crack zone porous medium.

[0034] Exemplarily, for the crack zone , the seepage control equation should ignore the influence of the volume expansion of the rock mass on seepage, and the seepage control equation can be expressed as:

[0035] Wherein, Is the fluid density, Represents the storage coefficient of the crack zone, Denotes the fluid flow velocity, Denotes the fluid source or sink flux. For the crack region The fluid flow velocity can be calculated by the following formula:

[0036] Where, Denotes the permeability tensor in the crack region And Denotes the viscosity of the fluid.

[0037] For the crack region , the effects of deformation and damage on the permeability coefficient need to be considered simultaneously. Assuming that the fluid flow velocity in the crack region Complies with Darcy's law, the fluid in the crack region can be assumed to be Darcy-Poiseuille flow. If the crack is in an open state, the permeability tensor can be derived from the cubic law, and the relationship between the first permeability tensor and the crack opening can be expressed as:

[0038] Where, Denotes the crack opening. Under tensile stress, the crack opening can be expressed as:

[0039] Where, Denotes the characteristic length of the element, Denotes taking the trace operation on the strain tensor, The calculation method of Is defined as d And only when in a tensile state will the crack have an opening; Denotes the phase field, Denotes a threshold, and , which means that the crack opening state is considered only when the phase field is greater than

[0040] Where, Denotes the unit normal vector of the crack surface, C denotes the right Cauchy-Green strain tensor. For simplicity of calculation, the unit vector of the gradient of the phase field is approximated as the normal vector of the crack surface, so as to describe the deformation on both sides of the crack.

[0041] When the crack is in a closed state, the residual permeability coefficient of crack closure needs to be considered, so the second permeability tensor is expressed as:

[0042] Where, ​Denotes the permeability coefficient when the crack is closed and the normal stress is equal to 0. is a bracket operator, when then When then ; Denotes the tensile length along the normal direction of the crack surface. When the crack is in an open state, and when the crack is in a closed state. Calculating this parameter can determine the opening and closing state of the crack, combined with the normal direction of the crack surface:

[0043] Among them, Denotes the unit vector perpendicular to the crack surface in the initial configuration. Denotes the unit vector perpendicular to the crack surface in the current configuration, and Can be expressed as:

[0044] Among them, F represents the deformation gradient matrix, and d represents the phase field. Combining the above formulas, we can get:

[0045] Among them, C represents the right Cauchy-Green strain tensor.

[0046] In actual engineering rock masses, the roughness of the structural plane will affect the seepage performance. It is necessary to consider the influence of the crack surface roughness on seepage. Combining the formulas, the permeability tensor in the crack area can be expressed as:

[0047] Among them, Is a quantity related to roughness, and the relationship with roughness can be expressed as , if the value of JRC is not clear, the default value is 1.

[0048] On this basis, before correlating the seepage control equations of all calculation regions with the phase field, including: interpolating the permeability tensors of the crack area and the matrix area to obtain the permeability tensor of the entire calculation region of the damaged rock mass.

[0049] Exemplarily, for the crack area and the matrix area Interpolating the permeability tensors can obtain the expression of the permeability tensor of the entire calculation region:

[0050] In order to uniformly express the seepage control equations of the three regions during the numerical simulation process, it is necessary to perform interpolation processing on the three divided calculation domains, so an interpolation function is introduced and Differences are taken for the seepage parameters of the three computational domains, and the formula of the interpolation function is as follows:

[0051]

[0052] When it indicates that the computational domain is the matrix region and when it indicates that the computational domain is the transition region and when it indicates that the computational domain is the crack region . The seepage equilibrium equation introducing the interpolation function can be uniformly expressed as:

[0053] In the above formula, , , , , , , the Biot coefficient of the crack region is equal to 1, and the porosity is equal to 0.

[0054] In some of the embodiments, in step S103, the seepage control equations of all computational domains are associated with the phase field to obtain the target seepage control equation, including: determining the total internal energy of the system where the damaged rock mass is located; combining the total internal energy of the system, according to the variational principle and the principle of minimum energy, obtaining the phase field control equation of seepage-stress coupling; introducing the historical energy field and determining the weak forms of each physical field in the seepage-stress coupling model; the seepage-stress coupling model includes a stress field, a phase field, and a seepage field; performing finite element discretization on the weak forms of each physical field to obtain the residuals of the corresponding physical fields; generating the target seepage control equation by combining the residuals of the physical fields.

[0055] Exemplarily, the total internal energy of the system can be expressed as:

[0056] Among them, represents the total internal energy of the system, represents the elastic energy of the system, represents the surface energy of the crack, represents the work done by the water pressure on the rock mass, represents the external work. The elastic energy is obtained by spectral decomposition, and the specific formula is as follows:

[0057] Among them, represents the stretching energy, Represents the energy of the compressed part. The surface energy of the crack can be obtained by multiplying the fracture energy release rate and the surface area. The specific formula is as follows:

[0058] in, represents the critical energy release rate, which is a material property related to the material. Represents the crack surface. Convert the surface integral into volume integral and substitute it into the density function of the crack. The surface energy of the crack can be calculated by the following formula:

[0059] Among them, t represents the external force, u represents the displacement field, represents the boundary of node N, represents the surface integral area of ​​the crack surface, represents the fluid pressure in the crack zone, represents the unit vector perpendicular to the crack surface in the current configuration, Represents the entire computational domain.

[0060] Considering the fluid pressure on the crack surface Continuity, so the fluid pressure in the crack area satisfy , in the matrix area satisfy , on the crack surface, it can be considered that the matrix area and crack zone The fluid pressure is equal, , through the divergence theorem and the degradation function, the formula can be transformed into:

[0061] The combined formulas, the total internal energy expression of the system can be expressed as:

[0062] According to the variational principle and the energy minimization principle, the phase field governing equation of seepage-stress coupling is derived as follows:

[0063] in, represents the gradient operator, represents the stress tensor, α represents the effective stress coefficient, g (·) represents the degradation function, I represents the unit tensor, b represents the physical force, g c represents the critical energy release rate, l 0 represents the characteristic width of the phase field model, Represents the first derivative of the degradation function.

[0064] To ensure the irreversibility of the phase field, consistent with the traditional phase field, a historical energy field is introduced to ensure that the phase field remains irreversible even when the boundary conditions change. The expression of the historical energy field of the seepage-stress coupling phase field damage model is:

[0065] where, represents the historical field variable, S represents time, and [0, T] represents the time period. The historical field variable is introduced, and according to the divergence theorem, the formula can be reformulated as:

[0066] The seepage-stress coupling model contains three physical fields, namely the stress field, the phase field, and the seepage field. The regularized crack evolution problem can be regarded as a minimization problem of the energy functional. According to the physical field control equations, the weak form of the finite element is derived, and corresponding algorithms are needed to solve the three physical fields. The weak forms corresponding to the three physical fields are obtained by applying the principle of virtual work. The weak forms of the displacement field, the phase field, and the seepage field are respectively:

[0067]

[0068]

[0069] where, represents the increment, represents the known external force, k represents the specific spatial direction of the displacement field, and K represents the permeability tensor. And, the differential term of the formula can be expressed by the difference method:

[0070] where, represents the time step size, represents the iteration step in the Newton iteration, represents the time step. According to the derived weak form of the finite element, finite element discretization can be performed to obtain the residuals of the three physical fields:

[0071]

[0072]

[0073] where, represents the stress field residual, represents the displacement field residual caused by the external force, represents the displacement field gradient matrix, represents the shape function matrix of the displacement field, represents the phase field residual, represents the shape function matrix of the phase field, represents the phase field gradient matrix, represents the seepage field residual, represents the fluid pressure field residual caused by external forces, represents the gradient matrix of the seepage field, represents the shape function matrix of the seepage field.

[0074] The gradient matrix of the seepage field and the shape function matrix are consistent with those of the phase field. Its gradient matrix and the shape function matrix are expressed as:

[0075]

[0076] and are expressed as:

[0077]

[0078] wherein, represents the shape function matrix of the displacement field, represents the shape function matrix of the seepage field.

[0079] Based on the above embodiments, in step S104, based on the target seepage control equation, the displacements, fluid pressures, and phase fields of the damaged rock mass are determined, including: solving for the displacements and fluid pressures of each physical field in the seepage-stress coupling model through an alternating iteration method, and determining the displacement vector of the nodes through the fluid pressures of the fixed nodes until the preset convergence criterion is met; determining the field variables of the historical energy field through the displacements and fluid pressures of the nodes; updating the field variables of the historical energy field, and obtaining the phase field of the damaged rock mass through the Newton-Raphson method.

[0080] Exemplarily, the seepage-stress coupling model is iterated using the Newton-Raphson method and the alternating algorithm is used to solve for the displacements, phase fields, and fluid pressures. The stiffness matrices of the stress field, phase field, and seepage field are respectively:

[0081]

[0082]

[0083] Among them, represents the stress field stiffness matrix, where i and j represent the matrix sequence numbers, represents the phase field stiffness matrix, represents the seepage field stiffness matrix, and S represents the storage coefficient of the fluid.

[0084] The displacement, fluid pressure, and phase field are solved sequentially using the staggered iteration algorithm. In an increment step ( ), the algorithm includes three steps: 1. In the i-th iteration, the decoupled displacement and fluid pressure are solved using the staggered iteration scheme. First, the fluid pressure at the nodes is solved by fixing the displacement, and then the displacement vector at the nodes is solved by fixing the fluid pressure at the nodes. These two sub-steps are repeated until the convergence criteria and are satisfied and then the loop is exited; among them, represents the volumetric strain; 2. Calculate the historical field variables using the displacement and fluid pressure of the nodes obtained from step 1; 3. Update the historical field variables , and solve the phase field using the Newton-Raphson method.

[0085] To verify the effectiveness of this method, the following experiments were conducted: 1. An impermeable rectangular plate was used for the seepage test. The seepage parameters of the impermeable rectangular plate model are shown in Table 1: Table 1 Setting table of seepage parameters of impermeable rectangular plate model

[0086] On the left is a rectangular plate with a constant water pressure boundary condition. A fluid pressure P 0 is applied to the left boundary of the model. The displacements around are fixed. The material of the plate is impermeable. Its geometric shape and boundary conditions are as Figure 2 shown, Figure 2 which is a schematic diagram of the boundary and geometric shape of the impermeable rectangular plate in the embodiment of the present invention. The length of the rectangular plate is 1 m, the height is 0.4 m, and there is a prefabricated crack running through from left to right in the middle of the rectangular plate. Its width is set to 4.0×10 -5 m. A suddenly applied fluid pressure P 0 is applied to the left side of the model, and the water pressure distribution of the prefabricated crack from left to right is calculated. The analytical solution of the fluid pressure along the prefabricated crack is:

[0087] Among them, P is the water pressure distribution of the prefabricated crack along the positive X-axis direction, , , where is the time when the left end of the model is subjected to fluid pressure, is the bulk modulus of the fluid, is the viscosity of the fluid.

[0088] To verify the correctness of the model, the calculation parameters are set as shown in Table 1. The energy release rate of the rectangular plate is set to 100 N / m. To study the influence of different time magnitudes on the water pressure distribution, values are taken as 0.01, 0.05, 0.1, and 0.2 respectively, and the results of the theoretical solution and the numerical solution are compared and analyzed. Figure 3 is a schematic diagram of the comparison results of the theoretical and numerical solutions for different time lengths in the embodiment of the present invention. As shown in Figure 3 , for the distribution curves of water pressure at different times, the results of the numerical solution and the theoretical solution are basically the same, indicating the accuracy of the seepage-stress coupling model; and as time increases, the water pressure at the right end of the prefabricated crack shows a gradually increasing trend, indicating the influence of time on the water pressure distribution in unsteady seepage. The longer the time, the more stable the water pressure. Furthermore, introducing the influence of JRC on the seepage parameters, to verify the influence of JRC on seepage, is fixed at 0.1, and the JRC parameters are taken as 1, 2, 4, and 6 respectively. The results of the numerical simulation are as shown in Figure 4 , Figure 4 is a schematic diagram of the influence of different roughness coefficients (Joint Roughness Coefficient, JRC) on the water pressure distribution when T is taken as 0.1 in the embodiment of the present invention. The calculation results show that as JRC increases, the fluid pressure along the prefabricated crack gradually decreases, indicating that the rougher the crack, the stronger the hindering effect of the crack surface on water penetration.

[0089] 2. Conduct an analytical experiment using a pressurized crack The geometric dimensions and boundary conditions of the model are as shown in Figure 5 , Figure 5 is a schematic diagram of the geometric dimensions and boundary conditions of the Sneddon analytical solution model in the embodiment of the present invention. The model is a square plate with a width of 4 m, and a prefabricated crack with a length of 0.2 m is prefabricated in the center. A constant fluid pressure acts inside the prefabricated crack, and the value of the constant fluid pressure is P = 1.0×10 4 Pa, and displacements will occur along the Y direction on the upper and lower surfaces of the crack. The deformations of the upper and lower surfaces of the prefabricated crack are symmetric. The analytical solution of the vertical displacement perpendicular to the X direction on the upper surface of the crack is:

[0090] Where, is the length of the prefabricated crack, is the X-axis coordinate of an arbitrarily selected node on the upper surface of the prefabricated crack (the coordinate origin is the center point of the model).

[0091] To verify the correctness of the seepage-stress coupling model, the numerical simulation results are compared and analyzed with the Sneddon analytical solution of the pressurized crack. The prefabricated crack of the model is introduced through the initial historical field variables, and the phase field characteristic length is selected as = 0.02 m. To study the influence of the mesh size on the numerical results, three groups of numerical results with different element sizes are set for comparison. The model is discretized using quadrilateral elements, and the element size is set to 0.004 m, 0.006 m, and 0.01 m respectively. The elastic modulus of the model is 1 GPa, and the Poisson's ratio = 0.2, the energy release rate = 10 N / m, = 0.1 s. The seepage parameters are the same as those in Example 1, as shown in Table 1. Figure 6 is a schematic diagram showing the influence of the element size in the Sneddon analytical solution on the numerical simulation results in the embodiment of the present invention. As Figure 6 shown, first, the numerical simulation results are close to the analytical solution results, and the displacement trends are basically the same; second, as the element size increases, the displacement at the center point position gradually increases, and the displacements on both sides are less different. The element size that best matches the analytical solution is = 0.004 m, and the maximum deviation is = 0.01 m, indicating that when performing numerical simulation, the phase field method is sensitive to the mesh size, and certain requirements for the mesh size are needed when dividing the mesh to ensure the calculation accuracy.

[0092] In addition, to further verify the promoting effect of water on crack propagation in the seepage-stress coupling model, the model is discretized using quadrilateral elements, and the element size is = 0.01 m, the phase field characteristic lengths are 0.02 m, 0.03 m, and 0.04 m respectively, and the time step is taken as = 0.1 s, and the water injection rate is 1×10-4 m³ / s. Figure 7 is a schematic diagram of the relationship curve between the fluid pressure and the fluid volume at the injection point in the embodiment of the present invention. From Figure 7 it can be seen that the fluid pressure is less affected by the characteristic length, and the peak fluid pressure slightly increases with the increase of the characteristic length, and the trends are basically the same. Figure 8 is a schematic diagram showing the phase field, water pressure, and Y-direction displacement distributions of fluid injection-induced crack propagation in the embodiment of the present invention. Specifically,[[]] Figure 8 is the characteristic length When [[ID=]] = 0.02 m, the phase field contour maps, fluid pressure contour maps, and displacement contour maps at 20 s, 120 s, and 160 s are shown. It can be seen that the phase field is consistent with the fluid pressure distribution. The larger the phase field value, the more severe the rock mass damage, the stronger the diversion ability, and the greater the fluid pressure.

[0093] To study the relationship between crack evolution and crack opening, the characteristic length was extracted. The crack opening variation curve when [[ID=]] = 0.02 m. Figure 9 This is a schematic diagram of the change process of crack opening in the embodiment of the present invention. As Figure 9 shown, before point A, the crack did not expand and was in the energy accumulation stage. The crack opening changed linearly and increased with the increase of fluid pressure. After point A, the crack started to initiate and expand, the fluid pressure reached the maximum value, and the crack opening increased the fastest. After point B, the fluid pressure dropped to the lowest point, the crack propagation speed slowed down, so the growth rate of the crack opening also slowed down, and the curve slope decreased. This shows that the crack opening evolution of the seepage-stress coupling model with opening independence is basically consistent with the crack evolution law.

[0094] 3. Conduct a multi-layer hydraulic fracture propagation test. The rock mass strength parameters are shown in Table 2: Table 2 Rock mass strength parameter table

[0095] In a multi-layer geological body, when the hydraulic fracture reaches the interface between hard and soft layers, penetration or deflection will occur. To further study the influence of rock mass strength on crack propagation, a numerical model test was designed. The geometric dimensions and boundary conditions of the model are as Figure 10 shown. Figure 10 This is a schematic diagram of the multi-layer hydraulic fracture propagation model in the embodiment of the present invention. The included angle between two rock layers is selected as [[ID=]] = 0°, 15°, and 30°. Among them, layer ① is the hard rock layer and layer ② is the soft rock layer. To study the influence of the hardness of the rock mass on the crack propagation direction, two sets of parameters are designed for each rock layer included angle. The material parameters of the rock layers are shown in Table 2. The elastic modulus ratio E1 / E2 of the hard rock and the soft rock is given as 2 and 4 respectively, and the energy release rate also changes with the change of the elastic modulus. The characteristic length of the phase field is [[ID=]] = 0.04 m, the water injection rate at the injection point is [[ID=]] = 0.005 m 3 / s. The quadrilateral elements are used for discretization, the minimum grid size is [[ID=]] = 0.015 m, the time interval is [[ID=]] = 0.1 s, the Poisson's ratio is [[ID=]] = 0.3, the permeability coefficient is [[ID=]] = 1×10-14 m 2 and is [[ID=]] = 1×10-6 m2 , and the remaining seepage parameters remain the same as those in Table 1.

[0096] When the ratio of the elastic modulus of the hard rock layer to the soft rock layer E1 / E2 is 2, the phase field contour map and the horizontal displacement contour map at t = 200 s are as Figure 11 shown. Figure 11 It is the phase field contour map and the horizontal displacement contour map of the rock layer with an elastic modulus ratio of 2 at t = 200 s in the embodiment of the present invention. When the rock layer included angle is 0° and 15°, the hydraulic fracture vertically penetrates the hard rock layer, as Figure 11 (a) and Figure 11 (b) shown. The main reason is that the hardness of the hard rock layer is insufficient to prevent the crack from penetrating the rock layer interface. As can be seen from Figure 11 (c), when the rock layer included angle is 30°, the crack propagates upward along the rock layer interface. Figure 12 It is the water pressure curve graph of the rock layer with an elastic modulus ratio E1 / E2 = 2 in the embodiment of the present invention. As Figure 12 shown, the crack starts to propagate at point A and reaches the rock layer boundary at point B. When the rock layer included angle is too large, such as when the rock layer included angle reaches 30°, more energy accumulation is required for the crack to deflect. Therefore, the water pressure at which the crack deflects when the rock layer included angle is 30° is greater than the water pressure when the rock layer included angles are 0° and 15°; Figure 13 is the schematic diagram of the horizontal displacement along the crack surface in the Parh-1 path when the elastic modulus ratio E1 / E2 = 2. As can be seen from Figure 13 it, because the water pressure when the rock layer included angle is 30° is greater than the water pressure when the rock layer included angles are 0° and 15°, the horizontal displacement of the rock mass crack surface when the rock layer included angle is 30° is greater than the displacement of the other included angles. In the latter half, because the crack does not penetrate the rock layer, the displacement of the crack surface is less than the horizontal displacement when the rock layer included angles are 0° and 15°. Because the water pressures when the rock layer included angles are 0° and 15° are close, their displacements are not much different.

[0097] Figure 14 It is the phase field contour map and the horizontal displacement contour map of the rock layer with an elastic modulus ratio of 4 at t = 200 s in the embodiment of the present invention. As Figure 14 shown, there are significant differences in the propagation pattern of the hydraulic fracture from the simulation results when the ratio of the elastic modulus of the rock layer E1 / E2 = 2. Especially for the models with rock layer included angles of 0° and 15°, because the rock layer strength increases, the interface between the rock layers is not penetrated, and the crack stops propagating when it reaches the interface. There is a slight deflection in the crack propagation of the model with a rock layer included angle of 15°, while there is no sign of deflection in the model with a rock layer included angle of 0°. Figure 15 It is the water pressure curve graph of the rock layer with an elastic modulus ratio E1 / E2 = 4 in the embodiment of the present invention. As Figure 15As shown, the crack starts to expand at point C and reaches the rock formation boundary at point D. When the rock formation angle is 0° and 15°, the crack no longer continues to expand, and the water pressure increases linearly after reaching point D; when the rock formation angle is 30°, the crack continues to expand, so its water pressure is lower than the other two angles, and because of the increase in rock formation strength, from Figure 11 (c) and Figure 14 (c), it can be seen that the propagation speed of the second group of cracks is slower than that of the first group. The propagation requires more energy and a higher water pressure. From Figure 12 and Figure 15 it can be seen that the water pressure value with the elastic modulus ratio E1 / E2 = 2 is significantly lower than the value with the elastic modulus ratio E1 / E2 = 4. Figure 16 is a schematic diagram of the horizontal displacement of the crack surface along the Parh-1 path when the elastic modulus ratio E1 / E2 = 2 in the embodiment of the present invention. From Figure 16 it can be seen that the horizontal displacement of the crack surface along the path Path-1 when the rock formation angle is 30° is less than the displacement when the rock formation angles are 0° and 15°. Because the model cracks with rock formation angles of 0° and 15° no longer continue to expand after reaching the crack surface, their water pressure gradually increases, resulting in a gradual increase in their displacement.

[0098] The present invention also provides a device for simulating the evolution of hydraulic fractures. The device for simulating the evolution of hydraulic fractures provided by the present invention will be described below. The device for simulating the evolution of hydraulic fractures described below can be mutually corresponding and referenced to the method for simulating the evolution of hydraulic fractures described above. Figure 17 is a structural block diagram of the device for simulating the evolution of hydraulic fractures provided by the present invention. As Figure 17 shown, the device includes: A division module 1701, configured to divide the damaged rock mass into several different calculation regions and determine the seepage control equation for each calculation region; A processing module 1702, configured to determine the permeability tensor of the porous medium seepage-stress coupling for each calculation region in combination with seepage anisotropy; An association module 1703, configured to associate the seepage control equations of all calculation regions with the phase field to obtain a target seepage control equation; A determination module 1704, configured to determine the displacement, fluid pressure, and phase field of the damaged rock mass based on the target seepage control equation.

[0099] When the device is in use, first, the division module 1701 divides the damaged rock mass into calculation regions according to the structural components of the damaged rock mass, and determines the seepage control equation for each calculation region. Then, the processing module 1702 derives the permeability tensor of the porous medium seepage-stress coupling considering seepage anisotropy. Then, the correlation module 1703 correlates the seepage control equations of all calculation regions with the phase field to generate a unified target seepage control equation. Finally, the determination module 1704 derives the target seepage control equation and implements the algorithm to obtain the displacement, fluid pressure, and phase field of the damaged rock mass. During this process, a seepage-stress coupling phase field model is derived, the strain energy is decomposed using spectral decomposition to ensure that crack propagation occurs under the energy generated by tensile strain, the influence of water work on crack propagation is considered, and the permeability coefficient changes with the crack opening and rock mass deformation, solving the problem of the limitations of existing research methods for hydraulic fractures.

[0100] Figure 18 An example of the physical structure diagram of an electronic device is shown as Figure 18 shown. The electronic device may include: a processor 1801, a communication interface 1802, a memory 1803, and a communication bus 1804. Among them, the processor 1801, the communication interface 1802, and the memory 1803 communicate with each other through the communication bus 1804. The processor 1801 can call the logical instructions in the memory 1803 to execute the hydraulic fracture evolution simulation method, which includes: Dividing the damaged rock mass into several different calculation regions and determining the seepage control equation for each calculation region; Combining seepage anisotropy to determine the permeability tensor of the porous medium seepage-stress coupling for each calculation region; Correlating the seepage control equations of all calculation regions with the phase field to obtain the target seepage control equation; Based on the target seepage control equation, determining the displacement, fluid pressure, and phase field of the damaged rock mass.

[0101] In addition, when the logical instructions in the above-mentioned memory 1803 are implemented in the form of software functional units and sold or used as independent products, they can be stored in a computer-readable storage medium. Based on such an understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for causing a computer device (which may be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the methods according to the various embodiments of the present invention. The foregoing storage medium includes: various media such as USB flash drives, mobile hard disks, read-only memories (ROMs), random access memories (RAMs), magnetic disks, or optical discs that can store program codes.

[0102] On the other hand, the present invention also provides a computer program product. The computer program product includes a computer program. The computer program can be stored on a non-transitory computer-readable storage medium. When the computer program is executed by a processor, the computer can execute the hydraulic fracture evolution simulation method provided by the above-mentioned various methods. The method includes: Dividing the damaged rock mass into several different calculation regions and determining the seepage control equation for each calculation region; Combining seepage anisotropy to determine the permeability tensor of porous media seepage-stress coupling for each calculation region; Associating the seepage control equations of all calculation regions with the phase field to obtain the target seepage control equation; Based on the target seepage control equation, determining the displacement, fluid pressure, and phase field of the damaged rock mass.

[0103] On another aspect, the present invention also provides a non-transitory computer-readable storage medium, on which a computer program is stored. When the computer program is executed by a processor, it implements the hydraulic fracture evolution simulation method provided by the above-mentioned various methods. The method includes: Dividing the damaged rock mass into several different calculation regions and determining the seepage control equation for each calculation region; Combining seepage anisotropy to determine the permeability tensor of porous media seepage-stress coupling for each calculation region; Associating the seepage control equations of all calculation regions with the phase field to obtain the target seepage control equation; Based on the target seepage control equation, determining the displacement, fluid pressure, and phase field of the damaged rock mass.

[0104] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separated, and the components shown as units may or may not be physical units, that is, they may be located in one place or distributed to multiple network units. Some or all of the modules can be selected according to actual needs to achieve the purpose of the solution of this embodiment. Those of ordinary skill in the art can understand and implement it without creative work.

[0105] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on this understanding, the essence of the above technical solution, or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., and includes several instructions for causing a computer device (which can be a personal computer, server, or network device, etc.) to execute the methods described in each embodiment or some parts of the embodiments.

[0106] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and are not intended to limit them. Although the present invention has been described in detail with reference to the foregoing embodiments, those of ordinary skill in the art should understand that they can still modify the technical solutions recorded in the foregoing embodiments, or perform equivalent replacements for some of the technical features. However, these modifications or replacements do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for simulating the evolution of hydraulic fractures, characterized in that: include: Dividing the damaged rock mass into a number of different calculation areas, and determining the seepage control equation for each of the calculation areas; Determine the permeability tensor of the porous medium seepage-stress coupling in each calculation area in combination with the seepage anisotropy; Associating the seepage control equations of all the calculation areas with the phase field to obtain a target seepage control equation; Based on the target seepage control equation, the displacement, fluid pressure and phase field of the damaged rock mass are determined.

2. The hydraulic fracture evolution simulation method according to claim 1, characterized in that: The damaged rock mass is divided into several different calculation areas, including: dividing the damaged rock mass into a matrix zone, a crack zone and a transition zone; The transition zone is a calculation region located between the matrix zone and the crack zone.

3. The hydraulic fracture evolution simulation method according to claim 2, characterized in that: Determining the seepage control equation of the matrix region includes: Based on Darcy's law and Biot's theory, the storage coefficient and the fluid flow rate of the matrix region are obtained; A seepage control equation of the matrix zone is generated according to the storage coefficient and the fluid flow rate of the matrix zone.

4. The hydraulic fracture evolution simulation method according to claim 2, characterized in that: Determining the permeability tensor of the porous medium seepage-stress coupling in the matrix region, including: obtaining an initial permeability coefficient of the matrix region; Based on the product of the initial permeability coefficient of the matrix zone and the right Cauchy-Green strain tensor, the permeability tensor of the porous medium seepage-stress coupling in the matrix zone is determined.

5. The hydraulic fracture evolution simulation method according to claim 2, characterized in that: Determine the permeability tensor of the porous medium seepage-stress coupling in the crack zone, include, determining a first permeability tensor of the crack zone in an open state and a second permeability tensor in a closed state; The first permeability tensor and the second permeability tensor are combined to obtain the permeability tensor of the seepage-stress coupling of the porous medium in the crack zone.

6. The hydraulic fracture evolution simulation method according to claim 1, characterized in that: Before associating the seepage control equations of all the computational domains with the phase field, including: The permeability tensors of the crack zone and the matrix zone are interpolated to obtain the permeability tensor of the entire calculation area of ​​the damaged rock mass.

7. The hydraulic fracture evolution simulation method according to claim 1, characterized in that: The seepage control equations of all the calculation areas are associated with the phase field to obtain the target seepage control equations, including: Determine the total internal energy of the system in which the damaged rock mass is located; Combined with the total internal energy of the system, according to the variational principle and the energy minimization principle, the phase field control equation of seepage-stress coupling is obtained; Introducing a historical energy field to determine the weak form of each physical field in a seepage-stress coupling model; the seepage-stress coupling model includes a stress field, a phase field and a seepage field; Performing finite element discretization on the weak form of each of the physical fields to obtain the residual of the corresponding physical field; The target seepage control equation is generated by combining the residual of the physical field.

8. The hydraulic fracture evolution simulation method according to claim 7, characterized in that: Based on the target seepage control equation, the displacement, fluid pressure and phase field of the damaged rock mass are determined, including: Solving the displacement and fluid pressure of each physical field in the seepage-stress coupling model by a staggered iteration method, and determining the displacement vector of the node by the fluid pressure of the fixed node until a preset convergence criterion is met; Determine the field variables of the historical energy field through the displacement of the node and the fluid pressure; The field variables of the historical energy field are updated, and the phase field of the damaged rock mass is obtained by the Newton-Raphson method.

9. A hydraulic fracture evolution simulation device, characterized in that: include: A division module, used to divide the damaged rock mass into a number of different calculation areas and determine the seepage control equation of each calculation area; A processing module, used for determining the permeability tensor of the porous medium seepage-stress coupling in each calculation area in combination with the seepage anisotropy; A correlation module, used for correlating the seepage control equations of all the calculation areas with the phase field to obtain a target seepage control equation; A determination module is used to determine the displacement, fluid pressure and phase field of the damaged rock mass based on the target seepage control equation.

10. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the hydraulic fracture evolution simulation method according to any one of claims 1 to 8 is implemented.

Images