A method and system for implementing polymer-flood percolation simulations
By combining Gaussian process regression and physical information neural network, the numerical dispersion problem caused by grid dependence in traditional reservoir numerical simulation is solved, and efficient and stable simulation of polymer flooding seepage process is achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- GONGHE IND SOFTWARE (SHANGHAI) CO LTD
- Filing Date
- 2026-05-11
- Publication Date
- 2026-07-31
AI Technical Summary
Traditional reservoir numerical simulation methods rely on rigid grid partitioning, which leads to severe numerical dispersion. The nonlinear iterative solution process is time-consuming and has a computational bottleneck, making it difficult to accurately predict polymer flooding flow processes.
A Gaussian process regression algorithm is used to construct a fluid apparent viscosity and relative permeability mapping model. Combined with a physical information neural network, the residuals of the partial differential equations are transformed into a loss function. The weight parameters are updated using a neural network optimization algorithm to achieve rapid meshless reconstruction of the seepage physical field.
It eliminates the grid dependency of traditional simulations, improves the input fidelity of the seepage control equations, ensures efficient and stable convergence of polymer-driven seepage simulations, and enables rapid reconstruction of complex physical fields.
Smart Images
Figure CN122491131A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of reservoir numerical simulation, and in particular to a method and system for simulating polymer flooding seepage. Background Technology
[0002] Polymer flooding is a core technology for enhancing oil recovery, and its seepage process involves a complex coupling of multiphase fluid transport and physicochemical changes in polymer macromolecules. In the field of reservoir numerical simulation, to accurately predict subsurface fluid dynamics, it is usually necessary to establish a set of nonlinear partial differential equations based on the law of mass conservation, which includes two-phase seepage of oil and water and transport of polymer components. Traditional numerical calculation methods mainly rely on discretization techniques such as finite difference to divide the continuous reservoir space into a massive discrete grid, and then establish and solve a huge set of algebraic equations to track the spatiotemporal evolution of pressure field, saturation field, and concentration.
[0003] However, existing technologies rely heavily on rigid mesh partitioning, which is prone to numerical dispersion when facing complex boundaries. When solving this strongly coupled nonlinear system, each step requires calculating and solving the Jacobian matrix with extremely high dimensions, resulting in large memory consumption and long processing time in the nonlinear iterative inversion process, creating a computational bottleneck. Summary of the Invention
[0004] To overcome the above shortcomings, this invention provides a method and system for simulating polymer-driven seepage flow, aiming to improve the numerical dispersion problem caused by the reliance on rigid grids in traditional simulations.
[0005] In a first aspect, the present invention provides the following technical solution: a method for simulating polymer-driven percolation, comprising: S1. Obtain core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity and relative permeability, and use Gaussian process regression algorithm to fit the experimental data to construct a mapping model of polymer concentration and shear rate with fluid apparent viscosity and relative permeability correction coefficients. S2. Based on the law of conservation of mass, establish multiphase flow equations and polymer component transport equations. Substitute the fluid apparent viscosity and relative permeability correction coefficients determined by the mapping model into the multiphase flow equations and the polymer component transport equations to construct a set of polymer-driven flow partial differential equations. S3. Construct a physical information neural network, using the reservoir spatial coordinates and time as input nodes of the physical information neural network, and the reservoir pressure field, water saturation field and polymer concentration field as output nodes. S4. The calculation residuals of the polymer-driven seepage flow partial differential equation system and the preset initial boundary conditions are transformed into the loss function of the physical information neural network, and the loss function is minimized by the neural network optimization algorithm to update the weight parameters of the physical information neural network. S5. When the loss function meets the preset convergence condition, the iteration stops, and the target physical field is calculated and output using the trained physical information neural network. The target physical field includes the target reservoir spatial coordinates and the reservoir pressure field, water saturation field and polymer concentration field at the target time.
[0006] Preferably, in step S1, the step of obtaining the core displacement experimental data includes: The dynamic time-series rheological data acquired by the core displacement sensor was smoothed using a moving average filtering algorithm to obtain smoothed dynamic time-series rheological data. The Z-score normalization method was used to perform dimensionless processing on the noise-reduced polymer concentration and shear rate data, and the physical parameters of different dimensions were uniformly transformed into a standard normal distribution space with a mean of zero and a variance of one. The standardized feature parameters are used as input vectors, and the corresponding fluid apparent viscosity measurements and the relative permeability correction coefficients calculated from the core pressure drop curves are used as output vectors to construct a sample dataset for model training.
[0007] Preferably, in step S1, the step of fitting the experimental data using a Gaussian process regression algorithm includes: A covariance matrix containing a radial basis kernel function is constructed, and a nonlinear function is established with polymer concentration and shear rate as independent variables and fluid apparent viscosity and relative permeability correction coefficients as dependent variables. A log-marginal likelihood objective function is established, and the gradient descent algorithm is used to iteratively optimize the length scale of the covariance matrix and the signal variance hyperparameter to obtain the target hyperparameter of the Gaussian process. Based on the Bayesian inference framework, the predicted mean of the test points is calculated using the optimized covariance matrix, and a mapping model is generated between the polymer concentration and shear rate and the correction coefficients of the fluid apparent viscosity and relative permeability.
[0008] Preferably, in step S2, the step of constructing the partial differential equations for polymer-driven percolation includes: Extract the predicted mean from the output of the mapping model to construct the fluid apparent viscosity field matrix and residual drag coefficient matrix that vary with spatial location; Substituting the fluid apparent viscosity field matrix into the water phase mobility calculation formula, the local seepage velocity of the water phase is calculated, and substituting the residual resistance coefficient matrix into the absolute permeability calculation formula, the corrected absolute permeability is calculated. By introducing the Langmuir isotherm model formula and the inaccessible pore volume parameter, and combining it with the seepage velocity field calculated above, a set of partial differential equations for polymer-driven seepage coupled with multiphase fluids is assembled to generate the system.
[0009] Preferably, in step S3, the step of constructing the physical information neural network includes: A feedforward neural network containing multiple fully connected layers is constructed as the basic network architecture, and the reservoir spatial coordinates and time are concatenated into an input feature vector and then fed into the basic network architecture. The hyperbolic tangent function is configured as the neuron activation function in the hidden layer of the basic network architecture; In the output layer of the basic network architecture, a Sigmoid activation function is configured for the output node of the water saturation field, and the output value is the predicted water saturation value in the range of 0 to 1.
[0010] Preferably, in step S4, the step of converting the loss function into the physical information neural network includes: Using the automatic differentiation technique in the deep learning framework, the first-order spatial partial derivative, the first-order temporal partial derivative, and the second-order spatial partial derivative of the output variable of the physical information neural network with respect to the input spatiotemporal coordinates are calculated analytically. Substitute the obtained partial derivatives into the polymer-driven seepage partial differential equation system to calculate the equation imbalance at the spatiotemporal configuration point and construct the residual loss term of the partial differential equation. The mean square error of the network's output at the initial time with respect to the initial geological conditions, and the mean square error of the output at the spatial boundary with respect to the preset boundary conditions are calculated respectively. The initial boundary value loss term is constructed and added to the residual loss term of the partial differential equation to obtain the joint loss function.
[0011] Preferably, in step S4, the step of minimizing the loss function using a neural network optimization algorithm includes: The loss function is divided into a partial differential equation residual loss term that characterizes the calculated residual and an initial boundary value loss term that characterizes the initial boundary conditions. Initial constant weight coefficients are set for the partial differential equation residual loss term and the initial boundary value loss term, respectively. In the initial stage of network training, the Adam optimization algorithm is used to globally iteratively update the weight parameters of the physical information neural network based on a set learning rate. When the rate of change of the loss function is lower than a set threshold or the preset number of iterations is reached, the algorithm switches to the L-BFGS optimization algorithm to perform local iterative optimization using the Hessian matrix approximation information.
[0012] Preferably, in step S5, the step of calculating and outputting the target physical field using the trained physical information neural network includes: Based on the set simulated geological region and time step, uniform sampling is performed in the spatial and temporal domains to generate a discrete coordinate lattice of the target spatiotemporal location. The discrete coordinate point array is used as an input tensor and fed into the trained physical information neural network for forward propagation calculation to obtain the network prediction value at the corresponding coordinate point. The predicted values output by the network are denormalized to calculate the physical field data with the original physical dimensions, and the output of the physical field data is saved as a three-dimensional physical flow field data volume.
[0013] Secondly, the present invention provides the following technical solution: a polymer-driven percolation simulation system, comprising: The mapping model construction module is used to acquire core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity and relative permeability, and to fit the experimental data using a Gaussian process regression algorithm to construct a mapping model of polymer concentration and shear rate with correction coefficients of fluid apparent viscosity and relative permeability. The equation system construction module is used to establish multiphase flow equations and polymer component transport equations based on the law of conservation of mass. The fluid apparent viscosity and relative permeability correction coefficients determined by the mapping model are substituted into the multiphase flow equations and the polymer component transport equations to construct a set of polymer-driven flow partial differential equations. The network construction module is used to construct a physical information neural network, taking the reservoir spatial coordinates and time as input nodes of the physical information neural network, and the reservoir pressure field, water saturation field and polymer concentration field as output nodes. The network training module is used to transform the computational residuals of the polymer-driven percolation partial differential equations and the preset initial boundary conditions into the loss function of the physical information neural network, and to minimize the loss function using a neural network optimization algorithm to update the weight parameters of the physical information neural network. The simulation output module is used to stop iterating when the loss function meets the preset convergence condition, and to calculate and output the target physical field using the trained physical information neural network. The target physical field includes the target reservoir spatial coordinates and the reservoir pressure field, water saturation field and polymer concentration field at the target time.
[0014] The present invention has the following beneficial effects: 1. In this invention, by constructing a physical information neural network and transforming the residuals of partial differential equations into loss functions, the computational bottleneck of traditional simulation relying on discrete grids and solving large Jacobian matrices is overcome. The time-consuming nonlinear iterative inversion is transformed into efficient forward tensor derivation, realizing the meshless rapid reconstruction of complex physical fields driven by polymers.
[0015] 2. In this invention, the Gaussian process regression algorithm is used to perform nonparametric fitting on the core displacement experimental data, which improves the problem that traditional empirical formulas are difficult to accurately characterize the abrupt changes in fluid rheology, characterizes the nonlinear effects of concentration and shear rate on apparent viscosity and drag coefficient, and improves the input fidelity of the bottom-layer physical property parameters of the seepage control equation.
[0016] 3. In this invention, an activation function with high-order continuous differentiability and numerical truncation characteristics is configured for the multiphase coupling characteristics of polymer-driven seepage simulation. A two-stage optimization algorithm combining global exploration and local Hessian matrix approximation optimization is introduced. This not only ensures the continuity of the flow field when automatically differentiating the spatial gradient, but also strictly constrains the physical boundary of water saturation, ensuring the efficient and stable convergence of the entire polymer-driven seepage simulation solution system. Attached Figure Description
[0017] Figure 1 This is a flowchart of a polymer-driven percolation simulation method proposed in this invention; Figure 2 This is an architecture diagram of a polymer-driven percolation simulation system proposed in this invention. Detailed Implementation
[0018] The technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.
[0019] Example 1 In a first embodiment of the present invention, the present invention provides a method for simulating polymer-driven percolation, such as... Figure 1 As shown, it includes: S1. Obtain core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity and relative permeability. Use Gaussian process regression algorithm to fit the experimental data and construct a mapping model between polymer concentration and shear rate and fluid apparent viscosity and relative permeability correction coefficients. Furthermore, step S1, the step of obtaining core displacement experimental data, includes: The dynamic time-series rheological data acquired by the core displacement sensor was smoothed using a moving average filtering algorithm to obtain smoothed dynamic time-series rheological data. The Z-score normalization method was used to perform dimensionless processing on the noise-reduced polymer concentration and shear rate data, and the physical parameters of different dimensions were uniformly transformed into a standard normal distribution space with a mean of zero and a variance of one. The standardized feature parameters are used as input vectors, and the corresponding fluid apparent viscosity measurements and the relative permeability correction coefficients calculated from the core pressure drop curves are used as output vectors to construct a sample dataset for model training.
[0020] Furthermore, step S1, which involves fitting the experimental data using a Gaussian process regression algorithm, includes: A covariance matrix containing a radial basis kernel function is constructed, and a nonlinear function is established with polymer concentration and shear rate as independent variables and fluid apparent viscosity and relative permeability correction coefficients as dependent variables. A log-marginal likelihood objective function is established, and the gradient descent algorithm is used to iteratively optimize the length scale of the covariance matrix and the signal variance hyperparameter to obtain the target hyperparameter of the Gaussian process. Based on the Bayesian inference framework, the predicted mean of the test points is calculated using the optimized covariance matrix, and a mapping model is generated between polymer concentration and shear rate and the correction coefficients of fluid apparent viscosity and relative permeability.
[0021] Specifically, core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity, and relative permeability were acquired. The dynamic time-series rheological data acquired by the core displacement sensor was smoothed using a moving average filtering algorithm to obtain smoothed dynamic time-series rheological data. Let the original dynamic time-series rheological data sequence be... The sliding window size is Smoothed data The calculation formula is: ; in Representing the The fluid apparent viscosity or pressure drop measurement at each sampling time. This represents the smoothed measurement value. The algorithm is used to perform point-by-point sliding calculations, eliminating outliers and high-frequency noise during the measurement process.
[0022] Z-score normalization was used to perform dimensionless processing on the denoised polymer concentration and shear rate data. Let the denoised polymer concentration or shear rate sequence be... The mean of this sequence is The standard deviation is Standardized feature parameters The calculation formula is: ; This calculation transforms physical parameters of different dimensions into a standard normal distribution space with a mean of zero and a variance of one. The standardized feature parameters are used as the input vector, and the corresponding fluid apparent viscosity measurements and the relative permeability correction coefficient calculated from the core pressure drop curve are used as the output vector. These are combined to construct a sample dataset for model training. Let the input vector be... The output vector is Input vector The output vector includes two-dimensional features of normalized polymer concentration and shear rate. It includes two-dimensional features such as the fluid apparent viscosity and relative permeability correction factor.
[0023] Gaussian process regression algorithm was used to fit the experimental data. A covariance matrix including a radial basis function kernel was constructed, and a nonlinear function was established with polymer concentration and shear rate as independent variables and fluid apparent viscosity and relative permeability correction coefficients as dependent variables. Radial basis function kernel. The calculation formula is: ; in and They represent the first and second elements in the sample dataset, respectively. The and the first A polymer concentration and shear rate input vector This represents the squared Euclidean distance between two input vectors. This represents the signal variance hyperparameter. This represents the length scale hyperparameter.
[0024] A log-marginal likelihood objective function is established, and the gradient descent algorithm is used to iteratively optimize the length scale of the covariance matrix and the signal variance hyperparameters. Log-marginal likelihood objective function The calculation formula is: ; in This represents the set of output vectors from the training sample dataset. This represents the covariance matrix calculated by the radial basis function based on the input vector. Represents the variance of observation noise. Represents the identity matrix. This represents the total number of training samples. The objective hyperparameters of the Gaussian process are obtained by maximizing this objective function to update the signal variance hyperparameter and the length scaling hyperparameter.
[0025] Based on a Bayesian inference framework, the predicted mean of test points is calculated using the optimized covariance matrix, generating a mapping model between polymer concentration and shear rate and correction coefficients for fluid apparent viscosity and relative permeability. Let the new test point input vector be... The corresponding predicted mean The calculation formula is: ; in The covariance vector, calculated using the optimized radial basis function, represents the relationship between the polymer concentration and shear rate input vectors at the new test point and the input vectors of all training samples. This predicted mean is the predicted value of the fluid apparent viscosity or relative permeability correction coefficient output by the mapping model.
[0026] The above steps improve the accuracy of fitting fluid apparent viscosity and relative permeability parameters by filtering and denoising the core displacement rheological data and extracting dimensionless features, and by constructing a probabilistic prediction model based on a Gaussian process regression algorithm with radial basis kernel function.
[0027] S2. Based on the law of conservation of mass, establish multiphase flow equations and polymer component transport equations. Substitute the fluid apparent viscosity and relative permeability correction coefficients determined by the mapping model into the multiphase flow equations and polymer component transport equations to construct a set of polymer-driven flow partial differential equations. Further, step S2, the step of constructing the partial differential equations for polymer-driven percolation, includes: Extract the predicted mean from the output of the mapping model to construct the fluid apparent viscosity field matrix and residual drag coefficient matrix that vary with spatial location; Substituting the fluid apparent viscosity field matrix into the formula for calculating the water phase mobility, the local seepage velocity of the water phase is calculated. Substituting the residual resistance coefficient matrix into the formula for calculating the absolute permeability, the corrected absolute permeability is calculated. By introducing the Langmuir isotherm model formula and the inaccessible pore volume parameter, and combining it with the seepage velocity field calculated above, a set of partial differential equations for polymer-driven seepage coupled with multiphase fluids is assembled.
[0028] Specifically, the predicted mean value output by the mapping model is extracted to construct the fluid apparent viscosity field matrix and residual drag coefficient matrix that vary with spatial location. Based on the reservoir spatial grid partitioning results, the predicted fluid apparent viscosity scalar value and residual drag coefficient scalar value at each discrete grid node are arranged according to the spatial topology to generate a three-dimensional fluid apparent viscosity field matrix and residual drag coefficient matrix consistent with the dimension of the reservoir computational domain, which are used as inputs for heterogeneous physical property parameters at the current simulation time step.
[0029] The local seepage velocity of the water phase is calculated by substituting the fluid apparent viscosity field matrix into the formula for calculating the water phase mobility, and the corrected absolute permeability is calculated by substituting the residual drag coefficient matrix into the formula for calculating the absolute permeability. According to Darcy's law, the local seepage velocity of the water phase... The calculation formula is: ; The corrected absolute penetration rate The calculation formula is: ; in the formula Represents the initial absolute permeability tensor of the reservoir. Represents the residual drag coefficient matrix. Represents the relative permeability of the water phase. Represents the apparent viscosity field matrix of the fluid. This represents the potential gradient of the aqueous phase fluid. The above calculations quantitatively transform the dynamic plugging and thickening effect of the polymer on the pores into flow field velocity distribution characteristics.
[0030] By introducing the Langmuir isotherm model formula and inaccessible pore volume parameters, and combining them with the seepage velocity field calculated above, a set of partial differential equations for polymer-driven seepage coupled with multiphase fluids is generated. The Langmuir isotherm model calculation formula is as follows: ; in the formula This represents the concentration of polymers adsorbed on the rock surface. This represents the polymer concentration in the aqueous phase. and This represents the adsorption constant calibrated by core experiments. A polymer component transport equation is constructed using the law of conservation of mass, incorporating convection, dispersion, adsorption, and inaccessible pore correction terms: ; in the formula Represents the dispersion tensor of polymers in porous media. Represents the source and sink of aqueous phase injection or production. This represents the concentration of the injected polymer. Represents simulation time. Represents the initial porosity of the reservoir. Represents the inaccessible pore volume parameter. Represents water saturation. This represents the density of the reservoir rock. This equation is then nonlinearly combined with the traditional oil-water two-phase flow control equations to generate a set of partial differential equations for polymer-driven flow that incorporates multiphase flow and component transport.
[0031] The above steps deeply couple the dynamic rheological parameters obtained through data-driven acquisition with the physical laws of multiphase and multicomponent seepage, and construct a closed set of partial differential equations that conform to the laws of microscopic displacement mechanics.
[0032] S3. Construct a physical information neural network, using reservoir spatial coordinates and time as input nodes, and reservoir pressure field, water saturation field, and polymer concentration field as output nodes. Furthermore, step S3, the step of constructing the physical information neural network, includes: A feedforward neural network containing multiple fully connected layers is constructed as the basic network architecture, and the reservoir spatial coordinates and time are concatenated into an input feature vector and then fed into the basic network architecture. Configure the hyperbolic tangent function as the neuron activation function in the hidden layers of the basic network architecture; In the output layer of the basic network architecture, a Sigmoid activation function is configured for the output nodes of the water saturation field, and the output value is the predicted water saturation value in the range of 0 to 1.
[0033] Specifically, a feedforward neural network containing multiple fully connected layers is constructed as the basic network architecture of the physical information neural network. The reservoir spatial coordinates of the target computational block and the current simulation time are extracted as input features to the network. Let the three-dimensional components of the reservoir spatial coordinates be the horizontal axis... ordinate and vertical coordinates The simulation time is These four physical scalars are concatenated in order to form the input feature vector. The reservoir pressure field, water saturation field, and polymer concentration field are used as output nodes in the basic network architecture to form the output feature vector. The output feature vector fully contains the predicted reservoir pressure value corresponding to the input spatiotemporal location. and predicted water saturation value and polymer concentration prediction values .
[0034] A deep neural network (DNN) with multiple fully connected layers is constructed as the basic network architecture for a physical information neural network. The specific architecture of this deep neural network is configured to include... The number of fully connected hidden layers is preferably [missing information]. Layers, each hidden layer contains The first neuron node is used. The concatenated input feature vector is fed into this deep neural network for forward propagation. The hyperbolic tangent function is configured as the neuron activation function in each hidden layer of the deep neural network. Let the first neuron node be... The input feature vectors of the hidden layers are The weight matrix of the current layer is The bias vector is Then the first The output vector of each hidden layer The calculation formula is: ; The hyperbolic tangent function For any input scalar The formula for nonlinear mapping is: ; In the initialization phase of the forward propagation algorithm, the input state of the first hidden layer is set as follows: The application of the hyperbolic tangent function in this computational process ensures that the gradient propagation flow of the computational graph remains continuous when the neural network subsequently solves for the higher-order spatial partial derivatives of the reservoir pressure field and polymer concentration field using an automatic differentiation engine.
[0035] In the output layer of the basic network architecture, a Sigmoid activation function is configured for numerical truncation of the water saturation field output nodes. Let the linear weighted calculation result corresponding to the water saturation physical field nodes in the last layer of the feedforward neural network be... The final extracted predicted water saturation value The calculation formula is: ; Taking advantage of the non-negativity of the natural exponential function, this formula algorithm forces the output value to strictly fall within a certain range. to The predicted water saturation value within the numerical range. For reservoir pressure field and polymer concentration field output nodes that do not have strict upper and lower bound physical constraints, the linear calculation result of the last layer is directly output.
[0036] The above steps construct a deep neural network forward derivation channel that maps the macroscopic spatiotemporal coordinates of the reservoir to the microscopic multiphase fluid physical field, providing numerical approximation basis functions for partial differential equations that have high-order smooth differentiation properties and conform to the physical constraints of pore fluid volume ratio.
[0037] S4. The calculation residuals of the polymer-driven seepage flow partial differential equation system and the preset initial boundary conditions are transformed into the loss function of the physical information neural network, and the loss function is minimized by the neural network optimization algorithm to update the weight parameters of the physical information neural network. Further, step S4, the step of converting the loss function into a physical information neural network, includes: Using the automatic differentiation technique in the deep learning framework, we analyze the first-order spatial partial derivative, first-order temporal partial derivative, and second-order spatial partial derivative of the output variable of the computational physical information neural network with respect to the input spatiotemporal coordinates. Substitute the obtained partial derivatives into the polymer-driven seepage partial differential equation system, calculate the equation imbalance at the spatiotemporal configuration point, and construct the residual loss term of the partial differential equation. The mean square error of the network output at the initial time and the initial geological conditions, and the mean square error of the output at the spatial boundary and the preset boundary conditions are calculated respectively. The initial boundary value loss term is constructed and added to the partial differential equation residual loss term to obtain the joint loss function.
[0038] Furthermore, step S4, which involves minimizing the loss function using a neural network optimization algorithm, includes: The loss function is divided into a partial differential equation residual loss term that characterizes the computation of residuals and an initial boundary value loss term that characterizes the initial boundary conditions. Initial constant weight coefficients are set for the partial differential equation residual loss term and the initial boundary value loss term, respectively. In the initial stage of network training, the Adam optimization algorithm is used to globally iteratively update the weight parameters of the physical information neural network based on a set learning rate. When the rate of change of the loss function is lower than the set threshold or the preset number of iterations is reached, the algorithm switches to the L-BFGS optimization algorithm to perform local iterative optimization using the Hessian matrix approximation information.
[0039] Specifically, the automatic differentiation technique within a deep learning framework is used to analytically calculate the partial derivatives of the output variables of the physical information neural network with respect to the input spatiotemporal coordinates. A chain rule for differentiation is constructed based on the computational graph to accurately extract the reservoir pressure prediction values output by the neural network. and predicted water saturation value and polymer concentration prediction values For input spatial coordinates , , The first and second spatial partial derivatives of these output variables are calculated simultaneously with respect to the input simulation time. The first-order time partial derivatives are used. These partial derivatives are then substituted into the polymer-driven seepage partial differential equations to calculate the unbalance at the spatiotemporal configuration point, constructing the residual loss term of the partial differential equations. Let the residual loss term of the partial differential equations be... The calculation formula is as follows: ; in This represents the total number of spatiotemporal configuration points randomly sampled within the reservoir computational domain. , , These represent the residual values calculated by substituting the network predictions and their partial derivatives into the oil phase flow equation, the water phase flow equation, and the polymer component transport equation, respectively. and Representing the Spatial coordinates and time of each configuration point.
[0040] The initial-boundary value loss term is constructed by calculating the mean square error of the network's output with respect to the initial geological conditions at the initial time, and the mean square error of the output with respect to the preset boundary conditions at the spatial boundary. Let the initial-boundary value loss term be... The calculation formula is as follows: ; in This represents the total number of initial configuration points sampled at zero point in the simulation time. This represents the network's predicted physics vector at the initial point. A vector representing the true values of the initial pressure saturation and concentration of the reservoir; This represents the total number of boundary configuration points sampled at the boundary of the computational domain. This represents the network's predicted vector at the boundary points. This represents the setpoint of the physical quantity corresponding to constant pressure and constant production or closed boundary conditions. The joint loss function is obtained by adding the residual loss term and the initial boundary value loss term of the partial differential equation.
[0041] The loss function is divided into a partial differential equation residual loss term characterizing the computation of residuals and an initial boundary value loss term characterizing the initial boundary conditions, and initial constant weight coefficients are set for each term. Let the joint loss function be... The calculation formula is as follows: ; in and and represent the weight coefficients of the residual loss term and the initial boundary value loss term in the partial differential equation, respectively. In the initial stage of network training, the Adam optimization algorithm is used to globally iteratively update the weight parameters of the physical information neural network based on a set learning rate. Let the set of all weight parameters of the network be . The learning rate is , No. The gradient vector of the next iteration is The decay rates of the first-order momentum and the second-order momentum are respectively and The parameter update formula is: ; in and This represents the first and second moment estimates after bias correction. To prevent smooth terms with a denominator of zero.
[0042] When the rate of change of the loss function falls below a set threshold or reaches a preset number of iterations, the algorithm switches to the L-BFGS optimization algorithm, which uses the Hessian matrix approximation information for local iterative optimization. The L-BFGS algorithm determines the search direction by approximating the second derivative information of the loss function. And perform a line search to find the optimal step size. Update the parameters: ; The formula for calculating the search direction is approximately: , It is the inverse approximation matrix of the quasi-Newton-Hessen matrix based on a small number of previous iterations of gradients.
[0043] The above steps utilize automatic differentiation to precisely embed the physical laws of seepage into the neural network loss function, and achieve rapid and stable convergence of network parameters through a two-stage hybrid optimization strategy.
[0044] S5. When the loss function meets the preset convergence condition, stop the iteration, use the trained physical information neural network to calculate and output the target physical field, which includes the target reservoir spatial coordinates and the reservoir pressure field, water saturation field and polymer concentration field at the target time.
[0045] Further, step S5, which involves calculating and outputting the target physical field using the trained physical information neural network, includes: Based on the set simulated geological region and time step, uniform sampling is performed in the spatial and temporal domains to generate a discrete coordinate lattice of the target spatiotemporal location. The discrete coordinate point array is used as an input tensor and fed into the trained physical information neural network for forward propagation calculation to obtain the network prediction value at the corresponding coordinate point. The predicted values output by the network are denormalized to calculate the physical field data with the original physical dimensions, and the output of the physical field data is saved as a three-dimensional physical flow field data volume.
[0046] Specifically, the numerical evolution of the joint loss function is monitored. Iteration stops when the relative rate of change of the joint loss function falls below a preset convergence tolerance within a set number of iterations. At this point, the physical information neural network completes training and fixes the weight parameters of all layers. Based on the set reservoir simulation geological region and time step, discrete coordinate lattices of the target spatiotemporal location are generated by uniformly sampling in both the spatial and temporal domains. The spatial boundary of the target geological region is defined as a three-dimensional Cartesian coordinate interval, and the total simulation time is defined as a one-dimensional time interval. This multi-dimensional interval is divided according to the set spatial grid step size and time step size, generating a test set input tensor composed of discrete spatiotemporal coordinate nodes. The discrete points in the test set input tensor are represented as a four-dimensional spatiotemporal vector containing three-dimensional spatial horizontal, vertical, and horizontal coordinates and one-dimensional time. .
[0047] A discrete coordinate matrix is input as an input tensor to a trained physical information neural network for forward propagation calculations to obtain the network prediction values at the corresponding coordinate points. Let the forward mapping function of the converged neural network be... ,in This represents the set of optimal weight matrices and bias vectors for the fixed fully connected layers of the network. The test set is input as a tensor. Input this forward mapping function to obtain the network prediction tensor in a dimensionless state. The calculation formula is as follows: ; The dimensionless prediction tensor obtained by this calculation It fully includes a dimensionless reservoir pressure matrix, a dimensionless water saturation matrix, and a dimensionless polymer concentration matrix that are strictly corresponding to the input discrete coordinate lattice.
[0048] Denormalizing the dimensionless predicted tensor output by the network yields physical field data with the original physical dimensions. Let a dimensionless physical quantity of one dimension in the network output be... The mean of the original samples obtained during the model training data preprocessing stage for this type of physical quantity is... The original sample standard deviation is Then the predicted value of the physical quantity with the original physical dimensions is restored. The calculation formula is: ; Through this inverse linear mapping calculation, the abstract neural network output node feature tensor is accurately restored to a reservoir pressure field with megapascals, a water saturation field with percentage values, and a polymer concentration field with milligrams per liter in real oil reservoir engineering. The aforementioned seepage physical field data with real physical dimensions are then subjected to three-dimensional spatial topological reconstruction according to a time evolution sequence, and output and saved as a three-dimensional physical flow field data volume.
[0049] The above steps replace the time-consuming implicit inversion operation of large sparse matrices in traditional numerical simulation with pure forward tensor multiplication derivation in neural networks, realizing rapid analysis of the physical field of complex seepage in polymer-driven flow and high-fidelity reconstruction of flow field data.
[0050] Example 2 Traditional techniques heavily rely on rigid mesh generation, which easily leads to numerical dispersion when dealing with complex boundaries. When solving this strongly coupled nonlinear system, each step requires calculating and solving a high-dimensional Jacobian matrix, resulting in high memory consumption and long processing time during the nonlinear iterative inversion process, creating a computational bottleneck. To address these issues, this invention provides a polymer-driven seepage simulation system, the structure of which is as follows: Figure 2 As shown. The specific implementation process of this system is as follows: The mapping model construction module is used to acquire core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity and relative permeability, and to fit the experimental data using a Gaussian process regression algorithm to construct a mapping model of polymer concentration and shear rate with correction coefficients of fluid apparent viscosity and relative permeability. The equation system construction module is used to establish multiphase flow equations and polymer component transport equations based on the law of conservation of mass. The fluid apparent viscosity and relative permeability correction coefficients determined by the mapping model are substituted into the multiphase flow equations and the polymer component transport equations to construct a set of polymer-driven flow partial differential equations. The network construction module is used to construct a physical information neural network, taking the reservoir spatial coordinates and time as input nodes of the physical information neural network, and the reservoir pressure field, water saturation field and polymer concentration field as output nodes. The network training module is used to transform the computational residuals of the polymer-driven percolation partial differential equations and the preset initial boundary conditions into the loss function of the physical information neural network, and to minimize the loss function using a neural network optimization algorithm to update the weight parameters of the physical information neural network. The simulation output module is used to stop iterating when the loss function meets the preset convergence condition, and to calculate and output the target physical field using the trained physical information neural network. The target physical field includes the target reservoir spatial coordinates and the reservoir pressure field, water saturation field and polymer concentration field at the target time.
[0051] Specifically, the mapping model construction module first filters, smooths, and dimensionlessly processes the acquired core displacement dynamic time-series rheological data. Then, using a covariance matrix with a radial basis kernel function and a log-marginal likelihood maximization objective function, it obtains the Gaussian process target hyperparameters to complete the construction of the nonlinear mapping model. The equation system construction module extracts the mean of the physical property predictions output by the model, calculates the corrected local seepage velocity and absolute permeability of the aqueous phase, and introduces the Langmuir isotherm model formula and inaccessible pore volume parameters to assemble and generate a set of polymer-driven seepage partial differential equations with strong coupling of multiphase fluids. Subsequently, the network construction module builds a multi-layer feedforward neural network, using the spliced spatiotemporal coordinates as input features, and configuring hyperbolic tangent and sigmoid functions in the hidden layers and some output layers respectively to ensure the continuity of gradient propagation and the rationality of the physical values of water saturation. Further... The network training module uses automatic differentiation technology to analytically calculate the partial derivatives of the network output with respect to the input spatiotemporal coordinates, substitutes them into a system of partial differential equations to calculate the imbalance quantity to construct the residual loss term, and combines them with the initial boundary value loss term to form a joint loss function. A two-stage hybrid optimization algorithm of Adam and L-BFGS is used to iteratively optimize the network weight parameters globally and locally. Finally, when the iteration meets the convergence condition, the simulation output module uniformly samples and generates a discrete coordinate matrix within the target geological region and time span, inputs it as a tensor into the shaped physical information neural network for pure forward propagation, and performs denormalization processing on the output network prediction values. This allows for the rapid calculation and storage of a three-dimensional physical flow field data volume with original physical dimensions, efficiently completing a high-fidelity simulation of complex polymer-driven seepage processes without the need for traditional mesh generation and large Jacobian matrix inversion.
[0052] Finally, it should be noted that the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.
Claims
1. A method for simulating polymer-driven percolation, characterized in that, include: S1. Obtain core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity and relative permeability, and use Gaussian process regression algorithm to fit the experimental data to construct a mapping model of polymer concentration and shear rate with fluid apparent viscosity and relative permeability correction coefficients. S2. Based on the law of conservation of mass, establish multiphase flow equations and polymer component transport equations. Substitute the fluid apparent viscosity and relative permeability correction coefficients determined by the mapping model into the multiphase flow equations and the polymer component transport equations to construct a set of polymer-driven flow partial differential equations. S3. Construct a physical information neural network, using the reservoir spatial coordinates and time as input nodes of the physical information neural network, and the reservoir pressure field, water saturation field and polymer concentration field as output nodes. S4. The calculation residuals of the polymer-driven seepage flow partial differential equation system and the preset initial boundary conditions are transformed into the loss function of the physical information neural network, and the loss function is minimized by the neural network optimization algorithm to update the weight parameters of the physical information neural network. S5. When the loss function meets the preset convergence condition, the iteration stops, and the target physical field is calculated and output using the trained physical information neural network. The target physical field includes the target reservoir spatial coordinates and the reservoir pressure field, water saturation field and polymer concentration field at the target time.
2. The method for simulating polymer-driven percolation according to claim 1, characterized in that, Step S1, the step of acquiring the core displacement experimental data, includes: The dynamic time-series rheological data acquired by the core displacement sensor was smoothed using a moving average filtering algorithm to obtain smoothed dynamic time-series rheological data. The Z-score normalization method was used to perform dimensionless processing on the noise-reduced polymer concentration and shear rate data, and the physical parameters of different dimensions were uniformly transformed into a standard normal distribution space with a mean of zero and a variance of one. The standardized feature parameters are used as input vectors, and the corresponding fluid apparent viscosity measurements and the relative permeability correction coefficients calculated from the core pressure drop curves are used as output vectors to construct a sample dataset for model training.
3. The method for simulating polymer-driven percolation according to claim 1, characterized in that, Step S1, the step of fitting the experimental data using the Gaussian process regression algorithm, includes: A covariance matrix containing a radial basis kernel function is constructed, and a nonlinear function is established with polymer concentration and shear rate as independent variables and fluid apparent viscosity and relative permeability correction coefficients as dependent variables. A log-marginal likelihood objective function is established, and the gradient descent algorithm is used to iteratively optimize the length scale of the covariance matrix and the signal variance hyperparameter to obtain the target hyperparameter of the Gaussian process. Based on the Bayesian inference framework, the predicted mean of the test points is calculated using the optimized covariance matrix, and a mapping model is generated between the polymer concentration and shear rate and the correction coefficients of the fluid apparent viscosity and relative permeability.
4. The method for simulating polymer-driven percolation according to claim 1, characterized in that, In step S2, the step of constructing the partial differential equations for polymer-driven percolation includes: Extract the predicted mean from the output of the mapping model to construct the fluid apparent viscosity field matrix and residual drag coefficient matrix that vary with spatial location; Substituting the fluid apparent viscosity field matrix into the water phase mobility calculation formula, the local seepage velocity of the water phase is calculated, and substituting the residual resistance coefficient matrix into the absolute permeability calculation formula, the corrected absolute permeability is calculated. By introducing the Langmuir isotherm model formula and the inaccessible pore volume parameter, and combining it with the seepage velocity field calculated above, a set of partial differential equations for polymer-driven seepage coupled with multiphase fluids is assembled to generate the equations.
5. The method for simulating polymer-driven percolation according to claim 1, characterized in that, In step S3, the step of constructing the physical information neural network includes: A feedforward neural network containing multiple fully connected layers is constructed as the basic network architecture, and the reservoir spatial coordinates and time are concatenated into an input feature vector and then fed into the basic network architecture. The hyperbolic tangent function is configured as the neuron activation function in the hidden layer of the basic network architecture; In the output layer of the basic network architecture, a Sigmoid activation function is configured for the output node of the water saturation field, and the output value is the predicted water saturation value in the range of 0 to 1.
6. The method for simulating polymer-driven percolation according to claim 1, characterized in that, Step S4, the step of converting the loss function into the physical information neural network, includes: Using the automatic differentiation technique in the deep learning framework, the first-order spatial partial derivative, the first-order temporal partial derivative, and the second-order spatial partial derivative of the output variable of the physical information neural network with respect to the input spatiotemporal coordinates are calculated analytically. Substitute the obtained partial derivatives into the polymer-driven seepage partial differential equation system to calculate the equation imbalance at the spatiotemporal configuration point and construct the residual loss term of the partial differential equation. The mean square error of the network's output at the initial time with respect to the initial geological conditions, and the mean square error of the output at the spatial boundary with respect to the preset boundary conditions are calculated respectively. The initial boundary value loss term is constructed and added to the residual loss term of the partial differential equation to obtain the joint loss function.
7. The method for simulating polymer-driven percolation according to claim 1, characterized in that, Step S4, the step of minimizing the loss function using a neural network optimization algorithm, includes: The loss function is divided into a partial differential equation residual loss term that characterizes the calculated residual and an initial boundary value loss term that characterizes the initial boundary conditions. Initial constant weight coefficients are set for the partial differential equation residual loss term and the initial boundary value loss term, respectively. In the initial stage of network training, the Adam optimization algorithm is used to globally iteratively update the weight parameters of the physical information neural network based on a set learning rate. When the rate of change of the loss function is lower than a set threshold or the preset number of iterations is reached, the algorithm switches to the L-BFGS optimization algorithm to perform local iterative optimization using the Hessian matrix approximation information.
8. The method for simulating polymer-driven percolation according to claim 1, characterized in that, Step S5, the step of calculating and outputting the target physical field using the trained physical information neural network, includes: Based on the set simulated geological region and time step, uniform sampling is performed in the spatial and temporal domains to generate a discrete coordinate lattice of the target spatiotemporal location. The discrete coordinate point array is used as an input tensor and fed into the trained physical information neural network for forward propagation calculation to obtain the network prediction value at the corresponding coordinate point. The predicted values output by the network are denormalized to calculate the physical field data with the original physical dimensions, and the output of the physical field data is saved as a three-dimensional physical flow field data volume.
9. A system for simulating polymer-driven percolation, characterized in that, A method for simulating polymer-driven percolation flow according to any one of claims 1-8, the system comprising: The mapping model construction module is used to acquire core displacement experimental data including polymer concentration, shear rate, fluid apparent viscosity and relative permeability, and to fit the experimental data using a Gaussian process regression algorithm to construct a mapping model of polymer concentration and shear rate with correction coefficients of fluid apparent viscosity and relative permeability. The equation system construction module is used to establish multiphase flow equations and polymer component transport equations based on the law of conservation of mass. The fluid apparent viscosity and relative permeability correction coefficients determined by the mapping model are substituted into the multiphase flow equations and the polymer component transport equations to construct a set of polymer-driven flow partial differential equations. The network construction module is used to construct a physical information neural network, taking the reservoir spatial coordinates and time as input nodes of the physical information neural network, and the reservoir pressure field, water saturation field and polymer concentration field as output nodes. The network training module is used to transform the computational residuals of the polymer-driven percolation partial differential equations and the preset initial boundary conditions into the loss function of the physical information neural network, and to minimize the loss function using a neural network optimization algorithm to update the weight parameters of the physical information neural network. The simulation output module is used to stop iterating when the loss function meets the preset convergence condition, and to calculate and output the target physical field using the trained physical information neural network. The target physical field includes the target reservoir spatial coordinates and the reservoir pressure field, water saturation field and polymer concentration field at the target time.