A method for parallel simulation of hydraulic fracturing fracture propagation

Abstract

To solve the problem of high computational complexity and low efficiency in numerical simulation of hydraulic fracturing, a parallel method for simulating hydraulic fracture propagation is proposed. Discontinuous displacement method (DDM) and finite volume method (FVM) are used to construct three-dimensional fracture grid and calculate stress and fluid pressure distribution in the fracture. By optimizing the DDM stress coefficient matrix generation and combining FVM to calculate the fracture tip pressure to drive the fracture propagation. Using OpenMP multi-thread scheduling and Intel MKL and Eigen library optimization, according to the matrix size, dynamic selection of LU or QR decomposition, so as to improve the calculation efficiency. The method realizes efficient fracture propagation simulation by initializing three-dimensional fracture grid, optimizing stress matrix, constructing fluid-solid coupling equation system, iterative solving, checking convergence and generating three-dimensional fracture propagation model, which provides an efficient and accurate numerical tool for fracturing design of unconventional oil and gas reservoirs.

Description

Technical Field

[0001] This invention relates to the field of parallel computing and numerical simulation technology, and in particular to a fluid-structure interaction parallel numerical simulation method for fracture propagation in hydraulic fracturing of unconventional oil and gas reservoirs. Background Technology

[0002] Hydraulic fracturing is a core technology for developing unconventional oil and gas reservoirs such as shale gas and tight gas. It significantly improves reservoir permeability and oil and gas production by injecting high-pressure fluids into the formation to create artificial fractures. In recent years, multi-well staged fracturing technology has effectively improved operational efficiency and reduced costs by simultaneously or alternately activating multiple fractures in horizontal wells. However, the fracture propagation process in hydraulic fracturing involves strong coupling effects from multiple physics fields, including rock deformation, fluid flow, stress disturbance, and dynamic competition among multiple fractures. Numerical simulation of this process faces challenges such as high computational complexity and long processing times.

[0003] In existing technologies, three-dimensional fluid-structure interaction models based on the Displacement Discontinuity Method (DDM) and the Finite Volume Method (FVM) can accurately describe the crack propagation process, but their computational efficiency is low. Furthermore, existing parallel computing research mostly focuses on simple scenarios and lacks systematic optimization for data interaction and load balancing in large-scale three-dimensional crack propagation, resulting in computational efficiency that fails to meet the needs of real-time engineering decision-making.

[0004] To address the aforementioned issues, this invention proposes a method for parallel simulation of hydraulic fracturing fracture propagation. By optimizing the generation of the elastic coefficient matrix, solving the fluid-structure interaction equations, and employing parallel computing strategies, the computational efficiency and numerical stability of large-scale three-dimensional fracture propagation simulation are significantly improved. Summary of the Invention

[0005] To simplify the cumbersome numerical simulation process of hydraulic fracturing and rapidly generate high-precision 3D fracture propagation models containing details of multi-fracture interactions, this invention combines the advantages of existing methods and designs a parallel simulation method for hydraulic fracturing fracture propagation. This method uses the Displacement Discontinuity Method (DDM) to extract global elastic response features from reservoir geological parameters and injection parameters, and preserves stress coefficient matrices at different scales during fracture mesh discretization. The stress coefficient matrices at different scales are optimized and recombined, and the stress influence coefficients between fracture elements are quickly mapped through queue traversal and mapping sequence tables, generating efficient global stress features. Starting from the most complex mesh of the generated stress coefficient matrices at different scales, matrix partitioning and parallel computation optimization are performed level by level, fusing with the stress features of the previous level mesh to obtain a high-precision stress distribution map during fracture propagation. Simultaneously, the fluid pressure distribution value of each fracture element is calculated using the Finite Volume Method (FVM). The global stress and local fluid characteristics of all crack elements were decoded using linear solvers based on the Intel MKL and Eigen libraries to obtain the global width distribution and local pressure distribution of the cracks, respectively. The stress interference effects of all elements were then added to obtain the implicit representation of the final crack propagation model (solution to the fluid-structure interaction equations). The crack front location was then extracted using the Implicit Level Set Method (ILSM) post-processing, transforming the implicit representation into an explicit three-dimensional crack propagation model.

[0006] Specifically, the method for parallel simulation of hydraulic fracturing fracture propagation provided by the present invention includes the following steps:

[0007] S1. By inputting reservoir geological parameters, fracturing parameters and mesh parameters, an initial three-dimensional fracture mesh is constructed to provide geometric data support for subsequent fluid-structure interaction calculations;

[0008] S2. After generating the initial three-dimensional crack mesh, the elastic stress response of the crack elements is calculated by the displacement discontinuity method (DDM) to generate the stress coefficient matrix and simulate the elastic deformation of the rock mass and the stress shadowing effect between cracks.

[0009] S3. Combining flow distribution and fluid pressure distribution, the pressure at the crack tip is calculated, driving crack propagation. This process simultaneously generates a fluid-structure interaction matrix and provides visualization support for subsequent iterative calculations.

[0010] S4. Solve the fluid-structure interaction equations using ILSM iteration, calculate the crack width and pressure distribution, and monitor convergence to optimize computational efficiency.

[0011] S5. Check the convergence of the fluid-structure interaction equation. If the set threshold is met, advance the time step and save the intermediate results, and continue the iterative calculation.

[0012] S6. Based on the convergence results, extract the fracture front position, generate a three-dimensional fracture propagation model, and visualize key data to provide intuitive support for fracturing design optimization.

[0013] Step S1 constructs an initial three-dimensional fracture mesh by inputting reservoir geological parameters, fracturing parameters, and mesh parameters, providing geometric data support for subsequent fluid-structure interaction calculations. In this process, the Fast Marching Method (FMM) is used to generate the initial three-dimensional fracture mesh, and each element of the fracture mesh is labeled with stress characteristics and location coordinates to ensure that subsequent calculations can accurately obtain the geometric and stress characteristics of the fractures.

[0014] Step S2, after generating the initial 3D crack mesh, calculates the elastic stress response of the crack elements using the displacement discontinuity method (DDM), generating a stress coefficient matrix to simulate the elastic deformation of the rock mass and the stress shadowing effect between cracks. To improve computational efficiency, the following method is employed... Figure 2 The optimization algorithm shown in the flowchart reduces redundant calculations by traversing the spatial locations of crack surfaces and crack elements, ensuring that the generation of the stress coefficient matrix is ​​efficient and accurate.

[0015] Step S3, based on the flow distribution and fluid pressure distribution in step S2, calculates the pressure at the crack tip and uses this pressure to drive crack propagation. The construction of the fluid-structure interaction matrix combines fluid pressure with crack stress characteristics, and multi-threaded parallel computation improves computational efficiency. The results are displayed through a visualization system, providing intuitive support for subsequent crack propagation and numerical solutions.

[0016] Step S4 uses ILSM iterative solution to solve the fluid-structure interaction equations, calculates the crack width and pressure distribution, dynamically adjusts the crack propagation path, and optimizes computational efficiency. This implicit method not only ensures computational accuracy but also improves the speed and stability of the solution process through multi-threaded parallel processing, enabling efficient operation of large-scale crack propagation simulations.

[0017] Step S5 checks the convergence of the fluid-structure interaction equations after each iteration. If the convergence condition meets a preset threshold (e.g., residual less than 1.5%), the crack width and pressure distribution are updated, and the time step is advanced to the next calculation. If convergence is not achieved, iterative calculation continues. Multi-threaded optimization ensures computational efficiency and accuracy for each iteration.

[0018] Step S6 extracts the fracture front location based on the convergence results, generates a three-dimensional fracture propagation model, and visualizes key data. The model provides a dynamic evolution of fracture propagation, helping engineers intuitively understand the morphology and evolution process of fractures, and providing a reliable basis for fracturing design optimization. Attached Figure Description

[0019] Figure 1 Flowchart for the numerical model calculation of hydraulic fracture propagation;

[0020] Figure 2 Generate a stress coefficient matrix for crack propagation;

[0021] Figure 3 Optimize the flowchart for DDM calculation process;

[0022] Figure 4 Generate an optimization plot for the stress matrix;

[0023] Figure 5 Construct an optimization graph for the stress matrix;

[0024] Figure 6 A schematic diagram of a multi-well fracturing model and its mesh discretization;

[0025] Figure 7 The experimental results are shown in the figure. Detailed Implementation

[0026] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the embodiments and accompanying drawings. Here, the illustrative embodiments and their descriptions are used to explain the invention, but are not intended to limit the invention.

[0027] like Figures 1 to 7 As shown, the present invention provides a method for parallel simulation of hydraulic fracturing fracture propagation, comprising the following steps:

[0028] S1. The purpose of the initialization system is to construct an initial fracture mesh from reservoir geological parameters and fracturing parameters. By inputting reservoir geological parameters (Young's modulus 27000, Poisson's ratio 0.25, fracture toughness 1.0), fracturing parameters (injection rate 0.1, fracturing fluid viscosity 0.001), and mesh parameters (mesh size 5, computational domain size 2000×2000×30), mesh elements are generated. A three-dimensional fracture mesh is generated using the Fast Marching Method (FMM), and the mesh is traversed, with elements identified as non-fracture surfaces, fracture boundaries, adjacent fracture boundaries, or fracture interiors, thereby constructing an initial three-dimensional fracture computational domain model. The calculation steps are as follows: Figure 1The flowchart shown lays the foundation for subsequent iterative calculations of the crack mesh. The initial 3D crack mesh uses a geometric mapping algorithm to determine the crack element coordinates and initial crack width, providing accurate geometric data for flow distribution. OpenMP multi-threaded scheduling technology is utilized to significantly improve mesh generation efficiency.

[0029] S2. Based on the crack mesh generated in step S1, the displacement discontinuity method (DDM) is used to calculate the elastic stress response between crack elements, generating a stress coefficient matrix to efficiently simulate the elastic deformation of rock mass during multi-crack propagation. For example... Figure 2 As shown, the algorithm first traverses the crack propagation mesh structure, extracts the spatial location indices of crack surfaces, tip elements, and non-tip elements, establishes a crack interaction queue, and forms a preliminary partition of the stress coefficient matrix, including the coupling relationships between crack surfaces, tips and non-tips, and non-tips and non-tips. To avoid redundant calculations in traditional methods, this invention introduces an optimization process, such as... Figure 3 As shown in the flowchart. Figure 4 As shown, the optimization algorithm traverses the discrete mesh elements on the crack surface, adding elements to the mapping queue according to the expansion direction (e.g., from the crack front to the interior). During iteration, the first crack element in the queue is retrieved sequentially, and its status is checked: if unprocessed, its number is recorded at the end of the mapping order list, and its stress influence coefficient with all mapped elements is calculated based on the three-dimensional elasticity integral form (using the Green's function kernel and displacement discontinuity vector). This generates only one column in the upper triangular matrix, reducing redundant calculations by approximately 50%. If processed, it is skipped to continue processing the next element, looping until the queue is empty. After generating the upper triangular stress coefficient matrix, further fine classification and query optimization are performed. For example... Figure 5 As shown, by traversing the fracture structure surface, pointed elements (crack front), non-pointed internal elements, and complete crack surface elements are extracted and identified, establishing a classification queue, and combining it with a mapping sequence table to achieve efficient indexing. For example, the non-pointed-non-pointed coefficient (such as C...) 3,3 ) directly from the upper triangular matrix (e.g., A) using a sequential list 1,1 Mapping extraction is used, and the coefficients of the tip and non-tip are quickly located through position offset. This optimization reduces memory usage by 30%-50%, results in a compact data structure, and allows for independent computational tasks. The tasks are allocated to multi-core CPUs (such as a 6-core environment) via OpenMP multi-threaded scheduling, achieving parallel acceleration. Stress calculation is based on the following elasticity equation:

[0030]

[0031] Where C represents the elastic coefficient matrix, w is the fracture width vector of the discrete grid, p is the fluid pressure vector of the discrete grid, and σ is the far-field stress vector. This equation efficiently simulates the elastic response and stress shadowing effect of rock mass by calculating the stress influence between fracture elements. The inlet flow rate of each fracture cluster is dynamically adjusted using the Newton-Raphson iteration method to accurately simulate the interaction between multiple fractures. The fluid pressure calculation unit of the data processing system is based on the finite volume method (FVM) to calculate the initial fluid pressure distribution, providing reliable input for subsequent pressure iterations. Parallel optimization is coordinated by the computer CPU, utilizing OpenMP multi-threaded scheduling and the Intel MKL library to significantly improve matrix generation efficiency and provide efficient support for real-time fracturing design optimization.

[0032] S3. By traversing the fracture surface and its external adjacent boundaries, and combining the flow distribution results from step S2 with the fluid pressure distribution in the fracture tip region, the fracture tip pressure is calculated to drive fracture propagation. The viscous fluid flow within the fracture follows Poiseuille's law, and the fluid velocity is expressed as:

[0033]

[0034] in, v The flow velocity of the fracturing fluid within the fracture. w The width of the crack. p For fluid pressure, μ For liquid viscosity, The pressure gradient is along the crack path. The mass balance equation based on the finite volume method is:

[0035]

[0036] in, This represents the rate of change of crack width over time. For fluid flux divergence, q This section describes the fracturing fluid injection process. Lubrication theory equations are used to simulate fluid flow within the fracture, comprehensively considering fracturing fluid viscosity and leakage effects to provide accurate pressure input for calculating the elastic response of the rock mass. The stress matrix generated by DDM and FVM is combined with the AVX instruction set and optimized based on the L1 buffer size to construct the fluid-structure interaction matrix. The calculation process is monitored in real time to ensure computational stability. A visualization of the pressure distribution at the fracture tip is generated, providing intuitive support for fracture propagation analysis.

[0037] S4. Using the Implicit Level Set Method (ILSM), the fluid-structure interaction equations constructed in parallel in step S3 are solved simultaneously, iteratively calculating the location, width, and fluid pressure of the crack front. The fluid-structure interaction matrix equations are expressed as follows:

[0038]

[0039] In this context, matrices A, B, C, and D represent the coupling relationship between the elastic stress coefficient, fluid pressure, and crack width among the crack elements, respectively. w Let the crack width vector be... p For fluid pressure vectors, f and g These are external stress and fluid injection terms, respectively. Through a gating mechanism, the computational efficiency is significantly improved by dynamically invoking Intel MKL's LU decomposition or Eigen's QR decomposition based on the matrix size and ill-conditioned nature, combined with OpenMP multi-threaded scheduling.

[0040] S5. The numerical solution unit in step S5 checks the convergence of the fluid-structure interaction equations based on the solution results of step S4. If the residual is less than a preset threshold (0.5%-1.5%), the crack width is updated and the process moves to the next time step; if convergence fails, the iterative solution continues. The system displays the convergence status in real time, facilitating monitoring of the computation process. Parallel optimization dynamically adjusts thread allocation, combined with OpenMP multi-threaded scheduling, to improve iterative computation efficiency. Intermediate results such as crack width and fluid pressure are saved through data input and output, preparing for the final result output and analysis.

[0041] S6. Using the MATLAB image library, and combining the Fast Marching Model (FMM) and tip asymptotic solution, extract the fracture front position based on the convergence results of step S5, generating a three-dimensional fracture propagation model. Output key data such as fracture width, pressure distribution, and propagation path. The display unit generates a visualization of the three-dimensional fracture propagation, intuitively presenting the fracture morphology and dynamic evolution. Parallel optimization of data storage and output efficiency improves large-scale data processing performance. Coordinated result verification ensures a high degree of consistency between simulation results and actual operating conditions (fracture front offset difference controlled within 5%), providing a reliable basis for fracturing design optimization.

[0042] In a specific embodiment of this invention, the data used originates from a shale gas block in Southwest China. Using preset geometric mapping parameters, a fracture model is projected and rendered to obtain corresponding two-dimensional stress and pressure distribution maps. Simultaneously, the position coordinates and width values ​​of spatial points near the fracture surface are saved as reference data. In this specific embodiment, the computational model employs a C++-based numerical simulation program, integrating Intel MKL and Eigen libraries. The optimizer is OpenMP multi-threaded scheduling, and the computational environment is a high-performance server with a 6-core CPU.

[0043] The corresponding data, calculation process diagrams, and simulation result diagrams are as follows: Figure 7As shown, the color bars on the right side of the figure represent the width of the crack, with different colors corresponding to different crack width values; the red boundary line clearly marks the boundary position of the crack. Compared with previous methods, the simulation method of this invention has a better performance in the detailed parts of the crack surface and the stress interference parts.

[0044] The technical means disclosed in this invention are not limited to those disclosed in the above embodiments, but also include technical solutions composed of any combination of the above technical features. It should be noted that those skilled in the art can make various improvements and modifications without departing from the principles of this invention, and these improvements and modifications are also considered within the scope of protection of this invention.

Claims

1. A method for parallel simulation of hydraulic fracturing fracture propagation, characterized in that, Includes the following steps: S1: By inputting reservoir geological parameters, fracturing parameters and mesh parameters, an initial three-dimensional fracture mesh is constructed to provide basic geometric data for subsequent fluid-structure interaction calculations; S2: Based on the initial 3D fracture mesh generated in step S1, the queue traversal method effectively avoids redundant mesh calculations, the mapping sequence table is used to realize the orderly construction of the stress coefficient triangular matrix, and the overall generation process of the stress coefficient matrix is ​​accelerated by the team query method, thereby significantly improving the computational efficiency and storage utilization. The inlet flow rate of each cluster of fractures is dynamically adjusted by the Newton iteration method to accurately simulate the stress shadow effect between multiple fractures, and the boundary element displacement discontinuity method (DDM) is used to calculate the elastic response of the rock mass and the initial fluid pressure distribution, providing data support for subsequent fluid-structure interaction calculations; S3: Based on the flow distribution and initial fluid pressure distribution in step S2, the fluid pressure calculation unit calculates the pressure at the fracture tip to drive fracture propagation by combining the finite volume method (FVM) with the lubrication theory equation, taking into account the fracturing fluid viscosity and leakage effect. It constructs the fluid-structure interaction matrix and outputs the pressure distribution visualization results to support subsequent iterative calculations. S4: Based on the fluid-structure interaction matrix constructed in step S3, the implicit level set method (ILSM) is used for iterative solution to calculate the crack width and pressure distribution, optimize computational efficiency and monitor convergence in real time; S5: Based on the solution results of the fluid-structure interaction equations in step S4, check the convergence. If the residual is less than the preset threshold of 0.5%-1.5%, update the crack width and advance the time step iteration; otherwise, continue the iteration. S6: Based on the convergence results of step S5, extract the fracture front position, generate a three-dimensional fracture propagation model, output the fracture width, pressure distribution and propagation path, and visualize them to support fracturing design optimization.

2. The method for parallel simulation of hydraulic fracturing fracture propagation according to claim 1, characterized in that, The reservoir geological parameters include Young's modulus of 27,000, Poisson's ratio of 0.25, and fracture toughness of 1.0; the fracturing parameters include injection rate of 0.1 and fracturing fluid viscosity of 0.001; the grid parameters include grid size of 5 and computational domain range of 2000×2000×30.

3. The method for parallel simulation of hydraulic fracturing fracture propagation according to claim 1, characterized in that, The initial three-dimensional crack mesh determines the crack element coordinates and initial crack width through a geometric mapping algorithm, and identifies the crack mesh elements as non-crack surfaces, crack boundaries, crack adjacent boundaries, or crack interiors, thus constructing an initial three-dimensional crack computational domain model.

4. The method for parallel simulation of hydraulic fracturing fracture propagation according to claim 1, characterized in that, The construction of the fluid-structure interaction matrix utilizes the AVX instruction set and optimizes the construction of the fluid-structure interaction matrix based on the L1 cache size, thereby improving the efficiency of large-scale matrix operations.

5. The method for parallel simulation of hydraulic fracturing fracture propagation according to claim 1, characterized in that, In step S4, the gating mechanism dynamically calls the LU decomposition of Intel MKL or the QR decomposition of Eigen library according to the matrix size and ill-conditioned nature, and combines OpenMP multi-threaded scheduling to improve the efficiency of iterative solution.

6. The method for parallel simulation of hydraulic fracturing fracture propagation according to claim 1, characterized in that, In step S6, the display unit calls the MATLAB image library and combines the Fast Marching Method (FMM) and tip asymptotic solution to generate a three-dimensional crack propagation model, ensuring that the difference in crack leading edge offset is controlled within 5%.

7. The method for parallel simulation of hydraulic fracturing fracture propagation according to claim 1, characterized in that, The method is applicable to the fracturing optimization design of unconventional oil and gas reservoirs such as shale gas and tight gas.

Images