Method for parallel simulation of hydraulic fracturing crack propagation

Abstract

In order to solve the problems of high calculation complexity and low efficiency in hydraulic fracturing numerical simulation, the invention provides a method for parallel simulation of hydraulic fracturing crack propagation. A three-dimensional fracture grid is constructed by adopting a displacement discontinuity method (DDM) and a finite volume method (FVM), and stress and fluid pressure distribution in fractures are calculated. Crack expansion is driven by optimizing DDM stress coefficient matrix generation and combining FVM to calculate crack tip pressure. OpenMP multi-thread scheduling and Intel MKL and Eigen library optimization are adopted, LU or QR decomposition is dynamically selected according to the matrix scale, and therefore the calculation efficiency is improved. According to the method, by initializing a three-dimensional fracture grid, optimizing a stress matrix, constructing a fluid-solid coupling equation set, carrying out iterative solution, checking convergence and generating a three-dimensional fracture expansion model, efficient fracture expansion simulation is achieved, and an efficient and accurate numerical tool is provided for fracture design of unconventional oil and gas reservoirs.

Description

TECHNICAL FIELD

[0001] The present application relates to the field of parallel computing and numerical simulation, and particularly relates to a fluid-solid coupling parallel numerical simulation method for hydraulic fracture propagation of unconventional oil and gas reservoirs. BACKGROUND

[0002] Hydraulic fracturing is a core technology for developing unconventional oil and gas reservoirs such as shale gas and tight gas, which forms artificial fractures by injecting high-pressure fluid into the formation, significantly improving reservoir permeability and oil and gas production. In recent years, multi-well staged fracturing technology has effectively improved operation efficiency and reduced costs by simultaneously or alternately activating multiple fractures in horizontal wells. However, the hydraulic fracture propagation process involves strong coupling of multiple physical fields such as rock deformation, fluid flow, stress interference, and multi-fracture dynamic competition, and its numerical simulation faces challenges such as high computational complexity and long time consumption.

[0003] In the prior art, three-dimensional fluid-solid coupling models based on the displacement discontinuity method (DDM) and the finite volume method (FVM) can accurately describe the fracture propagation process, but the computational efficiency is low. In addition, existing parallel computing researches are mostly focused on simple scenarios, and lack of systematic optimization of data interaction and load balancing in large-scale three-dimensional fracture propagation, resulting in difficulty in meeting the real-time decision-making needs of engineering with the computational efficiency.

[0004] To solve the above problems, the present application proposes a parallel simulation method for hydraulic fracture propagation, which significantly improves the computational efficiency and numerical stability of large-scale three-dimensional fracture propagation simulation by optimizing the generation of elastic coefficient matrix, the solution of fluid-solid coupling equation system and parallel computing strategy. SUMMARY

[0005] In order to simplify the cumbersome numerical simulation process of hydraulic fracturing and quickly generate a high-precision three-dimensional fracture propagation model containing the details of the interaction of multiple fractures, the present application combines the advantages of existing methods and designs a method for parallel simulation of hydraulic fracturing fracture propagation. The method uses the displacement discontinuity method (DDM) to extract the global elastic response characteristics from the reservoir geological parameters and injection parameters, and saves the stress coefficient matrices of different scales in the fracture grid discretization process. The stress influence coefficient between the fracture elements is quickly mapped through queue traversal and mapping order table by optimizing and recombining the stress coefficient matrices of different scales, so as to generate efficient global stress characteristics. Starting from the most complex grid of the generated stress coefficient matrix of different scales, the matrix block and parallel calculation optimization are carried out step by step, and the stress characteristics of the previous grid are fused to obtain the high-precision stress distribution map in the fracture propagation process. At the same time, the fluid pressure distribution value of each fracture element is calculated by the finite volume method (FVM). The global stress characteristics and local fluid characteristics of all fracture elements are decoded into the global width distribution and local pressure distribution of the fracture by using the linear solver based on Intel MKL and Eigen library, and the stress interference effect of all elements is added to obtain the implicit representation (fluid-structure coupling equation set) of the final fracture propagation model. The fracture front position is extracted through the implicit level set method (ILSM) post-processing, and the implicit representation is converted into an explicit three-dimensional fracture propagation model.

[0006] Specifically, the method for parallel simulation of hydraulic fracturing fracture propagation provided by the present application comprises the following steps: S1, by inputting reservoir geological parameters, fracturing parameters and grid parameters, an initial three-dimensional fracture grid is constructed to provide geometric data support for subsequent fluid-structure coupling calculation; S2, after generating the initial three-dimensional fracture grid, the elastic stress response of the fracture element is calculated by the displacement discontinuity method (DDM) to generate a stress coefficient matrix, simulating the elastic deformation of the rock mass and the stress shadow effect between the fractures; S3, combined with flow distribution and fluid pressure distribution, the crack tip pressure is calculated, and the crack propagation is driven. This process generates a fluid-structure coupling matrix, and provides support for subsequent iterative calculation through visual results; S4, the fluid-structure coupling equation is solved iteratively by ILSM to calculate the crack width and pressure distribution, and the convergence is monitored to optimize the calculation efficiency; S5, the convergence of the fluid-structure coupling equation is checked, if the set threshold is met, the time step is advanced and the intermediate results are saved, and the iterative calculation is continued; S6, extracting the crack front position based on the convergence result, generating a three-dimensional crack propagation model, and visualizing key data to provide intuitive support for fracturing design optimization.

[0007] Step S1 builds an initial three-dimensional crack grid by inputting reservoir geological parameters, fracturing parameters and grid parameters, providing geometric data support for subsequent fluid-solid coupling calculation. In this process, the initial three-dimensional crack grid is generated using the Fast Marching Method (FMM), and each unit of the crack grid is identified with stress characteristics and position coordinates to ensure that the subsequent calculation can accurately obtain the geometric and stress characteristic information of the crack.

[0008] Step S2, after generating the initial three-dimensional crack grid, calculates the elastic stress response of the crack unit by the displacement discontinuity method (DDM) to generate a stress coefficient matrix, simulating the elastic deformation of the rock mass and the stress shadow effect between cracks. To improve the calculation efficiency, the optimization algorithm shown in the flow chart is used to reduce redundant calculations by traversing the spatial positions of the crack surface and crack unit, ensuring that the stress coefficient matrix is generated efficiently and accurately. Figure 2 The flow chart shows the optimization algorithm, which reduces redundant calculations by traversing the spatial positions of the crack surface and crack unit, ensuring that the stress coefficient matrix is generated efficiently and accurately.

[0009] Step S3 calculates the pressure at the crack tip based on the flow distribution and fluid pressure distribution in step S2, and drives the crack propagation through these pressures. The construction of the fluid-solid coupling matrix combines fluid pressure with crack stress characteristics, and improves the calculation efficiency through multi-thread parallel calculation. The results are displayed through a visualization system to provide intuitive support for subsequent crack propagation and numerical solution.

[0010] Step S4 solves the fluid-solid coupling equation iteratively through ILSM to calculate the crack width and pressure distribution, dynamically adjusts the crack propagation path and optimizes the calculation efficiency. The implicit method not only ensures the accuracy of the calculation, but also improves the speed and stability of the solution process through multi-thread parallel processing, enabling large-scale crack propagation simulation to run efficiently.

[0011] Step S5 checks the convergence of the fluid-solid coupling equation set after each iteration calculation. When the convergence condition meets the preset threshold (such as residual less than 1.5%), the crack width and pressure distribution will be updated, and the time step will be advanced for the next calculation. If the convergence does not meet, continue the iteration calculation. Through multi-thread optimization, the calculation efficiency and accuracy of each iteration are ensured.

[0012] Step S6 extracts the crack front position based on the convergence result, generates a three-dimensional crack propagation model, and visualizes key data. The model provides a dynamic evolution of crack propagation, helping engineers intuitively understand the morphology and evolution process of the crack, and providing reliable basis for fracturing design optimization. BRIEF DESCRIPTION OF DRAWINGS

[0013] Figure 1 Flow chart for hydraulic fracture propagation numerical model calculation; Figure 2 Stress coefficient matrix generation for fracture propagation; Figure 3 Optimization flow chart for DDM calculation process; Figure 4 Optimization chart for stress matrix generation; Figure 5 Optimization chart for stress matrix construction; Figure 6 Schematic diagram for multi-well fracturing model and grid discretization; Figure 7 Experimental results chart; DETAILED DESCRIPTION

[0014] To make the objectives, technical solutions, and advantages of the present application clearer, further detailed explanations of the present application are provided below in conjunction with embodiments and drawings. Here, the illustrative embodiments of the present application and their explanations are used to explain the present application, but are not intended to limit the present application.

[0015] As shown in Figures 1 to 7 The present application provides a method for parallel simulation of hydraulic fracturing fracture propagation, comprising the following steps: S1, the purpose of initializing the system is to construct an initial fracture grid 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 grid parameters (grid size 5, calculation domain size 2000x2000x30), a grid cell is generated. A fast marching method (FMM) is used to generate a three-dimensional fracture grid, and the grid is traversed to identify the cell as a non-fracture surface, a fracture boundary, a fracture adjacent boundary or a fracture interior, thereby constructing an initial three-dimensional fracture calculation domain model, the calculation steps of which are shown in the flow chart Figure 1 The initial three-dimensional fracture grid uses a geometric mapping algorithm to determine the fracture cell coordinates and the initial fracture width, providing accurate geometric data for flow distribution. The use of OpenMP multi-thread scheduling technology significantly improves the efficiency of grid generation.

[0016] S2, based on the fracture grid generated in step S1, the displacement discontinuity method (DDM) is used to calculate the elastic stress response between fracture cells, generate a stress coefficient matrix, and efficiently simulate the elastic deformation of rock mass in multi-fracture propagation. As shown in Figure 2As shown, the algorithm first traverses the crack propagation grid structure, extracts the spatial position index of the crack surface, tip element and non-tip element, establishes the interaction queue between cracks, and forms the preliminary division of the stress coefficient matrix, including the coupling relationship between crack surface and crack surface, tip and non-tip, and non-tip and non-tip. In order to avoid redundant calculation in the traditional method, the present application introduces an optimization process, such as Figure 3 As shown in the flow chart. As Figure 4 shown, the optimization algorithm traverses the discrete grid elements on the crack surface, and adds the elements to the mapping queue according to the propagation direction (such as from the crack front to the interior); in the iteration, the first crack element in the queue is taken out in turn, and it is judged whether it has been processed: if not, record its number to the tail of the mapping order table, calculate its stress influence coefficient with all the mapped elements based on the integral form of three-dimensional elasticity mechanics (using Green function kernel and displacement discontinuity vector), only generate a new column of upper triangular matrix, reduce about 50% of repeated calculation; if it has been processed, skip to continue processing the next element, and loop until the queue is empty. After generating the upper triangular stress coefficient matrix, further fine classification and query optimization are carried out. As Figure 5 shown, by traversing the crack structure surface, extracting and identifying tip elements (crack front), non-tip internal elements and complete crack surface elements, a classification queue is established, and efficient indexing is realized combined with the mapping order table. For example, non-tip-non-tip coefficient (such as C 3,3 ) is directly extracted from the upper triangular matrix (such as A 1,1 ) through the order table, and the tip and non-tip coefficient is quickly located through the position offset. This optimization reduces the memory occupation by 30%-50%, the data structure is compact, and the calculation task space is independent, which is allocated to multi-core CPU (such as 6-core environment) through OpenMP multi-thread scheduling mechanism, realizing parallel acceleration. The stress calculation is based on the following elastic equation: Where C represents the elastic coefficient matrix, w is the crack 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 shadow effect of rock mass by calculating the stress influence between crack elements. Newton iteration method is used to dynamically adjust the inlet flow of each cluster of cracks, accurately simulating the interaction between multiple cracks. The fluid pressure calculation unit of the data processing system is based on the finite volume method (FVM), which calculates the initial fluid pressure distribution to provide reliable input for subsequent pressure iteration. The computer CPU coordinates parallel optimization, calls OpenMP multi-thread scheduling and Intel MKL library, significantly improves the efficiency of matrix generation, and provides efficient support for real-time fracturing design optimization.

[0017] S3, calculate the crack tip pressure by traversing the crack surface and its external adjacent boundary, combining the flow distribution result of step S2 and the fluid pressure distribution of the crack tip region, and driving the crack propagation. The viscous fluid flow in the crack follows the Poiseuille law, and the fluid flow rate is represented as: wherein, v is the in-crack flow rate of the fracturing fluid, w is the crack width, p is the fluid pressure, μ is the liquid viscosity, is the pressure gradient along the crack path. The mass balance equation based on the finite volume method is: wherein, represents the time-varying rate of crack width, is the fluid flux divergence, q is the fracturing fluid injection term. The lubrication theory equation is used to simulate the fluid flow in the crack, considering the fracturing fluid viscosity and leakage effect, and providing accurate pressure input for rock mass elastic response calculation. The stress matrix generated by combining DDM and FVM is used to construct the fluid-solid coupling matrix using AVX instruction set and according to L1 cache size optimization. Real-time monitoring of the calculation process ensures calculation stability. The visualization result of the crack tip pressure distribution is generated, providing intuitive support for crack propagation analysis.

[0018] S4, use the implicit level set method (ILSM) to solve the fluid-solid coupling equation system constructed in parallel in step S3, and iteratively calculate the crack front position, width and fluid pressure. The fluid-solid coupling matrix equation system is represented as: wherein, matrices A, B, C, D represent the elastic stress coefficients between crack elements, the coupling relationship between fluid pressure and crack width, w is the crack width vector, p is the fluid pressure vector, f and g are the external stress and fluid injection term, respectively. Through the gating mechanism, according to the matrix size and ill-conditionedness, dynamically call the LU decomposition of Intel MKL or the QR decomposition of Eigen library, combined with OpenMP multi-thread scheduling, significantly improve the calculation efficiency.

[0019] S5, the numerical solving unit of step S5 checks the convergence of the fluid-structure coupling equation set based on the solving result of step S4. If the residual is less than a preset threshold (0.5%-1.5%), the crack width is updated and advanced to the next time step; if it is not converged, the iterative solving is continued. The convergence state is prompted in real time through a display system, facilitating the monitoring of the calculation process. The thread allocation is dynamically adjusted in parallel optimization, combined with OpenMP multi-thread scheduling, to improve the efficiency of iterative calculation. The intermediate results such as crack width and fluid pressure are saved through data input and output, preparing for the final result output and analysis.

[0020] S6, the MATLAB image library is called, the crack front position is extracted based on the convergence result of step S5 by combining the fast marching method (FMM) and the tip asymptotic solution, and a three-dimensional crack propagation model is generated. The key data such as crack width, pressure distribution and propagation path are output. The display unit generates a visual result of three-dimensional crack propagation, which intuitively presents the crack morphology and dynamic evolution. The data storage and output efficiency are optimized in parallel, improving the performance of large-scale data processing. The result is verified to ensure that the simulation result is highly consistent with the actual working condition (the crack front displacement difference is controlled within 5%), providing a reliable basis for fracturing design optimization.

[0021] In a specific implementation case of the present application, the data used is derived from a shale gas block in the southwest region. By presetting the geometric mapping parameters, the corresponding two-dimensional stress distribution map and pressure distribution map are obtained by projecting and rendering the crack model, and the position coordinates and width values of the space points near the crack surface are saved as reference data. In the specific implementation case of the present application, the calculation model uses a numerical simulation program based on C++, integrates Intel MKL and Eigen library, the optimizer is OpenMP multi-thread scheduling, and the calculation environment is a 6-core CPU high-performance server.

[0022] The corresponding data, calculation process diagram and simulation result diagram are shown in Figure 7 The right color bar in the figure represents the width of the crack, and the change of color corresponds to different crack width values; the red boundary line clearly indicates the boundary position of the crack. Compared with the previous method, the simulation method of the present application has more excellent performance effect in the detail part of the crack surface and the stress interference part.

[0023] The technical means disclosed in the present application scheme is not limited to the technical means disclosed in the above-mentioned embodiments, but also includes the technical solutions composed of any combination of the above technical features. It should be noted that, for ordinary skilled persons in the art, without departing from the principle of the present application, some improvements and refinements can be made, which are also considered as the protection scope of the present application.

Claims

1. A method for parallel simulation of hydraulic fracturing crack expansion, characterized in that: The following steps are involved: S1: By inputting reservoir geological parameters, fracturing parameters and grid parameters, an initial 3D fracture grid is constructed to provide basic geometric data for subsequent fluid-solid coupling calculations; S2: Based on the initial three-dimensional fracture grid generated in step S1, optimize the stress coefficient matrix generation, calculate the rock mass elastic response and initial fluid pressure distribution, simulate the stress shadow effect between multiple fractures, and provide data support for subsequent fluid-solid coupling calculations; S3: Based on the flow distribution and initial fluid pressure distribution in step S2, the fracture tip pressure is calculated to drive fracture propagation, a fluid-solid coupling matrix is ​​constructed, and a pressure distribution visualization result is output 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 the calculation efficiency and monitor the convergence in real time; S5: Based on the solution results of the fluid-structure coupling 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, the fracture front position is extracted, a three-dimensional fracture propagation model is generated, and the fracture width, pressure distribution, and propagation path are output and visualized to provide support for fracturing design optimization.

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

3. The method for parallel simulation of hydraulic fracturing crack propagation according to claim 1, characterized in that: The initial three-dimensional crack grid determines the crack unit coordinates and initial crack width through a geometric mapping algorithm, and identifies the crack grid units as non-crack surfaces, crack boundaries, crack adjacent boundaries or crack interiors to construct an initial three-dimensional crack calculation domain model.

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

5. The method for parallel simulation of hydraulic fracturing crack propagation according to claim 1, characterized in that: In step S2, queue traversal and mapping sequence table are used to optimize the generation of stress coefficient matrix, and the inlet flow of each cluster of cracks is dynamically adjusted through Newton iteration method to accurately simulate the stress shadow effect between multiple cracks.

6. The method for parallel simulation of hydraulic fracturing crack propagation according to claim 1, characterized in that: In step S3, the fluid pressure calculation unit calculates the fracture tip pressure based on the finite volume method (FVM) combined with the lubrication theory equation, taking into account the viscosity of the fracturing fluid and the leakage effect.

7. The method for parallel simulation of hydraulic fracturing crack propagation according to claim 1, characterized in that: In step S4, the LU decomposition of Intel MKL or the QR decomposition of the Eigen library is dynamically called according to the matrix size and pathology through a gating mechanism, combined with OpenMP multi-threaded scheduling to improve the iterative solution efficiency.

8. The method for parallel simulation of hydraulic fracturing crack 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 the tip asymptotic solution to generate a three-dimensional crack propagation model to ensure that the offset difference of the crack front is controlled within 5%.

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

Images