Method for solving unsteady seepage of fractured oil and gas reservoir based on physical information neural network

Abstract

This invention discloses a method for solving unsteady seepage in fractured oil and gas reservoirs based on a physical information neural network, belonging to the field of reservoir seepage mechanics technology. The method includes the following steps: S1. Establishing a dual-medium seepage mathematical model: Constructing the governing equations, initial pressure conditions, wellbore flow boundary conditions, and far-field boundary conditions for the fracture system and the matrix system. This method constructs a deep neural network with residual connections and embeds the governing equations, initial conditions, and boundary conditions of the dual-pore, dual-seepage model into the neural network loss function. Through physical constraint learning, it achieves efficient prediction and solution of the seepage field in the fracture system and the matrix system. Compared with traditional numerical solution methods based on grid discretization, this method reduces the dependence on grid partitioning and local refinement, alleviates the problems of high computational cost and low solution efficiency under complex seepage conditions in fractured reservoirs, improves the accuracy of seepage field prediction and computational efficiency, and has good engineering applicability.

Description

Technical Field

[0001] This invention relates to the field of reservoir seepage mechanics, specifically to a method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks. Background Technology

[0002] Fractured oil and gas reservoirs occupy an important position in global oil and gas resources. Their reservoir-permeability system typically exhibits a dual medium consisting of low-permeability matrix pores and a high-conductivity fracture network. In the development of fractured oil and gas reservoirs, accurately and efficiently solving the unsteady flow field in the dual medium of fractures and matrix is ​​crucial for dynamic analysis, production capacity prediction, and development decisions.

[0003] However, in existing technologies, traditional grid-based numerical methods heavily rely on complex grid generation and local refinement techniques to characterize fracture networks and drastic pressure changes in the near-wellbore region when solving such problems. This results in cumbersome preprocessing, a surge in computational degrees of freedom, high computational costs, and low solution efficiency, making it difficult to meet the needs of rapid on-site simulation and real-time analysis. Furthermore, the strong coupling and transient flow between fractures and the matrix system, the nonlinearity of parameters, and the multi-scale nature of the problem easily lead to numerical oscillations and convergence difficulties in traditional numerical solutions. Solution stability is extremely sensitive to the discretization scheme and iteration parameters.

[0004] Furthermore, while Physical Information Neural Networks (PINNs) offer a novel approach to solving partial differential equations without meshes, their direct application to specific engineering problems such as the dual-pore, dual-permeability model of fractured reservoirs—characterized by strong nonlinearity, multi-scale complexity, and intricate boundary conditions—faces significant technical bottlenecks. For instance, deep neural networks are prone to gradient anomalies when training high-dimensional, strongly coupled physical field problems, leading to instability or even training failure. Simultaneously, the highly uneven spatial distribution of physical residuals makes it difficult for conventional sampling strategies to accurately capture complex physical field changes in critical areas such as the near-wellbore region, resulting in insufficient local prediction accuracy. Moreover, balancing multiple physical constraints in the loss function is challenging, making training prone to local optima and resulting in slow overall convergence, failing to meet the accuracy and efficiency requirements for engineering applications.

[0005] Therefore, we propose a solution method for unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks, in order to alleviate or solve the above problems.

[0006] The information disclosed above in this background section is only for enhancing the understanding of the background section of this invention, and therefore may include prior art that is not known to those skilled in the art. Summary of the Invention

[0007] To address the aforementioned technical problems, this invention provides a method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks. This method solves the problems of strong dependence on computational grids, low solution efficiency, poor numerical stability, and unstable training, insufficient local fitting accuracy, and slow convergence when directly applying physical information neural networks in the prior art.

[0008] To achieve the above objectives, this invention provides a method for solving unsteady seepage in fractured reservoirs based on a physical information neural network, comprising the following steps:

[0009] S1. Establish a mathematical model for dual-medium seepage: construct the governing equations, initial pressure conditions, wellbore flow boundary conditions, and far-field boundary conditions for the fracture system and matrix system.

[0010] S2. Construct a multi-output physical information neural network: Build a deep neural network with a residual connection structure, whose inputs are spatial coordinates and time variables, and whose outputs are the pressure of the fracture system and the pressure of the matrix system;

[0011] S3. Construct the physical information loss function: Based on the control equation, initial conditions and boundary conditions, define the partial differential equation residual loss term, initial condition loss term, wellbore boundary loss term and far-field boundary loss term, and then sum the weighted losses of each term to obtain the total loss function.

[0012] S4. Perform adaptive sampling: Different strategies are used to sample training sample points inside the computational domain, at the boundary, near the wellbore, and in the initial time plane. The near-wellbore region uses a logarithmic sampling method based on polar coordinates.

[0013] S5. Calculate higher-order derivatives: Using automatic differentiation technology, calculate the first and second-order partial derivatives of the pressure field with respect to spatial coordinates and time based on the output of the neural network.

[0014] S6. Optimize the training of the neural network: Use the adaptive moment estimation algorithm to minimize the total loss function, implement gradient pruning and dynamic learning rate scheduling during training, and monitor and adjust the weight coefficients of each loss term.

[0015] S7. Result Validation and Visualization: Generate a set of validation points in the computational domain, calculate the error between the pressure predicted by the neural network and the reference solution; plot the pressure distribution map of the crack system and the matrix system, the pressure profile curve along the radial direction, and the training loss change curve.

[0016] A computational system for solving unsteady seepage in fractured reservoirs based on a physical information neural network includes:

[0017] The data preprocessing module is used to generate training and validation sample points and perform adaptive stratified sampling.

[0018] The neural network building block is used to build and initialize multi-output deep neural networks with residual connections.

[0019] The loss function calculation module is used to calculate the components of the physical information loss function and the total loss.

[0020] The automatic differential calculation module is used to calculate the partial derivatives and radial derivatives of the pressure field.

[0021] An optimized training module is used to execute improved optimization algorithms to train neural network parameters;

[0022] The results verification module is used to evaluate the accuracy of the prediction results and generate a visualization chart of pressure distribution and a curve of loss change.

[0023] Compared with the prior art, the beneficial effects of the present invention are:

[0024] This invention employs a meshless physical information neural network solution paradigm, eliminating the need for complex preprocessing mesh generation and local refinement. The trained neural network model can then be directly used for rapid inference and prediction. Verification based on specific embodiments shows that, while achieving equivalent engineering accuracy, this method reduces computation time by approximately 12 times compared to the traditional finite element method, thus meeting the needs of rapid simulation and real-time analysis of seepage fields in oilfields.

[0025] This invention effectively alleviates the gradient vanishing or exploding problem during deep network training by introducing a deep neural network with residual connection structures, ensuring the stability of the optimization process. Simultaneously, considering the large pressure gradient in the near-wellbore region, an adaptive polar coordinate logarithmic sampling strategy and a dynamically weighted loss function are designed to strengthen the physical constraints in key areas, thereby improving the local accuracy of pressure field prediction, especially in the near-wellbore zone.

[0026] The method framework of this invention is specifically designed for the physical characteristics of a dual-pore, dual-permeability model for fractured reservoirs. Its multi-output network structure is naturally suitable for the synchronous prediction of coupled fields, ensuring the physical consistency between the pressure solutions of the fracture system and the matrix system. This method has good adaptability to complex geometric boundaries and is easily accelerated in parallel using hardware such as GPUs.

[0027] The above overview is for illustrative purposes only and is not intended to be limiting in any way. In addition to the illustrative aspects, embodiments, and features described above, further aspects, embodiments, and features of the invention will become readily apparent from the accompanying drawings and the following detailed description. Attached Figure Description

[0028] Figure 1 This is a schematic diagram comparing the pressure distribution of different media along the radial direction of the wellbore in an embodiment of the present invention.

[0029] Figure 2 This is an enlarged schematic diagram of the pressure distribution of different media within a 50-meter radius of the well point in an embodiment of the present invention.

[0030] Figure 3 This is an overall flowchart of the unsteady seepage solution method for fractured oil and gas reservoirs based on physical information neural networks, as implemented in this invention.

[0031] Figure 4 It is a contour map of the spatial distribution of pressure in a crack system obtained by applying the method of this invention.

[0032] Figure 5 It is a contour map of the spatial distribution of pressure in the matrix system obtained by applying the method of this invention.

[0033] Figure 6 It is a contour map of the spatial distribution of the total system pressure obtained by applying the method of this invention.

[0034] Figure 7 This is a graph showing the change of the total loss function value with the training cycle during the neural network training process of this invention.

[0035] Figure 8 This is a graph showing how the various loss components change with the training cycle during the neural network training process of this invention.

[0036] Figure 9 This is a graph showing the dynamic adjustment of the learning rate during the training process of the neural network in this invention. Detailed Implementation

[0037] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. It should be noted that the drawings are schematic and not illustrated to scale. For clarity and convenience, the relative sizes and proportions of the parts shown in the drawings have been exaggerated or reduced in size. Any size is only illustrative and not limiting.

[0038] The overall process of the unsteady seepage solution method for fractured reservoirs of the present invention is as follows: Figure 3 As shown, the steps include defining input parameters, establishing control equations and initial boundary conditions, constructing a PINN neural network, constructing a physical information loss function, adaptive sampling, automatic differentiation, optimization training, result verification, and application output. Each step is interconnected, forming an efficient solution paradigm without mesh dependency.

[0039] Example 1

[0040] A method for solving unsteady seepage in fractured hydrocarbon reservoirs based on physical information neural networks includes the following steps:

[0041] S1. Constructing a mathematical model for a dual-pore, dual-permeability reservoir.

[0042] Establish the governing equations for the fracture system and the matrix system. The governing equation for the fracture system is as follows: The governing equations for the matrix system are: ;in, Indicates the pressure in the fracture system. Indicates the pressure of the matrix system. , These are the compressive conductivity coefficients of the fracture and matrix systems, respectively. This is the crossflow coefficient;

[0043] Set initial and boundary conditions:

[0044] Initial conditions: ,in, This represents the initial formation pressure;

[0045] Wellbore boundary conditions: Based on Darcy's law, at the wellbore radius... satisfy ;

[0046] Far-field boundary conditions: When the distance approaches infinity, , ;

[0047] S2. Construct a multi-output physical information neural network

[0048] Input layer receives spatial coordinates The input features are 4 in total, including time t; the hidden layer contains 6-10 residual blocks, each residual block consists of two fully connected layers and The activation function is composed of a skip connection structure; the output layer contains two independent output heads, which output the pressure of the fracture system respectively. and matrix system pressure The network weights are initialized using the Xavier normal distribution initialization method.

[0049] The activation function of the residual block is the Tanh function; the output of the residual block... With input Satisfying Relationship: In the formula, , This is the weight matrix. , This is the bias vector.

[0050] S3. Construct the physical information loss function

[0051] Residual loss in partial differential equations:

[0052]

[0053] Boundary condition loss:

[0054] Initial condition loss:

[0055] Wellbore condition losses:

[0056] The total loss function is expressed as follows: In the formula, These are adjustable weighting coefficients;

[0057] The weighting coefficient is set to a value of 1. And adjust dynamically according to the training progress.

[0058] S4. Implement adaptive sampling strategy

[0059] Internal point sampling: in the computational domain N points are randomly sampled uniformly within the area. For spatial domain, Maximum simulation time;

[0060] Boundary point sampling: at the boundary of the computational domain The polar coordinate logarithmic sampling method is used, with a focus on sampling the near-wellbore area;

[0061] Wellbore point sampling: on the wellbore surface Enhancing sampling density at the location, with the sampling density decreasing as the distance from the wellbore increases;

[0062] Initial point sampling: Uniform sampling is performed on the initial time plane Ω×{0};

[0063] S5. Calculation of higher-order partial derivatives based on automatic differentiation

[0064] Using the automatic differentiation function of the neural network output, the first and second partial derivatives of the pressure field with respect to spatial coordinates and the first partial derivative with respect to time are calculated; the radial derivative is calculated using the chain rule.

[0065] The radial distance r in the polar coordinate logarithmic sampling method is sampled according to the following formula: In the formula, A random number uniformly distributed in the interval [0,1]. Where is the wellbore radius. The radius of the oil drain is denoted as .

[0066] Wellbore point sampling uses the radial density function: Sampling was performed, among which This is the distance from the wellbore.

[0067] The partial derivatives in automatic differential calculations specifically include:

[0068] First-order spatial partial derivative: , ,

[0069] First-order time partial derivative:

[0070] Second-order spatial partial derivative: , ,

[0071] Radial derivative: ,in

[0072] S6. Multi-strategy collaborative optimization training of neural networks

[0073] The Adam optimizer is used for parameter optimization, with the initial learning rate set to 0.001-0.01; gradient pruning is implemented to limit the gradient norm to no more than a preset threshold; a learning rate scheduling strategy is adopted to dynamically adjust the learning rate according to the training progress; and the loss components are monitored in real time and the loss weight coefficients are dynamically adjusted.

[0074] The learning rate scheduling strategy is one of the following two:

[0075] The StepLR scheduler is used, and the learning rate is multiplied by a decay factor of 0.8 every 1000 training cycles;

[0076] The ReduceLROnPlateau scheduler is used, and the learning rate is halved when the loss function plateaus for more than 500 cycles.

[0077] S7. Result Validation and Visualization

[0078] A test grid is generated within the computational domain, and verification points are uniformly sampled. The relative error between the predicted pressure field and the reference solution is calculated. Contour maps of pressure distribution in the fracture system and matrix system, pressure distribution curves along the radial line, and curves showing the changes in various losses during the training process are plotted.

[0079] Example 2

[0080] This embodiment addresses the unsteady seepage problem in a fractured oil and gas reservoir well. The reservoir exhibits characteristics of a dual-pore, dual-permeability medium. The basic parameters of the reservoir, model, and wellbore are as follows:

[0081] Reservoir physical parameters: initial formation pressure: 60 MPa; bottom hole flowing pressure: 45 MPa; fracture system conductivity: 0.9; matrix system conductivity: 0.01; channeling coefficient: 0.005.

[0082] Model domain parameters: X-direction length: 1000 m, Y-direction length: 1000 m; maximum time: 100 d.

[0083] Wellbore parameters, wellbore location coordinates: X well 500, Y well :500; shaft radius: 0.1m.

[0084] The specific implementation method of this embodiment is the same as that of Embodiment 1, setting training parameters, including 10,000 training generations, 2,000 internal sampling points, and 500 boundary sampling points. The network architecture is designed as follows: an 8-layer residual network is constructed with a hidden layer dimension of 50, initialized using Xavier.

[0085] Training data is generated within the computational domain. Internal points are randomly sampled in the spatial and temporal domains, boundary points are uniformly sampled at the boundaries of the computational domain, wellhead points are sampled in polar coordinates around the wellhead, and initial points are sampled in the initial time plane.

[0086] This embodiment employs a phased training strategy. The changes in loss and learning rate during the training process are shown in the attached figure. Figures 7-9 As shown in the attached document. Figure 7 As shown in the total training loss graph, the total loss value decreases exponentially with the increase of training rounds, and tends to stabilize after 10,000 training rounds, indicating that the model converges well; as shown in the attached figure. Figure 8 As shown in the loss diagrams of each component, the PDE loss, boundary loss, initial loss, and wellbore loss decrease synchronously and eventually remain in equilibrium, proving that the physical constraints are all effectively satisfied.

[0087] As attached Figure 9 As shown in the learning rate change graph, this embodiment uses the StepLR scheduler, and the learning rate gradually decreases with the training cycle, avoiding oscillations in the later stages of training and ensuring the stability of model optimization.

[0088] As attached Figure 1 As shown, the radial pressure distribution diagram clearly illustrates the pressure variations of the fracture system, matrix system, and overall system with distance from the well. The pressure decrease rate of the fracture system is faster, consistent with the dual-medium seepage characteristics of high-permeability fractures. (Attached) Figure 2 Further magnification of the pressure distribution within a 50-meter radius of the well point reveals a steep pressure gradient in the near-well region, verifying the invention's ability to accurately fit strong gradient regions.

[0089] As attached Figure 4 Figure 5 Figure 6 As shown, the pressure contour maps of the fracture system, matrix system, and overall system intuitively reflect the spatial distribution characteristics of pressure. The pressure field gradually recovers to the initial formation pressure from the wellbore to the far field, which is highly consistent with the physical laws of dual-pore dual-permeability flow.

[0090] Therefore, the traditional finite element method takes about 2 hours to solve this problem, while the method of this invention takes only about 10 minutes, improving computational efficiency by about 12 times. The relative error of pressure prediction for the fracture system is less than 2%, and the relative error of pressure prediction for the matrix system is less than 3%. The pressure distribution characteristics are highly consistent with the physical laws of dual-pore dual-permeability flow.

[0091] Example 3

[0092] The computational system for solving the unsteady seepage in fractured reservoirs based on a physical information neural network includes:

[0093] The data preprocessing module is used to generate training and validation sample points and perform adaptive stratified sampling.

[0094] The neural network building block is used to build and initialize multi-output deep neural networks with residual connections.

[0095] The loss function calculation module is used to calculate the components of the physical information loss function and the total loss.

[0096] The automatic differential calculation module is used to calculate the partial derivatives and radial derivatives of the pressure field.

[0097] An optimized training module is used to execute improved optimization algorithms to train neural network parameters;

[0098] The results verification module is used to evaluate the accuracy of the prediction results and generate a visualization chart of pressure distribution and a curve of loss change.

[0099] In summary, compared with traditional grid-based numerical solution methods, this approach eliminates the reliance on grid partitioning and local refinement, and can directly adapt to complex fracture distributions and reservoir geometries. After training, it can quickly predict seepage fields, meeting the real-time dynamic analysis and decision-making needs in engineering scenarios. The neural network possesses excellent interpolation and extrapolation capabilities, adapting to reservoir conditions with varying fracture development levels and matrix heterogeneity. Furthermore, it is suitable for solving multi-field coupled seepage problems and exhibits strong scalability. The highly parallel inference process of the neural network supports GPU acceleration, further improving the computational efficiency of large-scale reservoir models.

[0100] In the description of this specification, the references to terms such as "one embodiment," "some embodiments," "example," "specific example," or "some examples," etc., indicate that a specific feature, structure, material, or characteristic described in connection with that embodiment or example is included in at least one embodiment or example of the present invention. In this specification, the illustrative expressions of the above terms do not necessarily refer to the same embodiment or example. Furthermore, the specific features, structures, materials, or characteristics described may be combined in any suitable manner in one or more embodiments or examples. Moreover, without contradiction, those skilled in the art can combine and integrate the different embodiments or examples described in this specification, as well as the features of different embodiments or examples.

[0101] Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks, characterized in that, Includes the following steps: S1. Constructing a mathematical model for a dual-pore, dual-permeability reservoir. Establish the governing equations for the fracture system and the matrix system. The governing equation for the fracture system is as follows: The governing equations for the matrix system are: ;in, Indicates the pressure in the fracture system. Indicates the pressure of the matrix system. , These are the compressive conductivity coefficients of the fracture and matrix systems, respectively. This is the crossflow coefficient; S2. Construct a multi-output physical information neural network (PINN). Input layer receives spatial coordinates The input features are 4 in total, including time t; the hidden layer contains 6-10 residual blocks, each residual block consists of two fully connected layers and The activation function is composed of a skip connection structure; the output layer contains two independent output heads, which output the pressure of the fracture system respectively. and matrix system pressure The network weights are initialized using the Xavier normal distribution initialization method. S3. Construct the physical information loss function The total loss function is derived from the residual loss in the partial differential equation. Boundary condition loss Initial condition loss and wellbore condition loss Weighted composition, the expression is: In the formula, These are adjustable weighting coefficients; S4. Implement adaptive sampling strategy Internal point sampling: in the computational domain N points are randomly sampled uniformly within the area. For spatial domain, Maximum simulation time; Boundary point sampling: at the boundary of the computational domain The polar coordinate logarithmic sampling method is used, with a focus on sampling the near-wellbore area; Wellbore point sampling: on the wellbore surface Enhancing sampling density at the location, with the sampling density decreasing as the distance from the wellbore increases; Initial point sampling: Uniform sampling is performed on the initial time plane Ω×{0}; S5. Calculation of higher-order partial derivatives based on automatic differentiation Using the automatic differentiation function of the neural network output, the first and second partial derivatives of the pressure field with respect to spatial coordinates and the first partial derivative with respect to time are calculated, and the radial derivative is calculated using the chain rule. S6. Multi-strategy collaborative optimization training of neural networks The Adam optimizer is used for parameter optimization, with the initial learning rate set to 0.001-0.01; gradient pruning is implemented to limit the gradient norm to no more than a preset threshold; a learning rate scheduling strategy is adopted to dynamically adjust the learning rate according to the training progress; and the loss components are monitored in real time and the loss weight coefficients are dynamically adjusted.

2. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that: In step S2, the activation function of the residual block is the Tanh function; the output of the residual block... With input Satisfying Relationship: In the formula, , This is the weight matrix. , This is the bias vector.

3. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that: In step S3, the weighting coefficient is set to... And adjust dynamically according to the training progress.

4. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that: Step S1 further includes setting initial conditions and boundary conditions: Initial conditions: ,in, This represents the initial formation pressure; Wellbore boundary conditions: Based on Darcy's law, at the wellbore radius... satisfy ; Far-field boundary conditions: When the distance approaches infinity, , .

5. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that, In step S4, the radial distance r of the polar coordinate logarithmic sampling method is sampled according to the following formula: In the formula, A random number uniformly distributed in the interval [0,1]. Where is the wellbore radius. The oil drain radius; The wellbore point sampling uses a radial density function: Sampling was performed, among which This is the distance from the wellbore.

6. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that, In step S5, the partial derivatives calculated by the automatic differentiation specifically include: First-order spatial partial derivative: , , First-order time partial derivative: Second-order spatial partial derivative: , , Radial derivative: ,in .

7. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that, In step S6, the learning rate scheduling strategy is one of the following two: The StepLR scheduler is used, and the learning rate is multiplied by a decay factor of 0.8 every 1000 training cycles; The ReduceLROnPlateau scheduler is used, and the learning rate is halved when the loss function plateaus for more than 500 cycles.

8. The method for solving unsteady seepage in fractured oil and gas reservoirs based on physical information neural networks according to claim 1, characterized in that, In step S6, the preset threshold for gradient clipping is 1.0; specifically, the gradient norm of all parameters is calculated, and when it exceeds the threshold of 1.0, the gradient vector is scaled proportionally to make its norm equal to 1.

0.

9. A computational system for implementing the solution method for unsteady seepage in fractured reservoirs based on a physical information neural network as described in any one of claims 1-8, characterized in that, include: The data preprocessing module is used to generate training and validation sample points and perform adaptive stratified sampling. The neural network building block is used to build and initialize multi-output deep neural networks with residual connections. The loss function calculation module is used to calculate the components of the physical information loss function and the total loss. The automatic differential calculation module is used to calculate the partial derivatives and radial derivatives of the pressure field. An optimized training module is used to execute improved optimization algorithms to train neural network parameters; The results verification module is used to evaluate the accuracy of the prediction results and generate a visualization chart of pressure distribution and a curve of loss change.

Images