Heterogeneous reservoir seepage simulation method based on PIKAN

Abstract

The invention provides a heterogeneous reservoir seepage simulation method based on PIKAN, and belongs to the technical field of heterogeneous reservoir development, and the method comprises the following steps: S1, using three independent KANs to respectively approximate PIKAN of pressure p, velocity component vx and velocity component vy; and S2, constructing a loss function of a mixed pressure-speed formula including the residual error of the mass conservation equation, the residual error of the Darcy law and the residual error of the related boundary condition, and optimizing the loss function. According to the heterogeneous reservoir seepage simulation method based on the PIKAN, KAN is introduced into porous medium flow simulation for the first time, and a preliminary reference is provided for developing a general seepage simulation tool based on the PIKAN.

Description

Technical Field

[0001] The invention relates to the technical field of heterogeneous oil reservoir development, and in particular to a heterogeneous reservoir seepage simulation method based on PIKAN. Background Art

[0002] Multilayer Perceptron (MLP) is a deep learning model based on a feedforward neural network. The multi-layer nested internal structure of MLP and the interaction between weights and nodes enhance its expressiveness and ability to adapt to complex tasks, but this also makes the model lack of interpretability and may over-fit to the noise and details of the training data, resulting in overfitting problems.

[0003] PINN (Physics-Informed Neural Network) based on MLP is a semi-supervised or unsupervised learning neural network method. It can not only learn labeled training data samples like traditional neural networks, but also directly obtain approximate solutions to differential equations based on the information of the differential equation problem without labeled data. PINN has been proven to be effective and interpretable in solving complex real-world problems with noisy data and involving physical laws, making it widely used in various scientific and engineering fields.

[0004] At present, PINN has also been widely used in flow simulation in porous media. In the current study, the fluid pressure is substituted into the continuity equation as a state variable through Darcy's law to simulate the flow in porous media. The continuity equation, boundary conditions and initial conditions are embedded in the loss function as physical constraints to ensure that the neural network learns a solution that conforms to the law of seepage. The gradient of the loss function to the neural network parameters relies on the automatic differentiation technique to calculate the spatial derivative of the permeability of the medium. However, porous media are usually heterogeneous. When the permeability is discontinuous in space, the use of automatic differentiation techniques will lead to non-physical and non-conservative solutions. To address this problem, Tartakovsky et al. in 2020 and He and Tartakovsky in 2021 used analytical permeability expressions to represent heterogeneity, but it may be very difficult to obtain such analytical expressions in practical applications; Zhang et al. in 2023 and Yan et al. in 2024 used the finite volume method to discretize the seepage equation in the heterogeneous computational domain, thereby replacing the original continuous partial differential equation with a discrete finite volume discretization format as the physical part of the error function. However, this implementation method may cause PINNs to lose the advantage of directly using differential equations to obtain approximate solutions, and its accuracy will also be limited by the error of the numerical discretization format itself, reducing the necessity of its application, because when the numerical discretization format is available, directly calculating the numerical solution is generally much more convenient than iterative training of neural networks that consumes a lot of computing resources. In addition, the domain decomposition technology can also be used to divide the entire heterogeneous computational domain, so that each subdomain can be processed separately. However, this also depends on the actual heterogeneous distribution. In the case of strong heterogeneity, it is still difficult to directly apply or the robustness of its method is poor. In 2023, Lehmann et al. proposed using a mixed pressure-velocity formula to replace the pure pressure model that only uses pressure as the main unknown for PINN modeling of heterogeneous porous media flow, and found that the mixed formula can achieve higher calculation accuracy and better convergence. This method can directly avoid the derivation of discontinuous heterogeneous functions and is a robust method for dealing with strongly heterogeneous reservoir properties.

[0005] The successful application of PINN is due to its inclusion of MLP as the main neural network structure to approximate complex nonlinear mapping relationships. The universal approximation theorem shows that when MLP has enough neurons and layers, it can approximate any continuous function with arbitrary accuracy. Although MLP is widely used as the cornerstone of today's deep learning models, it also has significant shortcomings. In 1987, Hecht-Nielsen proposed an alternative to MLP - the Kolmogorov network, which is built on the Kolmogorov-Arnold representation theorem, which asserts that any sufficiently smooth multivariable function can be represented as a superposition of univariate functions. However, the properties of univariate functions in the original Kolmogorov-Arnold representation theorem are usually poor, which makes the Kolmogorov network with shallow structure almost ineffective in practical applications.

[0006] Recently, in 2024, Liu et al. proposed KAN (Kolmogrov-Arnold Network), which deepens the shallow Kolmogorov network to construct univariate functions with good properties. The core feature of KAN is to have learnable activation functions on the edges (i.e. weights) of the network, instead of using fixed activation functions on the nodes (i.e. neurons) like traditional MLPs. Therefore, KANs do not have a linear weight matrix. Instead, each weight parameter is replaced by a learnable one-dimensional function, parameterized as a spline function, and the nodes of KANs simply sum the incoming signals without applying any nonlinearity. Compared with MLPs, KANs achieve higher accuracy with fewer parameters in various tasks, bring greater performance improvements as the model size increases, and provide interpretability and interactivity that MLPs do not have. Although each MLP weight parameter becomes a spline function of KAN, KAN generally allows for smaller computational graphs than MLP, thus not introducing additional computational costs. In the original KAN, B-splines are used for activation function construction due to their excellent fitting ability. Numerical analysis includes a variety of orthogonal polynomials (Hermite polynomials, Legendre polynomials, and Chebyshev polynomials), including other traditional approximation functions (such as radial basis functions), which can replace the B-spline function in the original KAN. Therefore, the essential core of KAN is a composite function framework composed of learnable activation functions.

[0007] KAN is currently also used to solve differential equations and operator learning. In 2024, Wang et al. systematically compared the computational performance of KAN-based PIKNN and MLP-based PINNs in different forms of PDEs, proving the superiority of PIKAN in dealing with many mechanical problems.

[0008] However, PIKAN has not been fully applied to the field of porous media seepage. Therefore, the present invention intends to construct for the first time a physical information KAN (PIKAN, i.e., replacing the MLP-based neural network in PINN with KAN) for seepage simulation in heterogeneous porous media, and explore the computational performance of PIKAN and the advantages of the mixed pressure-velocity formula over the pure pressure formula, so as to promote PIKAN to receive more and more attention and application in the field of numerical simulation of porous media seepage. Summary of the invention

[0009] The purpose of the present invention is to provide a PIKAN-based heterogeneous reservoir seepage simulation method, to introduce the promising KAN into porous media flow simulation for the first time, and to provide a preliminary reference for the development of a general seepage simulation tool based on PIKAN.

[0010] To achieve the above object, the present invention provides a PIKAN-based heterogeneous reservoir seepage simulation method, comprising the following steps:

[0011] Step S1: Use three independent KAN networks to approximate the pressure p and velocity component v respectively x and the velocity component v y PIKAN;

[0012] Step S2: construct a loss function of a mixed pressure-velocity formula including the residual of the mass conservation equation, the residual of Darcy's law and the residual of related boundary conditions, and optimize the loss function.

[0013] Preferably, in step S1, three independent KAN networks are used to approximate the pressure p, the velocity component v x and the velocity component v y The specific operations of PIKAN are as follows:

[0014] For a smooth have:

[0015]

[0016] in, Φ q :→; x represents n independent variables x i The vector is composed of i=1,L,n; f(x) represents the n independent variables x i A function of x p represents the pth component of vector x; univariate function φ q,p Process the pth component of vector x and contribute a term to the sum of the qth external function; there are 2n+1 external functions, each of which has q is a univariate function;

[0017] The Chebyshev polynomial function is used to replace the B-spline function;

[0018] For a fully connected KAN with L+1 layers, the input layer is denoted as layer 0, n l represents the number of neurons in the lth layer, (l,j) represents the jth neuron in the lth layer, represents the activation value of the neuron; between layer l and layer l+1, there will be n l n l+1 parameterized univariate functions, where the univariate function connecting (l,j) and (l+1,i) is denoted by Then the activation value of neuron (l+1,i) is calculated as:

[0019]

[0020] in, represents the value of the jth neuron in the l+1th layer; represents the value of the i-th neuron in the l-th layer; n l represents the number of neurons in the lth layer;

[0021] Therefore, according to the activation value x of the neuron in the l+1th layer l Calculate the vector x consisting of the activation values ​​of the neurons in the l+1th layer l +1 The matrix form expression of is:

[0022] x l+1 =T l x l (3);

[0023] Among them, T l represents the function matrix corresponding to the lth KAN layer, T l The i-th row and j-th column of is a Chebyshev polynomial function;

[0024] All Chebyshev polynomial functions in KAN All learnable parameters are denoted as θ, l = 0, 1, L, L-1, j = 1, 2, L, n l ,i=1,2,L,n l+1 , we get: When the input layer data is x 0 When , the output layer neuron value of the fully connected KAN with L+1 layers is calculated as:

[0025]

[0026] in, Represents the composition of functions.

[0027] Preferably, in step S2, the loss function of the mixed pressure-velocity formula is:

[0028]

[0029] Among them, L(θ) represents the loss function; ω1, ω2, ω3, ω4 all represent weighting coefficients; represents the residual of the mass conservation equation; represents the residual of Darcy's law; L DB represents the residual of the Dirichlet boundary condition; L NB represents the residual of the Neumann boundary condition.

[0030] Preferably, the residual of the mass conservation equation is The expression is as follows:

[0031]

[0032] Among them, N in Indicates the number of points inside the computational domain; (x i ,y i ) represents the coordinates of a certain point; x and y represent the spatial independent variables of the function; v x 、v y They represent the components of velocity v in the x and y directions respectively.

[0033] Preferably, the residual of Darcy's law The expression is as follows:

[0034]

[0035] Here, k represents the permeability.

[0036] Preferably, the expression of the residual of the Dirichlet boundary condition is as follows:

[0037]

[0038] Among them, N DB represents the number of collocation points on the Dirichlet frontier; p represents the number of collocation points (x i ,y i ) is the PIKAN prediction value on ; g represents the Dirichlet boundary condition value.

[0039] Preferably, the residual of the Neumann boundary condition is expressed as follows:

[0040]

[0041] Among them, N NB represents the number of collocation points on the Neumann boundary; h xrepresents the x-direction velocity value given by the Neumann boundary condition; h y Represents the y-direction velocity value given by the Neumann boundary condition.

[0042] Therefore, the present invention adopts the above-mentioned PIKAN-based heterogeneous reservoir seepage simulation method, introduces the promising KAN into porous media flow simulation for the first time, and provides a preliminary reference for the development of a general seepage simulation tool based on PIKAN. BRIEF DESCRIPTION OF THE DRAWINGS

[0043] Figure 1 Comparison between the MLP on which PINN is based and the KAN on which PIKAN is based;

[0044] Figure 2 It is a schematic diagram of the PV-3-PIKAN structure;

[0045] Figure 3 is the reservoir domain in Example 1; wherein, Figure 3 (a) is the computational domain; Figure 3 (b) in the figure is the matching point;

[0046] Figure 4 are the loss functions and their sub-items of different methods in Example 1; wherein, Figure 4 (a) is the total loss; Figure 4 (b) in the figure is the PDE loss; Figure 4 (c) in the equation is the Dirichlet boundary condition loss; Figure 4 (d) in the figure is the Neumann boundary condition loss;

[0047] Figure 5 is the pressure distribution predicted by different methods in Example 1; wherein, Figure 5 (a) is the exact solution; Figure 5 (b) in the figure is P-PIKAN; Figure 5 (c) in the figure is PV-3-PINN; Figure 5 (d) in the figure is PV-3-PIKAN;

[0048] Figure 6 The v predicted by different methods in Example 1 is x Distribution; among them, Figure 6 (a) is the exact solution; Figure 6 (b) in the figure is PV-3-PINN; Figure 6 (c) in the figure is PV-3-PIKAN;

[0049] Figure 7 is the v predicted by different methods in Example 1 y Distribution; among them, Figure 7 (a) is the exact solution; Figure 7 (b) in the figure is PV-3-PINN; Figure 7 (c) in the figure is PV-3-PIKAN;

[0050] Figure 8 is the reservoir domain in Example 2;

[0051] Fig. 9 are the loss functions and their sub-items of different methods in Example 2; wherein, Fig. 9 (a) is the total loss; Fig. 9 (b) in the figure is the PDE loss; Fig. 9 (c) in the equation is the Dirichlet boundary condition loss; Fig. 9 (d) in the figure is the Neumann boundary condition loss;

[0052] Fig.10 is the pressure distribution predicted by different methods in Example 2; wherein, Fig.10 (a) is the analytical solution; Fig.10 (b) in the figure is P-PIKAN; Fig.10 (c) in the figure is PV-3-PINN; Fig.10 (d) in the figure is PV-3-PIKAN;

[0053] Fig.11 The v predicted by different methods in Example 2 is x Distribution; among them, Fig.11 (a) is the exact solution; Fig.11 (b) in the figure is PV-3-PINN; Fig.11 (c) in the figure is PV-3-PIKAN;

[0054] Fig.12 is the v predicted by different methods in Example 2 y Distribution; among them, Fig.12 (a) is the exact solution; Fig.12 (b) in the figure is PV-3-PINN; Fig.12 (c) in the figure is PV-3-PIKAN;

[0055] Fig.13 Reservoir domains and collocation points generated for implementing three Latin hypercube sampling strategies; among them, Fig.13 (a) in the figure is the reservoir domain; Fig.13 (b) in the figure is the points used for training;

[0056] Fig.14 are the loss functions and their sub-items of different methods in Example 3; wherein, Fig.14 (a) is the total loss; Fig.14 (b) in the figure is the PDE loss; Fig.14(c) in the equation is the Dirichlet boundary condition loss; Fig.14 (d) in the figure is the Neumann boundary condition loss;

[0057] Fig.15 is the pressure distribution predicted by different methods in Example 3; wherein, Fig.15 (a) is the exact solution; Fig.15 (b) in the figure is P-PIKAN; Fig.15 (c) in the figure is PV-3-PINN; Fig.15 (d) in the figure is PV-3-PIKAN;

[0058] Fig.16 The v predicted by different methods in Example 3 is x Distribution; among them, Fig.16 (a) is the exact solution; Fig.16 (b) in the figure is P-PIKAN; Fig.16 (c) in the figure is PV-3-PINN; Fig.16 (d) in the figure is PV-3-PIKAN;

[0059] Fig.17 is the v predicted by different methods in Example 3 y Distribution; among them, Fig.17 (a) is the exact solution; Fig.17 (b) in the figure is P-PIKAN; Fig.17 (c) in the figure is PV-3-PINN; Fig.17 (d) in the figure is PV-3-PIKAN. DETAILED DESCRIPTION

[0060] The technical solution of the present invention is further described below through the accompanying drawings and embodiments.

[0061] Unless otherwise defined, technical or scientific terms used in the present invention shall have the common meanings understood by one having ordinary skills in the field to which the present invention belongs.

[0062] The theory used in the present invention is first introduced below.

[0063] Governing equations.

[0064] The governing equations for seepage in heterogeneous porous media mainly include the fluid mass conservation equation and Darcy's law obtained based on physical experiments. Assuming steady-state flow and no source or sink terms, the mass conservation equation for the fluid is:

[0065]

[0066] in, represents the divergence operator; v represents the seepage velocity of the fluid. In the two-dimensional case, v=(v x ,v y ) T ;v x 、v y They represent the components of velocity v in the x and y directions respectively.

[0067] Darcy's law is described as:

[0068]

[0069] Where k(x,y) represents the heterogeneous permeability; μ represents the fluid viscosity; p represents the fluid pressure;

[0070] Substituting equation (2) into equation (1), we can obtain an equation only about pressure p:

[0071]

[0072] There are two formulas for simulating the seepage flow. One is a mixed pressure-velocity formula consisting of formula (1) and formula (2) that takes v and p as the main unknowns to solve. The other is formula (3) that only takes p as the main unknown to solve, which is called the pure pressure formula. In the field of reservoir numerical simulation, the pure pressure formula is the most common.

[0073] In the absence of internal source and sink terms, as shown in equation (4), the seepage problem generally has both Dirichlet boundary conditions representing constant pressure boundaries and Neumann boundary conditions representing closed or constant flow.

[0074] Dirichlet boundary condition:

[0075] Neumann boundary condition:

[0076] Among them, Γ DB represents the Dirichlet boundary condition; Γ NB represents the Neumann boundary; h(x,y)=(h x ,h y ); p(x,y) represents the pressure function to be calculated; g(x,y) represents the pressure function given by the Dirichlet boundary condition at the boundary Γ DB The value on the boundary; h(x,y) represents the velocity function given by the Neumann boundary at the boundary Γ NB The value of h x 、h y They represent the components of h(x,y) in the x and y directions respectively.

[0077] Physical Information Neural Network (PINN).

[0078] The core idea of ​​PINN is to use a neural network based on MLP to approximate the solution of partial differential equations. Taking the pure pressure formula in equation (3) as an example, inspired by the universal approximation theorem of MLP, PINN is constructed by Figure 1 The neural network shown has spatial coordinates x and y as input and pressure p as output to approximate the solution p(x,y) to the well-defined problem posed by the boundary conditions in equations (3) and (4).

[0079] Assume that there are L+1 fully connected MLP layers in PINN, that is, there are L-1 hidden layers, the input layer is denoted as layer 0, and n l represents the number of neurons in the lth layer, (l,j) represents the jth neuron in the lth layer, represents the value of the neuron, x l represents the vector of values ​​of all neurons in layer l. Then, we can get:

[0080]

[0081] Among them, the output x of the last layer L Used to approximate the true solution. l-1 and b l Represent the weight matrix and bias vector of the lth layer respectively; x L represents the value of the neuron on the Lth layer; W L-1 represents the weight matrix between the L-1th layer neurons and the Lth layer neurons; b l represents the bias vector of the lth layer; x L-1 represents the value of the neuron on the L-1 layer; b L represents the bias vector of the Lth layer; σ(·) represents a nonlinear activation function. All trainable model parameters, such as weights and biases, are denoted as θ in the following text.

[0082] The above network structure is extended to the general MLP representation form:

[0083]

[0084] In the formula, Represents the composition between functions; N l x=W l-1 x+b l .

[0085] By arranging points inside and on the boundary of the computational domain, assuming that the number of points inside the computational domain, on the Dirichlet boundary, and on the Neumann boundary are N in 、N DB 、N NB, then PINN trains the parameters of the neural network by making the points of the neural network inside the computational domain satisfy the differential equation as much as possible, and making the points of the neural network on the boundary of the computational domain satisfy the boundary conditions as much as possible. Therefore, during the training process, the loss function L(θ) of the neural network is defined as:

[0086] L(θ)=ω1L PDE +ω2L DB +ω3L NB (7);

[0087]

[0088] Among them, ω1, ω2, ω3 represent the weighting coefficients of different loss terms; L PDE represents the residual term of the control equation, L DB and L NB They represent the residual terms of the Dirichlet boundary condition and the Neumann boundary condition, respectively, as shown in Equation (8), Equation (9) and Equation (10). Generally speaking, the mean square error of the L2 norm of the corresponding sample points is used to calculate the three types of losses in Equation (7).

[0089] Therefore, PINN converts the well-defined problem formed by equations (3) and (4) into a problem that is solved by training the neural network and iteratively updating the neural network parameters θ to minimize the loss function L(θ). Usually, the ADAM (Adaptive Moment Estimation) optimizer, a gradient-based first-order optimization adaptive algorithm, is used to optimize the model parameters θ.

[0090] One of the keys to calculating equations (8) and (10) is to calculate partial derivatives through automatic differentiation (AD). AD tracks the influence of the output of each layer of neurons on the output of the previous layer through the chain rule, and directly calculates the derivative of the network output relative to the input in the computational graph. Partial derivatives can be calculated using explicit expressions to obtain exact values, thus avoiding the introduction of truncation errors in traditional numerical approximations. At present, AD is widely used in machine learning and optimization algorithms, which facilitates the development of PINN.

[0091] Introduction to the basic theory of KAN.

[0092] The idea of ​​the KAN network structure comes from the Klmogorov-Arnold representation theorem, which was proposed by mathematician Kolmogorov. The theorem states that if f is any multivariable continuous function defined on a bounded domain, then the function f can be expressed as a two-layer nested addition of a finite number of single-variable, continuous functions, that is, any multivariable continuous function can be expressed as a combination of single-variable continuous functions and addition operations. That is, for a smooth have:

[0093]

[0094] in, Φ q :→,x p Represented as the pth component of the vector x, the univariate function φ q,p Process the pth component of vector x and contribute a term to the sum of the qth external function. There are 2n+1 external functions, each with q is a unary function.

[0095] When the Kolmogorov network is constructed strictly following the Kolmogorov-Arnold representation theorem, the corresponding φ q,p and Φ q The properties of are usually very poor, which makes the Kolmogorov network with shallow structure almost ineffective in practical applications. KAN deepens the shallow Kolmogorov network to construct a single variable function with good properties. In the original KAN, KAN parameterizes each one-dimensional function as a B-spline function with learnable coefficients of local B-spline basis functions. In the present invention, Chebyshev polynomial functions are used instead of B-spline functions.

[0096] Similarly, suppose we have a fully connected KAN with L+1 layers, with the input layer being layer 0 and n l represents the number of neurons in the lth layer, (l,j) represents the jth neuron in the lth layer, Represents the activation value of the neuron. Between layer l and layer l+1, there will be n l n l+1 parameterized univariate functions, where the univariate function connecting (l,j) and (l+1,i) is denoted by Then the activation value of neuron (l+1,i) is calculated as:

[0097]

[0098] Therefore, according to the activation value x of the neuron in the l+1th layer l Calculate the vector x consisting of the activation values ​​of the neurons in the l+1th layerl +1 The matrix form expression of is:

[0099]

[0100] Where, T l is the function matrix corresponding to the lth KAN layer, with T l The i-th row and j-th column of In the present invention, is a Chebyshev polynomial function.

[0101] All Chebyshev polynomial functions in KAN are still (l=0,1,L,L-1,j=1,2,L,n l ,i=1,2,L,n l+1 ) are denoted as θ, we can get: when the input layer data is x 0 When , the output layer neuron value of the fully connected KAN with L+1 layers is calculated as:

[0102]

[0103] The Chebyshev polynomial function has the same good continuity as the B-spline function, so similar to PINN, the derivative of the pressure or velocity of the output layer with respect to the spatial coordinates of the input layer can be obtained through automatic differentiation, and error back propagation can also be used to train KAN.

[0104] PIKAN constructed in the present invention.

[0105] like Figure 2 As shown in the figure, PIKAN essentially replaces the MLP based on PINN with KAN. Therefore, the essential difference between PIKAN and PINN is that the basic neural network architecture is different, but the loss function, optimization algorithm, etc. are the same. In 2023, Lehmann et al. discussed how the use of the mixed pressure-velocity formula in PINN can achieve higher computational accuracy and better convergence for the simulation of seepage in heterogeneous porous media than the pure pressure formula, and found that the use of three independent neural networks to approximate the pressure (p) and two velocity components (v x 、v y ) of PINN (PV-3-PINN, such as Figure 3 Compared to using only one neural network to approximate the pressure (p) and the two velocity components (v x 、v y) can achieve better computational performance. Therefore, the present invention will focus on testing whether the PIKAN framework has similar conclusions as those in the PINN framework in subsequent embodiments, and compare the computational performance of PIKAN and PINN.

[0106] Still Figure 2 As an example, the PV-3-PIKAN using the mixed pressure-velocity formula is shown in Equations (15-18) and (9). The loss function when the mixed pressure-velocity formula is used is the residual of the mass conservation equation: Residuals from Darcy's Law and the residual L of the relevant boundary conditions DB , L NB In the present invention, the PIKAN model in the present invention is implemented by the ChebyKAN solution framework based on Pytorch, and the two-step strategy commonly used in PINN is used to train the neural network: first, the Adam algorithm and the specified number of iterations are used to minimize the loss function; then, the L-BFGS algorithm is used to achieve the final convergence of the loss function.

[0107]

[0108] The following is a comparison of the computational performance of the three methods through specific examples. The first method uses a mixed pressure-velocity formula and uses three neural networks to approximate pressure, v x and v y The second method is the PINN, abbreviated as PV-3-PINN, the second method is the PIKAN using the pure pressure formula, abbreviated as P-PIKAN, and the third method is Figure 2 The following intends to test whether the mixed pressure-velocity formula in the PIKAN framework has advantages over the pure pressure formula in dealing with heterogeneous domains by comparing the computational performance of PV-3-PIKAN and PIKAN, and to test whether KAN has certain advantages over MLP in the problem studied in the present invention by comparing the computational performance of PV-3-PIKAN and PV-3-PINN.

[0109] Embodiment 1

[0110] This example is a single-phase flow in a two-dimensional heterogeneous domain, such as Figure 3 As shown in (a), the rectangular area is 100 meters long and 50 meters wide. The left and right sides are at constant pressures of 10 MPa and 5 MPa respectively, and the upper and lower sides are closed. The analytical expression of heterogeneous permeability is k(x, y) = e x50 , we can know that the exact solution of the pressure distribution in this embodiment is:

[0111]

[0112] Among them, v x (x, y) is the velocity function in the x direction in this embodiment; v y (x, y) is the velocity function in the y direction in this embodiment.

[0113] Table 1 lists the relevant parameters of PV-3-PINN, P-PINN, and PV-3-PIKAN constructed in this embodiment. Figure 3 As shown in (b), the Latin hypercube sampling strategy is used to randomly generate 214 and 210 collocation points inside and on the boundary of the computational domain, respectively. The weights of each loss term in equation (7) or equation (15) are set to 1. For the training process, the loss function is first minimized by the Adam optimizer in 15,000 iterations, and then further optimized by the L-BFGS optimizer in 5,000 iterations. Figure 4 The total loss function and its sub-items (including PDE loss, Neumann boundary condition loss and Dirichlet boundary condition loss) of the three methods during training are compared with the number of optimization steps. Cartesian points are selected in the computational domain (1m spacing in the x direction, 1m spacing in the y direction, a total of 101×51=5151 points). Figure 5 , Figure 6 and Figure 7 The pressure and v on the Cartesian collocation points predicted by different methods are shown respectively. x and v y The distribution and the exact solution in Eq. (19) are quantitatively calculated in Table 2. Figure 5 , Figure 6 and Figure 7 The pressure and v predicted by the three methods x and v y The L2 error of the distribution. From these pictures and tables, we can see that in this embodiment, the values ​​of PV-3-PIKAN in each item of the loss function are the lowest among the three methods. Among them, PV-3-PIKAN and PV-3-PINN using the mixed pressure-velocity formula can make each item in the loss function drop significantly faster than P-PINN using the pure pressure formula, especially the loss of Neumann boundary conditions, indicating that the mixed formula can significantly improve the calculation performance of PIKAN compared to the pure pressure formula. In addition, the v calculated by PV-3-PIKAN x and v y The distribution is compared with the v calculated by PV-3-PINN which also uses the hybrid formula x and v yThe distribution is closer to the exact solution, especially the vx distribution, which reveals that the KAN-based PIKAN framework is superior to the MLP-based PINN framework in this embodiment. The L2 error values ​​calculated in Table 2 further quantitatively demonstrate the conclusions of the above analysis, indicating that the computational performance of PV-3-PIKAN in this embodiment is the best among the three methods. Compared with PV-3-PIKAN and P-PIKAN, the L2 error of the hybrid formula is only 1 / 10 of the L2 error of the pure pressure formula on average, indicating that the hybrid formula can be significantly better than the pure pressure formula. Compared with PV-3-PIKAN and PV-3-PINN, when both use the hybrid formula, the L2 error of the KAN-based PIKAN framework is only 0.4279 of that of the MLP-based PINN framework, indicating that KAN is better than MLP in this example.

[0114] Table 1 Parameters of the neural network and optimizer of the three methods in Example 1 and Example 2

[0115]

[0116] Table 2 L2 error of pressure and velocity distribution predicted by different methods in Example 1

[0117]

[0118]

[0119] Embodiment 2

[0120] like Figure 8 As shown, the only difference between this embodiment and the first embodiment is that the heterogeneous permeability in this embodiment no longer has an analytical continuous function form, but is partitioned. In the middle blue area, the permeability is 0.1mD, and the permeability of the red area is 1mD. The analytical solution to this problem is as follows:

[0121]

[0122] The three methods of this embodiment are the same as those of Example 1 in terms of the structure of the neural network, the generation of matching points, the parameters of the optimizer, etc. Fig. 9 The total loss function and its sub-items (including PDE loss, Neumann boundary condition loss and Dirichlet boundary condition loss) of the three methods during the training process are compared with the change curves of the optimization steps. The same Cartesian points as in Example 1 are selected in the computational domain to test the accuracy of the three methods. Fig.10 , Fig.11 and Fig.12 The pressure, vx and vy distributions obtained by different methods and the exact solution in equation (20) are shown respectively. Table 3 quantitatively calculates Fig.10 , Fig.11 and Fig.12 The L2 errors of various distributions predicted by the three methods are shown in Figure 2. From these pictures and tables, we can see that:

[0123] First, in this example, PV-3-PIKAN decreases faster than PV-3-PINN in various loss terms, and the v predicted by PV-3-PIKAN x and v y The distribution is better than that predicted by PV-3-PINN x and v y The distribution is visually closer to the reference solution. The L2 error shown in Table 3 also quantitatively demonstrates that in the calculation of velocity distribution, the error of PV-3-PINN is more than twice that of PV-3-PIKAN. These all indicate that the KAN-based PIKAN framework is superior to the MLP-based PINN framework in this embodiment.

[0124] Second, although Fig. 9 (a) in Figure 2 shows that P-PIKAN using the pure pressure formula can achieve the best convergence speed and the lowest loss function value among the three methods. However, this does not mean that P-PIKAN has the highest calculation accuracy. On the contrary, Fig.10 The difference between the pressure distribution predicted by P-PIKAN in (b) and the exact solution is very obvious, and the calculation error is very significant, reaching 1000 times the error of PV-3-PIKAN or PV-3-PINN using the mixed formula. This is because when the pure pressure formula is used, for the discontinuously distributed permeability in this embodiment, when a sample point is determined to be located in a certain area, the permeability in the area near the point will be a constant, so that k(x, y) in equation (3) at this point will be taken as a constant without the need for differentiation, making the equations at all sample points in the entire calculation domain essentially the same, so that P-PIKAN in this embodiment actually describes the seepage problem in the homogeneous calculation domain, and does not really consider the heterogeneous permeability distribution in this embodiment. Fig.10 (b) The pressure distribution predicted by P-PIKAN shows an obvious linear distribution with constant slope in the x direction of the entire computational domain, which also verifies that P-PIKAN actually solves the problem of the reservoir domain with homogeneous permeability distribution combined with the boundary conditions of this embodiment. This result illustrates the importance of using the mixed pressure-velocity formula instead of the pure pressure formula to deal with this kind of regional heterogeneous reservoir model.

[0125] Table 3 L2 error of pressure and velocity distribution predicted by different methods in Example 2

[0126]

[0127]

[0128] Embodiment 3

[0129] This example is a two-dimensional flow problem on a non-homogeneous domain, such as Fig.13 As shown, the computational domain can be divided into 4 regions, where the red region has a permeability of 9mD, the bright blue region has a permeability of 0.9mD, the light blue region has a permeability of 0.1mD, the lower left corner and the upper right corner are set at 10MPa and 5MPa respectively, and the other boundaries are closed. Table 4 lists the relevant parameters of PV-3-PINN, P-PINN, and PV-3-PIKAN constructed in this embodiment. It can be seen that for this more complex problem, this embodiment increases the number of neurons in each hidden layer, thereby appropriately increasing the total number of parameters of each neural network.

[0130] Fig.14 The total loss function and its sub-items (including PDE loss, Neumann boundary condition loss and Dirichlet boundary condition loss) of the three methods during the training process are compared with the change curves of the optimization steps. Cartesian points are arranged on the computational domain to test the accuracy of the three methods, and the corresponding reference solutions are obtained using COMSOL Multiphsyics. Fig.15 , Fig.16 and Fig.17 The pressure and v obtained by different methods are shown respectively. x and v y Table 5 quantitatively calculates the distribution and exact solution. Fig.15 , Fig.16 and Fig.17 From these pictures and tables, we can still see that in the two-dimensional flow problem of this embodiment, PV-3-PIKAN obviously decreases faster than PV-3-PINN in various loss terms, and the v predicted by PV-3-PIKAN is x and v y The distribution is better than that predicted by PV-3-PINN x and v y The distribution is visually closer to the reference solution. The L2 error shown in Table 5 also quantitatively demonstrates that the error of PV-3-PINN is higher than that of PV-3-PIKAN in the calculation of velocity distribution. These all indicate that the KAN-based PIKAN framework is better than the MLP-based PINN framework in this embodiment. In addition, Fig.15The pressure distribution predicted by P-PIKAN in (b) is significantly different from the exact solution, reaching 100 times the error of PV-3-PIKAN using the hybrid formula and 50 times the error of PV-3-PINN. This is because, as analyzed above, P-PIKAN using the pure pressure formula actually solves the problem of the reservoir domain with homogeneous permeability distribution combined with the boundary conditions of this embodiment, which once again supports the necessity of using the hybrid pressure-velocity formula instead of the pure pressure formula to deal with this kind of regional heterogeneous reservoir model.

[0131] Table 4 Parameters of the neural network and optimizer of the three methods in Example 3

[0132]

[0133] Table 5 L2 error of pressure and velocity distribution predicted by different methods in Example 3

[0134] method P(MPa) <![CDATA[v x (cm / s)]]> <![CDATA[v y (cm / s)]]> PV-3-PIKAN <![CDATA[1.7123×10 -2 ]]> <![CDATA[2.6216×10 -3 ]]> <![CDATA[2.3088×10 -3 ]]> PV-3-PINN <![CDATA[4.7411×10 -2 ]]> <![CDATA[2.8385×10 -3 ]]> <![CDATA[2.9015×10 -3 ]]> P-PIKAN 2.0459 <![CDATA[1.2763×10 -1 ]]> <![CDATA[1.2813×10 -1 ]]>

[0135] It is worth noting that the contents not elaborated in detail in the present invention are all prior art and are well known to those skilled in the art.

[0136] Therefore, the present invention adopts the above-mentioned PIKAN-based heterogeneous reservoir seepage simulation method, introduces the promising KAN into porous media flow simulation for the first time, and provides a preliminary reference for developing a general seepage simulation tool based on PIKAN.

[0137] Finally, it should be noted that the above embodiments are only used to illustrate the technical solution of the present invention rather than to limit it. Although the present invention has been described in detail with reference to the preferred embodiments, those skilled in the art should understand that they can still modify or replace the technical solution of the present invention with equivalents, and these modifications or equivalent replacements cannot cause the modified technical solution to deviate from the spirit and scope of the technical solution of the present invention.

Claims

1. A PIKAN-based heterogeneous reservoir seepage simulation method, characterized in that: The following steps are involved: Step S1: Use three independent KAN networks to approximate the pressure p and velocity component v respectively x and the velocity component v y PIKAN; Step S2: construct a loss function of a mixed pressure-velocity formula including the residual of the mass conservation equation, the residual of Darcy's law and the residual of related boundary conditions, and optimize the loss function.

2. The method for simulating heterogeneous reservoir seepage based on PIKAN according to claim 1, characterized in that: In step S1, three independent KAN networks are used to approximate the pressure p and velocity component v respectively. x and the velocity component v y The specific operations of PIKAN are as follows: For a smooth have: in, Φ q :→; x represents n independent variables x i The vector is composed of i=1,L,n; f(x) represents the n independent variables x i A function of x p represents the pth component of vector x; univariate function φ q,p Process the pth component of vector x and contribute a term to the sum of the qth external function; there are 2n+1 external functions, each of which has q is a univariate function; The Chebyshev polynomial function is used to replace the B-spline function; For a fully connected KAN with L+1 layers, the input layer is denoted as layer 0, n l represents the number of neurons in the lth layer, (l,j) represents the jth neuron in the lth layer, represents the activation value of the neuron; between layer l and layer l+1, there will be n l n l+1 parameterized univariate functions, where the univariate function connecting (l,j) and (l+1,i) is denoted by Then the activation value of neuron (l+1,i) is calculated as: in, represents the value of the jth neuron in the l+1th layer; represents the value of the i-th neuron in the l-th layer; n l represents the number of neurons in the lth layer; Therefore, according to the activation value x of the neuron in the l+1th layer l Calculate the vector x consisting of the activation values ​​of the neurons in the l+1th layer l+1 The matrix form expression of is: x l+1 =T l x l (3); Among them, T l represents the function matrix corresponding to the lth KAN layer, T l The i-th row and j-th column of is a Chebyshev polynomial function; All Chebyshev polynomial functions in KAN All learnable parameters are denoted as θ, l = 0, 1, L, L-1, j = 1, 2, L, n l ,i=1,2,L,n l+1 , we get: When the input layer data is x 0 When , the output layer neuron value of the fully connected KAN with L+1 layers is calculated as: in, Represents the composition of functions.

3. The method for simulating heterogeneous reservoir seepage based on PIKAN according to claim 2, characterized in that: In step S2, the loss function of the mixed pressure-velocity formula is: Among them, L(θ) represents the loss function; ω1, ω2, ω3, ω4 all represent weighting coefficients; represents the residual of the mass conservation equation; represents the residual of Darcy's law; L DB represents the residual of the Dirichlet boundary condition; L NB represents the residual of the Neumann boundary condition.

4. The method for simulating heterogeneous reservoir seepage based on PIKAN according to claim 3, characterized in that: Residual of the mass conservation equation The expression is as follows: Among them, N in Indicates the number of points inside the computational domain; (x i ,y i ) represents the coordinates of a certain point; x and y represent the spatial independent variables of the function.

5. The method for simulating heterogeneous reservoir seepage based on PIKAN according to claim 3, characterized in that: Residuals from Darcy's Law The expression is as follows: Here, k represents the permeability.

6. The method for simulating heterogeneous reservoir seepage based on PIKAN according to claim 3, characterized in that: The residual of the Dirichlet boundary condition is expressed as follows: Among them, N DB represents the number of collocation points on the Dirichlet frontier; p represents the number of collocation points (x i ,y i ) is the PIKAN prediction value on ; g represents the Dirichlet boundary condition value.

7. The method for simulating heterogeneous reservoir seepage based on PIKAN according to claim 3, characterized in that: The residual expression of the Neumann boundary condition is as follows: Among them, N NB represents the number of collocation points on the Neumann boundary; h x represents the x-direction velocity value given by the Neumann boundary condition; h y Represents the y-direction velocity value given by the Neumann boundary condition.

Images