Physical Information Neural Network Numerical Simulation Method Applicable to Low-Permeability Waterflood Reservoirs

Abstract

The present invention relates to the field of reservoir engineering technology, and discloses a physical information neural network numerical simulation method suitable for low permeability water drive reservoirs, the method comprising: using the finite volume method for full implicit solution, considering strong nonlinear seepage; constructing an expert experience loss function, limiting the pressure saturation output range, and improving the solution efficiency; using a fully connected neural network as a basis, adding adjacent position anchoring, adaptive activation function and jump connection gated update structure; using two independent networks to solve pressure and saturation respectively; introducing a well model to deal with large pressure gradient problems; substituting the pressure and saturation calculated by the network into the loss function to perform back propagation optimization of the network. The present invention has high prediction accuracy when dealing with complex problems such as large capillary pressure and nonlinear seepage, and effectively solves the problem that the existing physical information neural network is difficult to apply to the numerical simulation of low permeability water drive reservoirs, and the entire training process does not need to rely on any label data.

Description

Technical Field

[0001] The present invention relates to the technical field of reservoir engineering, and particularly relates to a physical information neural network numerical simulation method applicable to low-permeability waterflood reservoirs. Background Art

[0002] Low-permeability reservoirs have low permeability and insufficient natural energy. To improve the development level of low-permeability oil reservoirs, it is usually necessary to inject water into the formation to supplement energy for waterflood development. In the specific implementation process, to achieve the best development effect, it is necessary to optimize the injection-production plan by combining numerical simulation. Therefore, constructing a numerical simulation method that can accurately simulate the waterflood process of low-permeability oil reservoirs plays an important role in the effective implementation of waterflood development.

[0003] The essence of reservoir numerical simulation is to obtain the fluid state by numerically solving the control equations of underground fluid flow. This method has the problems of long simulation time and high calculation cost. In recent years, with the rapid development of machine learning technology, it has gradually emerged in the field of reservoir engineering. Among them, the physical information neural network numerical simulation method combines the physical knowledge in reservoir engineering with machine learning technology, improving the interpretability and reliability of the neural network on the basis of ensuring fast and efficient, and has good application prospects. However, most of the existing physical information neural network numerical simulation methods use the implicit pressure and explicit saturation method to solve the two-phase flow problem. Although this method has a high calculation speed, it is difficult to effectively handle special complex problems in low-permeability oil reservoirs such as large capillary pressure and nonlinear seepage. Therefore, it is necessary to establish a physical information neural network numerical simulation method applicable to low-permeability waterflood reservoirs. Summary of the Invention

[0004] The purpose of the present invention is to provide a physical information neural network numerical simulation method applicable to low-permeability waterflood reservoirs, aiming to solve the limitations of traditional reservoir numerical simulation methods based on physical information. The core innovations of this method include using the finite volume method for fully implicit solution, considering strong nonlinear seepage; constructing an expert experience loss function to limit the output range of pressure and saturation and improve the solution efficiency; using a fully connected neural network as the basis, adding adjacent position anchoring, adaptive activation functions and skip connection gating update structures; using two independent networks to solve pressure and saturation respectively; introducing a well model to handle large pressure gradient problems; substituting the pressure and saturation calculated by the network into the loss function for backpropagation to optimize the network.

[0005] To achieve the above purpose, the following technical solutions are adopted:

[0006] The present invention provides a physical information neural network numerical simulation method applicable to low-permeability waterflood reservoirs, and the method includes:

[0007] Step 1: Obtain the nonlinear seepage characteristics;

[0008] Step 2: Based on the non-linear seepage characteristics, establish a numerical simulation model for low-permeability waterflooding reservoirs, perform fully implicit numerical discretization using the finite volume method, and construct a physical information loss function applicable to low-permeability waterflooding reservoirs.

[0009] Step 3: According to the injection-production well data, construct an expert experience loss function, and combine the physical information loss function and the expert experience loss function to obtain a total loss function.

[0010] Step 4: Select a fully connected neural network as the basis, and improve it using adjacent position anchoring, an adaptive activation function, and a skip connection gated update structure to construct a network model. The network model includes two independent improved fully connected neural networks, which are respectively used to solve the pressure and saturation.

[0011] Step 5: Introduce a well model to calculate the pressure gradient of the network model.

[0012] Step 6: Input the pressure and saturation of the previous time step into the network model, forward-propagate to calculate the loss value, and use the gradient descent method to update the network parameters in reverse based on the well model and the total loss function to complete the solution of the pressure and saturation for one time step.

[0013] Step 7: Use the pressure, saturation, and training parameters output by the network model at each time step as the input for the next time step, and repeat Step 6 to obtain the pressure and saturation at any time step of the low-permeability waterflooding reservoir.

[0014] Furthermore, the Step 1 includes:

[0015] Establish a discrete grid based on the overall reservoir data; wherein, the overall reservoir data includes reservoir size, shape, and internal geological structure.

[0016] Input the initial reservoir data, rock properties, fluid data, characteristics data between phases, and / or production data into the discrete grid for numerical simulation to obtain non-linear seepage characteristics; wherein, the initial reservoir data includes initial reservoir pressure and / or initial water saturation, the rock properties include permeability distribution, porosity distribution, and / or pore volume compressibility, the fluid data includes oil phase density, oil phase viscosity, oil phase volume compressibility, oil phase viscosity compressibility, water phase density, water phase viscosity, water phase volume compressibility, and / or water phase viscosity compressibility, the characteristics data between phases includes oil-water relative permeability curve and / or oil-water capillary pressure curve, and the production data includes bottom hole pressure of production wells, bottom hole pressure of injection wells, and / or production time.

[0017] Furthermore, the Step 2 includes:

[0018] Let the control equations and auxiliary equations of the same form be satisfied among all phases in the reservoir:

[0019] (1)

[0020] (2)

[0021] (3)

[0022] where, is the density of phase is the seepage velocity of phase is the unit outer normal vector of surface is the mass flow rate of phase is the rock porosity; is the saturation of phase is the capillary force; is the non-wetting phase pressure; is the wetting phase pressure; V is the control volume; A is the control surface area;

[0023] Establish the velocity equation without considering gravity and considering strong non-linear seepage as follows:

[0024] (4)

[0025] where, is the relative permeability of phase is the absolute permeability of the rock; is the viscosity of phase is the pressure gradient of phase, which is expressed as the pressure difference between two adjacent connected grids in the discrete grid;

[0026] and represent the first coefficient and the second coefficient respectively, which are determined by the following equations:

[0027] (5)

[0028] (6)

[0029] where, is the scalar form permeability;

[0030] After discretizing equation (1) using the finite volume method, we get:

[0031] (7)

[0032] wherein, is the volume of the i th grid after discretization; is the i th and the j th grid contact surface unit normal vector pointing from the i th grid to the j th grid; is the i th and the j th grid contact surface area; t is time; is the i set of grids connected to the

[0033] Substitute the velocity equation (4) into formula (7) and take the first-order backward difference with respect to time for the cumulative term to obtain:

[0034] (8)

[0035] wherein, is the size of the time step; represents the current time step, represents the next time step; is the unit volume injection / production rate of the time step; is the phase fluid mass in the grid at the next time step, is the phase fluid mass in the grid at the current time step; is the j grid flow potential; is the i grid flow potential; e is the natural constant; is the product of the phase density and mobility; is the

[0036] is the phase mobility, expressed as:

[0037] (9)

[0038] wherein, is the phase relative permeability;

[0039] is the i th and thej The conductivity between grids, expressed as:

[0040] (10)

[0041] Wherein, is the average permeability between the i th and the j th grids, and are the distances from the center to the connection surface between the i th and the j th grids;

[0042] Construct a residual equation according to formula (8):

[0043] (11)

[0044] Wherein, is the phase residual;

[0045] Establish a loss function based on the residual equation:

[0046] (12)

[0047] Wherein, is the oil phase residual; is the water phase residual; is the oil phase residual of the i th grid; is the water phase residual of the i th grid; MES is the MSE loss function; loss int is the physical equation loss; N is the number of grids.

[0048] Furthermore, step 3 includes:

[0049] Taking the injection well pressure being greater than the reservoir formation pressure and the production well pressure being less than the reservoir formation pressure as limiting conditions, establish a first expert experience loss function, expressed as:

[0050] (13)

[0051] Wherein, is the bottom hole pressure of the injection well, is the bottom hole pressure of the production well, is the number of grids, loss exp,p is the pressure-related expert experience loss, max is the maximum value function, is the oil phase pressure at the next time step;

[0052] Limit the training range of the aqueous phase saturation to and to establish a second expert experience loss function; the second expert experience loss function is expressed as:

[0053] (14)

[0054] where loss exp,S is the saturation-related expert experience loss of saturation, is the aqueous phase saturation at the next time step, represents the irreducible water saturation, which refers to the lowest saturation of the aqueous phase in the reservoir, represents the residual oil saturation, which refers to the lowest saturation of the oil phase in the reservoir;

[0055] Based on the physical information loss function, the first expert experience loss function, and the second expert experience loss function, establish a total loss function expressed as:

[0056] (15)

[0057] where loss is the total loss, loss int is the physical equation loss, loss exp,p is the pressure-related expert experience loss, loss exp,S is the saturation-related expert experience loss.

[0058] Furthermore, the step 4 includes:

[0059] For the input -dimensional initial pressure matrix and the initial saturation matrix perform dimensionality increase in the input layer to obtain the first pressure matrix and the first saturation matrix ;

[0060] Multiply the first pressure matrix and the first saturation matrix by the learnable first parameter respectively and add the adjacency matrix for adjacency matrix anchoring. The adjacency matrix is a matrix of all 0 and 1 elements, with 1 at the connected positions and 0 at the unconnected positions. The calculation process is expressed as:

[0061] (16)

[0062] (17)

[0063] Among them, and are the second pressure matrix and the second saturation matrix respectively, is the adjacency matrix;

[0064] After anchoring and are normalized and input into the hidden layer and then multiplied by the learnable second parameter and the ReLU activation function is used:

[0065] (18)

[0066] (19)

[0067] Where and are the third pressure matrix and the third saturation matrix respectively, is the weight of the hidden layer, is the bias of the hidden layer;

[0068] By and input into the output layer for dimensionality reduction to obtain and , and the data before and after training are multiplied by the weight and added to achieve cross-layer connection:

[0069] (20)

[0070] (21)

[0071] Where is the pressure matrix after cross-layer connection processing, is the saturation matrix after cross-layer connection processing is the learnable third parameter.

[0072] Furthermore, in step 5, the well model is expressed as:

[0073] (22)

[0074] Where is the production; is the bottom-hole pressure; is the pressure of the grid where the well is located;

[0075] is the well index, expressed as:

[0076] (23)

[0077] Where is the oil layer thickness, which is the length of the well in the z - direction of the grid in numerical simulation; is the wellbore radius; is the skin factor; k is the permeability;

[0078] is the reservoir radius, which is expressed in numerical simulation as:

[0079] (24)

[0080] where is the length of the grid in the x - direction; is the length of the grid in the y - direction;

[0081] is the total mobility, which is expressed as:

[0082] (25)

[0083] where, is the oil - phase mobility, is the water - phase mobility.

[0084] Further, in the said step 6, the pressure and saturation of the current time step are respectively input into the network model constructed in step 4 for forward propagation to obtain the pressure and saturation data of the next time step, and the loss value is calculated according to the total loss function of step 3 using the pressure and saturation data of the current time step and the pressure and saturation data of the next time step for backpropagation to update the parameters of the two networks. This process is repeated until the loss value drops to the convergence range, and the pressure and saturation of the current next time step are then solved.

[0085] Further, in the said step 7, the training process of step 6 is repeated, and the pressure and saturation data of the solved next time step are used as the initial pressure and saturation data of the current time step, so as to solve the pressure and saturation data of the next time step again, obtaining the pressure and saturation data of any time step. Based on the pressure and saturation data of any time step, the pressure distribution, saturation distribution, production rate of the production well, injection volume of the injection well, and / or water production rate of the reservoir within the solved time step are solved.

[0086] The beneficial effects of the present invention are:

[0087] The present invention is applicable to low-permeability reservoirs; it considers strong non-linear seepage and uses the fully implicit solution of the finite volume method; constructs an expert experience loss function to limit the output range of pressure saturation to improve the solution efficiency; uses a fully connected neural network as the basis and adds adjacent position anchoring, an adaptive activation function and a skip connection gating update structure, and uses two independent networks to solve pressure and saturation respectively; introduces a well model to handle large pressure gradient problems; substitutes the pressure and saturation calculated by the network into the loss function for backpropagation to optimize the network until the loss is less than the preset value to complete the solution of the current step. Compared with traditional methods, this method not only improves the calculation accuracy, but also can better handle calculation problems under complex conditions such as low-permeability reservoirs, such as large capillary forces and strong non-linearity, on the premise of ensuring the solution efficiency. Description of the Drawings

[0088] Figure 1 It is a flow chart of the physical information neural network numerical simulation method for low-permeability waterflooding reservoirs provided by an embodiment of the present invention.

[0089] Figure 2 It is a capillary force and relative permeability curve graph of the reservoir provided by an embodiment of the present invention.

[0090] Figure 3 It is a reservoir permeability map provided by an embodiment of the present invention.

[0091] Figure 4 It is a specific calculation flow chart of the neural network model provided by an embodiment of the present invention.

[0092] Figure 5 It is a non-linear seepage coefficient graph under different mobilities provided by an embodiment of the present invention.

[0093] Figure 6 It is a true pressure field map, a predicted pressure field map and an error schematic diagram of the two at the last moment provided by an embodiment of the present invention, where (a) is the true pressure field map, (b) is the predicted pressure field map, and (c) is the error schematic diagram.

[0094] Figure 7 It is a true saturation field map, a predicted saturation field map and an error schematic diagram of the two at the last moment provided by an embodiment of the present invention, where (a) is the true saturation field map, (b) is the predicted saturation field map, and (c) is the error schematic diagram.

[0095] Figure 8 It is a comparison graph of the true value and the predicted value of the water production rate of the reservoir provided by an embodiment of the present invention.

[0096] Figure 9 It is a comparison graph of the true value and the predicted value of the cumulative liquid production of the reservoir provided by an embodiment of the present invention.

[0097] Figure 10The production comparison chart considering and not considering non - linear seepage flow provided by the embodiment of the present invention.

[0098] Figure 11 The water - cut comparison chart under different capillary pressures provided by the embodiment of the present invention. Specific embodiments

[0099] The following uses specific specific examples to illustrate the implementation manners of the present invention. Those skilled in the art can easily understand other advantages and effects of the present invention from the content disclosed in this specification. The present invention can also be implemented or applied through other different specific implementation manners. Various details in this specification can also be modified or changed based on different viewpoints and applications without departing from the spirit of the present invention. It should be noted that, without conflict, the following embodiments and the features in the embodiments can be combined with each other.

[0100] The following combines the drawings and embodiments to further describe in detail the specific implementation manners of the present invention.

[0101] The embodiment of the present invention provides a physical - information neural - network numerical simulation method applicable to low - permeability water - flooding reservoirs. This method has high prediction accuracy when dealing with complex problems such as large capillary pressure and non - linear seepage, effectively solves the problem that existing physical - information neural networks are difficult to apply to numerical simulation of low - permeability water - flooding reservoirs, and the entire training process does not require any labeled data. As Figure 1 shown, this method includes steps 1 to 7, which are introduced in detail as follows.

[0102] Step 1: Obtain non - linear seepage characteristics.

[0103] In this embodiment, the non - linear seepage characteristics can be characterized or reflected by rock properties (permeability, porosity), fluid properties (density, viscosity), relative permeability curve, capillary pressure curve, and seepage equation under low - permeability conditions. The non - linear seepage characteristics can be obtained through the following steps:

[0104] Step 1.1: Establish a discrete grid based on the overall reservoir data; wherein, the overall reservoir data includes reservoir size, shape, and internal geological structure.

[0105] Step 1.2: Input the initial reservoir data, rock properties, fluid data, characteristics data between phases and / or production data into the discrete grid for numerical simulation to obtain the non-linear seepage characteristics; wherein, the initial reservoir data includes the initial reservoir pressure and / or initial water saturation, the rock properties include the permeability distribution, porosity distribution and / or pore volume compressibility, the fluid data includes the oil phase density, oil phase viscosity, oil phase volume compressibility, oil phase viscosity compressibility, water phase density, water phase viscosity, water phase volume compressibility and / or water phase viscosity compressibility, the characteristics data between phases includes the relative permeability curve of oil-water phases and / or the capillary pressure curve of oil-water phases, and the production data includes the bottom hole pressure of production wells, the bottom hole pressure of injection wells and / or the production time.

[0106] Step 2: Based on the non-linear seepage characteristics, establish a numerical simulation model for low-permeability waterflooding reservoirs, perform fully implicit numerical discretization using the finite volume method, and construct a physical information loss function applicable to low-permeability waterflooding reservoirs.

[0107] In Step 2, the internal physical information loss function can be constructed through physical equations and implemented by the following steps:

[0108] Step 2.1: Let all phases in the reservoir satisfy control equations and auxiliary equations with the same form:

[0109] (1)

[0110] (2)

[0111] (3)

[0112] where is the density of phase, kg / m 3 ; is the seepage velocity of phase, m / s; is the unit outer normal vector of surface ; is the mass flow rate of phase, kg / s; is the rock porosity, %; is the saturation of phase; is the capillary pressure, Pa; is the non-wetting phase pressure, Pa; is the wetting phase pressure, Pa; V is the control volume; A is the control surface area;

[0113] Step 2.2: Without considering gravity and considering strong non-linear seepage, the velocity equation is as follows:

[0114] (4)

[0115] where is the seepage velocity of phase, m / s; is the relative permeability of phase; is the absolute permeability of the rock, m 2 ; is the viscosity of phase; is the pressure gradient of phase, which is expressed as the pressure difference between two adjacent connected grids in a discrete grid, Pa; and are the first coefficient and the second coefficient.

[0116] where the coefficients , are determined by the following equations:

[0117] (5)

[0118] (6)

[0119] It should be noted that in the above two formulas, is a scalar-form permeability in non-tensor form.

[0120] Step 2.3: After discretizing (1) using the finite volume method, we get:

[0121] (7)

[0122] where is the volume of the i th grid after discretization, m 3 ; is the unit normal vector pointing from the i th grid to the j th grid at the contact surface between the i th grid and the j th grid; is the area of the contact surface between the i th grid and the j th grid, m 2 ; t is time; is the set of grids connected to the i grid;

[0123] Step 2.4: Substitute the velocity equation (4) into (7) and take the first-order backward difference with respect to time for the accumulation term to obtain:

[0124] (8)

[0125] where, is the size of the time step; represents the current time step, represents the next time step; is the injection / production rate per unit volume at the time step; is the mass of the phase fluid in the grid at the next time step, is the mass of the phase fluid in the grid at the current time step; is j the grid flow potential; is i the grid flow potential; e is the natural constant; is the product of the phase density and mobility; is

[0126] where is called the mobility of the

[0127] (9)

[0128] where, is the relative permeability of the

[0129] where is called the conductivity between grids, defined by the following equation:

[0130] (10)

[0131] is the average permeability between grids, 、 are the distances from the grid center to the connection surface, m.

[0132] Step 2.5: According to (8), the residual equation can be constructed:

[0133] (11)

[0134] where, is the residual of the is the size of the time step, s; represents the current time step, represents the next time step, is the injection / production volume per unit volume of the time step, m 3 / (m 3 ·s).

[0135] Since there are only oil and water phases in the simulated reservoir, based on this residual and using the MSE loss function:

[0136] (12)

[0137] where is the oil-phase residual, kg / s; is the water-phase residual, kg / s; is the i oil-phase residual of the is the i water-phase residual of the MES is the MSE loss function; loss int is the physical equation loss; N is the number of grids; in this loss function, the label value is the zero matrix, indicating that the numerical simulation process of this neural network does not require any label data and only needs to satisfy the physical equation.

[0138] Step 3: According to the injection-production well data, construct an expert experience loss function, and combine the physical information loss function and the expert experience loss function to obtain the total loss function.

[0139] Step 3 is implemented through the following steps:

[0140] Step 3.1: Since the injection well pressure must be greater than the reservoir formation pressure and the production well pressure must be less than the reservoir formation pressure for practical production significance, the first expert experience loss function is constructed based on this characteristic:

[0141] (13)

[0142] where, is the bottom-hole pressure of the injection well, is the bottom-hole pressure of the production well, is the number of grids, loss exp,p is the pressure-related expert experience loss, and max is the maximum value function, is the oil-phase pressure at the next time step;

[0143] Step 3.2: For the water-phase saturation, its training range should also be within and in between, where represents the irreducible water saturation, which refers to the lowest saturation of the aqueous phase in the reservoir, represents the residual oil saturation, which refers to the lowest saturation of the oil phase in the reservoir. Based on this characteristic, the second expert experience loss function is constructed as follows:

[0144] (14)

[0145] where loss exp,S is the saturation-related expert experience loss, is the aqueous phase saturation at the next time step.

[0146] The significance of introducing the expert experience loss function is to limit the training range of pressure and saturation and improve the training efficiency.

[0147] Step 3.3: The total loss function is composed of three parts: the physical equation loss function and two expert experience loss functions:

[0148] (15)

[0149] where loss is the total loss, loss int is the physical equation loss, loss exp,p is the pressure-related expert experience loss, loss exp,S is the saturation-related expert experience loss.

[0150] Update the network parameters by backpropagation according to the total loss function until the convergence condition is reached, that is, the physical equation is satisfied and the expert experience is satisfied, and then stop the iteration and output the training result of the current time step.

[0151] Step 4: Based on the fully connected neural network, it is improved by using adjacent position anchoring, adaptive activation function and skip connection gating update structure to construct a network model. The network model includes two independent improved fully connected neural networks, which are used to solve pressure and saturation respectively.

[0152] In Step 4, the fully connected neural network is one of the most basic and common neural networks. It consists of an input layer, hidden layers, and an output layer. The role of the input layer is to receive the input data. The role of the hidden layers is to perform feature extraction and non-linear transformation. There can be multiple hidden layers, and each hidden layer processes the input information. The role of the output layer is to output the final result. The characteristic of the fully connected neural network is that each neuron is connected to all neurons in the previous layer. This fully connected structure enables the network to learn the complex relationship between the input and the output. The hidden layers of the fully connected neural network can be multiple to increase the expressive power of the network. By adjusting the weights of each connection in the network, the relationship between the input and the output can be learned. The fully connected neural network is widely used in various machine learning and deep learning tasks, such as image classification, speech recognition, natural language processing, etc. It is the basis for constructing more complex neural network models.

[0153] Based on this fully connected neural network, Step 4 is implemented through the following steps:

[0154] In Step 4.1, the network forms of the training pressure and saturation are the same. First, the input -dimensional data (the initial pressure matrix and the initial saturation matrix ) are upsampled to -dimensional data in the input layer: the first pressure matrix and the first saturation matrix ; then the upsampled data is multiplied by the learnable first parameter and added with the adjacency matrix for adjacency matrix anchoring. The physical meaning of this step is to improve the training efficiency by inputting the connection information between grids for the neural network. The adjacency matrix is a matrix with all 0 and 1 elements, where the value is 1 at the connected positions and 0 at the unconnected positions:

[0155] (16)

[0156] (17)

[0157] where , are the second pressure matrix and the second saturation matrix respectively, and is the adjacency matrix.

[0158] In Step 4.2, the anchored , are normalized and then input into the hidden layer, multiplied by the learnable second parameter and the ReLU activation function is used. Before the activation function, it is multiplied by a learnable second parameter to make the activation function adaptive:

[0159] (18)

[0160] (19)

[0161] Among them 、 are the third pressure matrix and the third saturation matrix respectively, is the weight of the hidden layer, is the bias of the hidden layer.

[0162] Step 4.3: Input 、 into the input-output layer, and after dimensionality reduction, obtain 、 . At this time, multiply the data before and after training by a certain weight and then add them to achieve cross-layer connection. The purpose of this step is to improve the training efficiency by referring to historical information:

[0163] (20)

[0164] (21)

[0165] Among them is the pressure after cross-layer connection processing, is the saturation after cross-layer connection processing. These two parameters are the finally output data after passing through the neural network. Among them is the third learnable parameter.

[0166] Step 5: Introduce a well model to calculate the pressure gradient of the network model.

[0167] In Step 5, due to the large pressure gradient around the well, it is difficult for training to converge in this case. Therefore, a well model is introduced to address the problem of large pressure gradients and be closer to the actual reservoir development situation. The definition of the well model is as follows:

[0168] (22)

[0169] Among them is the production rate, m 3 / s; is the bottom-hole pressure, Pa; is the pressure of the grid where the well is located, Pa;

[0170] Among them is the well index, defined as:

[0171] (23)

[0172] Among them is the oil layer thickness, in m, which is the length of the grid where the well is located in the z direction in numerical simulation; is the wellbore radius, in m; is the skin factor, k is the permeability.

[0173] Among them is the oil reservoir radius, which is in numerical simulation:

[0174] (24)

[0175] Among them is the length of the grid in the x direction, in m; is the length of the grid in the y direction, in m.

[0176] is the total mobility, 1 / (Pa·s), defined as:

[0177] (25)

[0178] Among them, is the oil-phase mobility, is the water-phase mobility.

[0179] The significance of introducing the well model is that directly using the source-sink term to calculate the grid pressure will lead to a very large pressure gradient in the grid where the well is located, and it is difficult for both traditional numerical simulation and neural networks to converge. After introducing the well model, the pressure gradient is reduced and it is closer to the actual situation of the oil reservoir.

[0180] Step 6: Input the pressure and saturation at the previous time step into the network model, calculate the loss value through forward propagation, and use the gradient descent method to update the network parameters backward based on the well model and the total loss function to complete the solution of the pressure and saturation at one time step.

[0181] In Step 6, the pressure and saturation at the current time step are respectively propagated forward in the two networks described in Step 4 to obtain the pressure and saturation data at the next time step, and the loss value is calculated according to the loss function described in Step 2 using the pressure and saturation data at the current time step and the pressure and saturation data at the next time step for backward propagation to update the parameters of the two networks. Repeat this process until the loss value drops to the convergence range, and the pressure and saturation at the current next time step are then solved.

[0182] Step 7: Use the pressure, saturation, and training parameters output by the network model at each time step as the input for the next time step, repeat Step 6, and obtain the pressure and saturation at any time step of the low-permeability waterflood oil reservoir.

[0183] In Step 7, repeat the training process of Step 6. Take the pressure and saturation data solved for the next time step as the initial pressure and saturation data for the current time step, and then the pressure and saturation data for the next time step can be solved. By repeating Step 6, the pressure and saturation data for any time step can be solved. With the pressure and saturation data for each time step, a series of data such as the pressure distribution, saturation distribution, production rate of production wells, injection volume of injection wells, and water production rate within the solved time steps of the reservoir can be obtained.

[0184] A conceptual example will be provided below to verify the feasibility and advancement of the method proposed in the present invention.

[0185] The reservoir model uses a discretized grid of 25×25×1, with a grid size of 4m×4m×10m. The initial reservoir pressure is 25MPa, the initial saturation is 0.2, the initial porosity is 0.05, and the permeability field is as Figure 3 shown. The rock pore compressibility is 1×10 -9 1 / Pa. The model includes a production well and an injection well, with constant bottom-hole flowing pressure for production and injection. The bottom-hole flowing pressure of the production well is 5MPa, and the bottom-hole flowing pressure of the injection well is 35MPa. The volume compressibility of the fluid is 1×10 -8 1 / Pa, the initial density of the oil is 850kg / m 3 , the initial viscosity is 2mPa·s, the initial density of the water is 1000kg / m 3 , the initial viscosity is 1mPa·s, the capillary pressure curve and relative permeability curve are as Figure 2 shown, and the seepage velocity equation is as shown in Equation (4). Figure 5 is the curve of the , coefficients in the non-linear seepage velocity equation varying with . The simulation time is 1000 days, and the time step is 0.1 day.

[0186] Construct a fully connected neural network model, as Figure 4 shown. The fully connected neural network includes an input layer, an output layer, and a hidden layer. The pressure and saturation are solved through two such networks respectively.

[0187] Figure 6 For the pressure distribution after 1000 days of simulation, the pressure distribution calculated by the traditional numerical simulator, and the comparison between the two, the maximum error between the two is only 0.16MPa. Figure 7 For the saturation distribution after 1000 days of simulation, the saturation distribution calculated by the traditional numerical simulator, and the comparison between the two, the maximum error between the two is only 0.008, indicating that the result of the neural network numerical simulation has a very high accuracy.

[0188] Figure 8For the comparison of the water cut of production wells between neural network numerical simulation and traditional numerical simulation, Figure 9 For the comparison of the cumulative liquid production between neural network numerical simulation and traditional numerical simulation, it can be obtained by comparison that the numerical simulation through the neural network has the same accuracy as the traditional numerical simulator.

[0189] Figure 10 For the comparison of the cumulative liquid production with and without considering non - linear seepage, it can be obtained by comparison that considering non - linear seepage will significantly reduce the cumulative liquid production, and large capillary force will also significantly reduce the cumulative liquid production.

[0190] Figure 11 For the comparison of the water cut of production wells under different capillary forces, it can be obtained by comparison that in the case of large capillary force, water breakthrough occurs earlier but the water cut curve rises more slowly.

[0191] The above embodiments are only used to illustrate the present invention and are not intended to limit the present invention. Those of ordinary skill in the relevant technical fields can also make various changes and modifications without departing from the spirit and scope of the present invention. Therefore, all equivalent technical solutions also belong to the scope of the present invention, and the patent protection scope of the present invention shall be defined by the claims.

Claims

1. A physical information neural network numerical simulation method suitable for low permeability water drive oil reservoirs, characterized in that: The method comprises: Step 1: Obtain nonlinear seepage characteristics; Step 2: A numerical simulation model of a low permeability water drive reservoir is established based on the nonlinear seepage characteristics, a finite volume method is used for full implicit numerical discretization, and a physical information loss function suitable for a low permeability water drive reservoir is constructed; Step 3: construct an expert experience loss function based on the injection and production well data, and combine the physical information loss function and the expert experience loss function to obtain a total loss function; Step 4: A fully connected neural network is selected as the basis, and an adjacent position anchoring, an adaptive activation function, and a skip connection gated update structure are used for improvement to construct a network model, wherein the network model includes two independent improved fully connected neural networks, which are used to solve pressure and saturation respectively; Step 5: Introduce the well model to calculate the pressure gradient of the network model; Step 6: Input the pressure and saturation of the previous time step into the network model, forward propagate and calculate the loss value, and use the gradient descent method to reversely update the network parameters based on the well model and the total loss function to complete the pressure and saturation solution for one time step; Step 7: Use the pressure, saturation and training parameters output by the network model at each time step as the input for the next time step, repeat step 6, and obtain the pressure and saturation of the low permeability water drive reservoir at any time step.

2. The method according to claim 1, characterized in that The step 1 comprises: Establishing a discrete grid based on the overall reservoir data; wherein the overall reservoir data includes the reservoir size, shape and internal geological structure; Inputting initial reservoir data, rock properties, fluid data, characteristic data between phases and / or production data into the discrete grid for numerical simulation to obtain nonlinear seepage characteristics; wherein the initial reservoir data includes initial reservoir pressure and / or initial water saturation, the rock properties include permeability distribution, porosity distribution and / or pore volume compression coefficient, the fluid data includes oil phase density, oil phase viscosity, oil phase volume compression coefficient, oil phase viscosity compression coefficient, water phase density, water phase viscosity, water phase volume compression coefficient and / or water phase viscosity compression coefficient, the characteristic data between phases include oil-water phase relative permeability curve and / or oil-water phase capillary force curve, and the production data includes production well bottom hole pressure, injection well bottom hole pressure and / or production time.

3. The method according to claim 1, characterized in that The step 2 comprises: Assume that all phases in the reservoir satisfy the same governing equations and auxiliary equations: (1) (2) (3) in, for Density of phase; for The percolation velocity of the phase; For face The unit external normal vector of ; for Phase mass flow rate; is the rock porosity; for Saturation of phase; is the capillary force; is the non-wetting phase pressure; is the wetting phase pressure; V is the volume of the control body; A To control body surface area; The velocity equation that ignores gravity and considers strong nonlinear seepage is established as follows: (4) in, for relative permeability of phases; is the absolute permeability of rock; for Phase viscosity; for Phase pressure gradient, expressed as the pressure difference between two adjacent connected grids in a discrete grid; and denote the first coefficient and the second coefficient respectively, which are determined by the following equations: (5) (6) in, is the permeability in scalar form; After discretizing equation (1) using the finite volume method, we get: (7) in, After discrete i The volume of the grid; For the i and j The mesh contact surface is composed of i The grid points to j The unit normal vector of the mesh; For the i and j The area of ​​the mesh contact surface; t For time; For i Grid A collection of connected grids; Substituting the velocity equation (4) into equation (7) and taking the first-order backward difference of the accumulation term with respect to time, we obtain: (8) in, is the size of the time step; represents the current time step, represents the next time step; for The amount of injection / production per unit volume in the time step; For the next time step, Phase fluid mass, The current time step grid Phase fluid quality; for j Grid flow potential; for i Grid flow potential; e is a natural constant; for The product of phase density and fluidity; for Phase pressure gradient; for The fluidity of the phase is expressed as: (9) in, for relative permeability of phases; For the i and j The conductivity between the grids is expressed as: (10) in, For the i and j The average permeability among the grids is and For the i and j The distance from the center of each grid to the connecting surface; According to formula (8), the residual equation is constructed: (11) in, for The residual of the phase; The loss function is established based on the residual equation: (12) in, is the oil phase residual; is the water phase residual; For the i The oil phase residual of the grid; For the i The water phase residual of each grid; MES is the MSE loss function; loss int is the loss of physical equation; N is the number of grids.

4. The method according to claim 1, characterized in that The step 3 comprises: Taking the injection well pressure greater than the reservoir formation pressure and the production well pressure less than the reservoir formation pressure as the restriction conditions, the first expert experience loss function is established, which is expressed as: (13) in, is the bottom hole pressure of the injection well, is the bottom hole pressure of the production well, is the number of grids, loss exp,p is the pressure-related expert experience loss, max is the maximum value function, is the oil phase pressure at the next time step; Limit the training range of water phase saturation to and A second expert experience loss function is established; the second expert experience loss function is expressed as: (14) in, loss exp,S is the saturation-related saturation expert experience loss, is the water phase saturation at the next time step, Represents irreducible water saturation, which refers to the lowest saturation of water phase in the reservoir. Residual oil saturation refers to the lowest saturation of oil phase in the reservoir; Based on the physical information loss function, the first expert experience loss function and the second expert experience loss function, the total loss function is established as: (15) in, loss is the total loss, loss int is the physical equation loss, loss exp,p For stress-related expert experience loss, loss exp,S is the loss of saturation-related expert experience.

5. The method according to claim 4, characterized in that The step 4 comprises: For input The initial pressure matrix and the initial saturation matrix Perform dimension increase at the input layer to obtain the first pressure matrix and the first saturation matrix ; The first pressure matrix and the first saturation matrix Multiply by the first learnable parameter And add the adjacency matrix to anchor the adjacency matrix. The adjacency matrix is ​​a matrix with all 0 and 1 elements. The connected position is 1 and the unconnected position is 0. The calculation process is expressed as: (16) (17) in, , are the second pressure matrix and the second saturation matrix respectively. is the adjacency matrix; After anchoring , Normalize and input into the hidden layer and multiply by the second learnable parameter And use the ReLU activation function: (18) (19) in , The third pressure matrix and the third saturation matrix are respectively, is the weight of the hidden layer, is the bias of the hidden layer; Will , After dimensionality reduction in the input and output layers, we get , , and multiply the data before and after training by the weights and then add them together to achieve cross-layer connection: (20) (21) in is the pressure matrix after cross-layer connection processing, is the saturation matrix after cross-layer connection processing is the third learnable parameter.

6. The method according to claim 1, characterized in that In step 5, the well model is expressed as: (22) in, for output; is the bottom hole pressure; is the pressure of the grid where the well is located; is the well index, expressed as: (23) in, is the thickness of the oil layer, which is the length of the grid in the z direction where the well is located in the numerical simulation; is the wellbore radius; is the skin factor; k is the permeability; is the reservoir radius, which is expressed in numerical simulation as: (24) in is the length of the grid in the x direction; is the length of the grid in the y direction; is the total flow rate, expressed as: (25) in, is the oil phase fluidity, is the water phase mobility.

7. The method according to claim 1, characterized in that In step 6, the pressure and saturation of the current time step are respectively input into the network model constructed in step 4 for forward propagation to obtain the pressure and saturation data of the next time step, and the pressure and saturation data of the current time step and the pressure and saturation data of the next time step are used to calculate the loss value according to the total loss function of step 3 for back propagation to update the parameters of the two networks, and this process is repeated until the loss value drops to the convergence range, and the pressure and saturation of the current next time step are solved.

8. The method according to claim 1, characterized in that In step 7, the training process of step 6 is repeated, and the solved pressure and saturation data of the next time step are used as the initial pressure and saturation data of the current time step, so as to solve the pressure and saturation data of the next time step again, and obtain the pressure and saturation data of any time step. Based on the pressure and saturation data of any time step, the pressure distribution, saturation distribution, production of production wells, injection volume of injection wells and / or water production rate of the oil reservoir in the solution time step are solved.

Images