A porous medium fluid-solid coupling iterative acceleration calculation method based on fourier neural operator

Abstract

The application provides a porous medium fluid-solid coupling iterative acceleration calculation method based on a Fourier neural operator, and relates to the technical field of porous medium fluid-solid coupling simulation. The method comprises the following steps: constructing a sample data set containing heterogeneous parameters and boundary conditions, and obtaining a pore pressure field and a volume strain field true value by using a traditional numerical simulation; training a Fourier neural operator mechanical proxy model based on the data set, wherein the Fourier neural operator mechanical proxy model takes a pressure field, a Young's modulus field, a Poisson's ratio and a boundary load as input, and takes a volume strain field as output. Subsequently, a mixed iterative calculation framework based on a constant stress splitting is constructed: a fluid module adopts a finite volume method to solve a pressure field, a solid module calls the trained proxy model to quickly infer a strain field, and a two-way coupling iteration is realized through a porosity evolution model until convergence. The method replaces the traditional finite element solution by using the proxy model, and significantly reduces the solution cost of a large linear system under the premise of maintaining physical consistency and calculation accuracy.

Description

Technical Field

[0001] This invention relates to the field of porous media fluid-structure interaction simulation technology, and in particular to an iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators. Background Technology

[0002] The phenomenon of fluid flow and solid deformation interaction (i.e., fluid-structure interaction) in porous media is widely present in underground energy and environmental engineering fields (such as oil and gas extraction, geothermal energy development, carbon dioxide sequestration (CCS), and groundwater resource management). The fundamental physical principle of fluid-structure interaction is the strong nonlinear bidirectional coupling between fluid flow and solid skeleton deformation. Changes in fluid pore pressure lead to deformation of the rock skeleton, which in turn alters physical properties such as porosity and permeability, thus affecting the fluid flow pattern. Therefore, accurate and efficient numerical simulation of this complex coupling process is of great practical engineering significance for assessing reservoir stability, predicting production capacity, and ensuring engineering safety.

[0003] Currently, numerical simulation methods for solving fluid-structure interaction problems are mainly divided into two categories: fully coupled methods and iterative coupled methods. Fully coupled methods combine the fluid flow equations and solid mechanics equations into a huge linear system of equations to simultaneously solve for flow and mechanical parameters. While this method has the advantages of unconditional stability and good convergence, its system matrix is ​​enormous and its condition number is poor, placing extremely high demands on memory and computational resources. Therefore, it is generally not used in the development of complex practical models. Iterative coupled methods, on the other hand, solve the flow and mechanical subproblems independently, repeatedly transferring pressure and volumetric strain information between the two solvers until the convergence criterion is met. Compared to fully coupled methods, iterative coupled methods have advantages such as high modularity, ease of utilizing existing mature solvers, and high flexibility, making them more common in practical engineering applications.

[0004] Despite the widespread application of iterative coupling methods, they still face significant computational challenges when dealing with large-scale, highly heterogeneous, and complex geological models. First, in classic iterative coupling processes, the finite volume method (FVM) is typically used to solve the fluid equations to ensure local mass conservation, while the finite element method (FEM) is used to solve the mechanical equations. Since the mechanical subproblems often involve high-dimensional stiffness matrix inversion, this process accounts for a large portion of the computational overhead. Second, within each time step, multiple calls to the mechanical solver are often required to achieve coupling convergence. Furthermore, as the mesh size increases, the computation time of traditional finite element mechanical solvers increases dramatically and non-linearly, resulting in low overall simulation efficiency. This makes it difficult to meet the demands of "multi-query" tasks requiring thousands of simulations, such as inversion analysis, uncertainty quantification, or real-time optimization.

[0005] To address the aforementioned computational efficiency bottlenecks, data-driven deep learning methods have offered new approaches to accelerating scientific computing in recent years. However, purely data-driven models often lack physical interpretability and struggle to guarantee conservation laws. Therefore, combining efficient deep learning models with accurate traditional numerical methods to significantly reduce the most time-consuming mechanical solution costs while preserving the fluid mass conservation characteristics has become a critical technical challenge in the field of fluid-structure interaction numerical simulation. Summary of the Invention

[0006] This invention provides an accelerated iterative computational method for fluid-structure interaction in porous media based on Fourier neural operators. It aims to address the technical bottlenecks of high computational cost and low efficiency in traditional iterative coupling methods for solid mechanics solutions, particularly for large-scale heterogeneous strata simulation and massive "multi-query" tasks. By deeply integrating a data-driven surrogate model based on Fourier neural operators with the finite volume method that strictly guarantees mass conservation, it replaces the traditional numerical solver within a constant stress splitting framework. This significantly improves computational efficiency while ensuring the physical consistency of the fluid-structure interaction process. Ultimately, it forms a rapid simulation method with high accuracy and strong generalization ability, providing real-time and reliable technical support for efficient reservoir development and geomechanical risk assessment.

[0007] In a first aspect, the present invention provides an iterative accelerated computation method for porous media fluid-structure interaction based on Fourier neural operators, comprising: Based on random sampling, a multi-condition input sample containing heterogeneous parameters and boundary conditions is constructed. The time-varying pressure field is obtained by solving the fluid equation using the finite volume method. The time-varying pressure field and the corresponding physical property boundary parameters are input into the mechanical solver to calculate the volumetric strain field as the true value label. The data is then normalized and partitioned to construct a sample dataset of the fluid-structure interaction physical field of heterogeneous porous media. A neural network architecture containing multiple Fourier transform layers was constructed as a mechanical surrogate model. Pore pressure, Young's modulus, Poisson's ratio and boundary load were used as inputs and volumetric strain was used as output. The mechanical surrogate model was trained offline using the sample dataset until it converged. A fluid-structure interaction (FSI) hybrid iterative acceleration computation framework based on constant stress splitting is constructed. The FSI hybrid iterative acceleration computation framework includes decoupling the FSI system into independent fluid solution and solid inference modules using a constant stress splitting strategy. The fluid module uses the finite volume method to solve the pore pressure field, while the solid module uses a trained mechanical surrogate model to infer the volumetric strain field. A porosity evolution model is established as a two-way coupling link to feed back the skeleton deformation to the fluid equations, driving the fluid solution and solid inference modules to alternately iterate in a closed loop until the physical field satisfies the convergence criterion.

[0008] Furthermore, based on random sampling, a multi-condition input sample containing heterogeneous parameters and boundary conditions is constructed. The time-varying pressure field is obtained by solving the fluid equation using the finite volume method. This time-varying pressure field and the corresponding physical property boundary parameters are input into a mechanical solver to calculate the volumetric strain field as a truth label. Data normalization and set partitioning are then performed to construct a sample dataset of the fluid-structure interaction physics field of heterogeneous porous media, including: The heterogeneity of underground reservoirs is simulated using a Gaussian random field generation algorithm to generate a heterogeneous Young's modulus field. At the same time, the key physical parameters are sampled using the Latin hypercube sampling method to construct an input parameter set covering a variety of physical conditions. The key physical parameters include boundary load, initial pore pressure, Poisson's ratio, and bottom hole flowing pressure. Unsteady seepage simulation was performed on the input parameter set using the finite volume method. Key time nodes were selected based on the logarithmic distribution strategy. The transient pore pressure field distribution covering the entire process from the violent fluctuations in the early stage of mining to the pseudo-steady state in the later stage was calculated and extracted. The transient pore pressure field distribution is mapped to the pore elastic force of the rock skeleton. Based on the heterogeneous mechanical parameters and boundary loads, the data is input into a solid mechanics solver to obtain the corresponding volumetric strain field as the ground truth label for neural network training. The input and output data are normalized, and the preprocessed pairs of data are randomly divided into training and test sets.

[0009] Furthermore, the heterogeneity of the subsurface reservoir is simulated using a Gaussian random field generation algorithm to generate a heterogeneous Young's modulus field, including: A random noise matrix following a standard normal distribution is generated on a two-dimensional grid as a spatially uncorrelated basic random field. The random noise matrix is ​​then convolved using a two-dimensional Gaussian kernel function to obtain a smooth field with spatial correlation. The smoothed field is Z-score normalized to obtain the normalized field; Based on the mean and standard deviation of Young's modulus set for the current sample, the standardized field is mapped to the physical parameter space using a linear transformation, and physical upper and lower limits are applied to truncate the field, thus obtaining the Young's modulus field. Simultaneously, by combining the Latin hypercube sampling method, key physical parameters are sampled to construct an input parameter set covering various physical conditions, including: Multiple sets of normalized samples are generated in a six-dimensional unit hypercube, where the six dimensions are the initial pore pressure. Vertical effective stress Lateral pressure coefficient Bottom hole pressure drop ratio and the mean Young's modulus Compared to Poisson ;in: The first-dimensional sample is linearly mapped to a preset range as the initial pore pressure. ; Map the second-dimensional sample to vertical effective stress. The upper boundary load is calculated using the effective stress principle. : Based on vertical effective stress The right boundary load is calculated using the following formula. : In the formula, Indicates the effective horizontal stress; Based on bottom hole pressure drop ratio The bottom hole flowing pressure is calculated using the following formula. : Map the 5th dimension sample to the mean of Young's modulus. , used for Gaussian random field sampling; The 6th dimension sample is linearly mapped to a preset range as the Poisson ratio. .

[0010] Furthermore, the finite volume method is used to perform unsteady seepage simulation on the input parameter set. Based on a logarithmic distribution strategy, key time nodes are selected, and the transient pore pressure field distribution covering the entire process from the initial violent fluctuations in mining to the later quasi-steady state is calculated and extracted, including: Based on the law of conservation of mass and Darcy's law, an unsteady seepage simulation was conducted to construct the pressure diffusion equation for a single-phase microcompressible fluid in a porous medium: in, Pore ​​pressure, For rock porosity, For penetration rate, For fluid viscosity, The overall compression coefficient is... For source and sink items, For gradient operators, t For time; For spatial discretization, a central difference scheme is used to discretize the continuous computational domain into a structured grid system; for temporal discretization, a fully implicit backward first-order difference scheme is used. in, For time step, For fluid density, For the mesh volume, A set of units, The mobility coefficient, For unit and conductivity, For source and sink items, For unit The pressure value at the current moment, For unit The pressure value at the current moment, For unit At the current moment, the source and sink items, For unit The fluid mass at the current moment, For unit The fluid mass at the previous moment, Adjacent units and The interface shows the current mass flow rate parameters. Set simulation start time With end time and the preset number of sampling frames. Time series are generated using geometric series: in, t l For the first l A time slice l It is a time series.

[0011] Furthermore, the transient pore pressure field distribution is mapped to the pore elastic force of the rock skeleton. Based on the heterogeneous mechanical parameters and boundary loads, this is input into a solid mechanics solver for solving, obtaining the corresponding volumetric strain field as the ground truth label for neural network training, including: If the rock skeleton is a linear elastic isotropic medium and the deformation is less than a set threshold, then the mechanical behavior of the rock skeleton is determined to follow the following quasi-static equilibrium equation: in, For gradient operators, For the total stress tensor, , For the effective stress tensor, , It is a fourth-order elastic stiffness tensor. For strain tensor, This represents the tensor contraction operation. For Biot coefficient, Pore ​​pressure, For unit tensors, This is the force vector; strain tensor It is calculated using the following formula: in, Let be the displacement vector of all nodes in the field. This is the transpose of the gradient operator; Substituting the formula for calculating the strain tensor into the constitutive equation, and then combining it with the effective stress principle into the quasi-static equilibrium equation, the intermediate variable total stress is eliminated. With strain The displacement is derived. The fluid-structure interaction deformation control equation for porous media with the main variable as the primary variable is: Based on the heterogeneous Young's modulus field The global stiffness matrix, which reflects the heterogeneity of the reservoir space, is obtained through unit integration assembly. The nodal pore pressure field calculated at the current time step It is mapped to an equivalent fluid nodal force vector through shape function interpolation. The boundary load conditions are then transformed into external load vectors. Establish the following system of linear algebraic equations: By solving the above system of linear equations, the total field nodal displacement vector at the current time step is obtained. , T For matrix transpose, u x and u y For nodes in x direction and y Displacement components in the direction; Obtaining the total field nodal displacement vector Then, the volumetric strain was calculated using Cauchy's strain geometry equations: in, For spatial location ( x , y ),time t The volumetric strain below, for u x right x Partial derivatives of coordinates, for u y right y Partial derivatives of coordinates; For each logarithmic time step of each sample, the above steps are performed to obtain the time-varying volumetric strain field corresponding to the input pressure field and the inhomogeneous parameter field, and the time-varying volumetric strain field is used as the standard truth label.

[0012] Furthermore, a neural network architecture containing multiple Fourier transform layers is constructed as a mechanical surrogate model. Pore pressure, Young's modulus, Poisson's ratio, and boundary load are used as inputs, and volumetric strain is used as the output. The mechanical surrogate model is trained offline using the aforementioned sample dataset until convergence, including: The normalized two-dimensional pore pressure field and Young's modulus field are extracted as the basic channels; Poisson's ratio, overburden pressure and horizontal stress are extended into full-field feature maps as constant value channels using spatial broadcasting technology; deterministic position coding is introduced to generate normalized coordinate grid channels; all features are spliced ​​in the channel dimension to form a 7-channel composite input tensor. The input data is mapped to a high-dimensional feature space through a fully connected layer, and then enters a multi-layer stacked Fourier layer. Each layer adopts a dual-path architecture, and the two outputs are added together and then activated by a GELU layer. Finally, the high-dimensional features are regressed and mapped to the target volume strain field through a decoding projection layer. We use mean squared error as the loss function, use the Adam optimizer to update parameters, introduce a dynamic scheduling strategy, monitor the training set loss in real time, automatically decay the learning rate, and save the model parameters with the lowest loss as the mechanical surrogate model.

[0013] Furthermore, the space broadcasting technology includes: for each scalar Build a size of The full-field constant matrix As input; where Poisson's ratio, For the upper boundary load; For right boundary load; Deterministic position encoding is introduced to generate normalized coordinate grid channels; all features are concatenated along the channel dimension to form a 7-channel composite input tensor, including: Generate two first coordinate grid matrices with the same size as the physical grid. Second coordinate grid matrix Wherein, the first coordinate grid matrix The value increases linearly from 0 to 1 along the X-axis, while remaining constant along the Y-axis; this is the second coordinate grid matrix. The value increases linearly from 0 to 1 along the Y-axis, while the value remains unchanged along the X-axis. Based on the first coordinate grid matrix Second coordinate grid matrix A 7-channel composite input tensor is generated using the following formula: in, P For pore pressure, Concat is the feature splicing operation. E For Young's modulus field, This is the global constant matrix corresponding to Poisson's ratio. This is the global constant matrix corresponding to the overburden pressure or upper boundary load. This is the global constant matrix corresponding to the right boundary load.

[0014] Furthermore, the dual-path architecture includes a global path and a local path. The global path performs spectral convolution operations, and the local path is connected in parallel with a spatial convolutional layer to directly perform linear combination of features in the spatial domain. The dynamic scheduling strategy includes: if the training loss does not decrease significantly within a set number of consecutive training rounds, a learning rate decay mechanism is triggered to halve the current learning rate.

[0015] Furthermore, a fluid-structure interaction (FSI) hybrid iterative acceleration computation framework based on constant stress splitting is constructed. This framework includes a constant stress splitting strategy to decouple the FSI system into independent fluid solution and solid inference modules. The fluid module uses the finite volume method to solve the pore pressure field, while the solid module uses a trained mechanical surrogate model to infer the volumetric strain field. A porosity evolution model is established as a two-way coupling link, feeding back the skeleton deformation to the fluid equations to drive the fluid solution and solid inference modules in alternating iterative closed-loop cycles until the physical field satisfies the convergence criterion, including: The fluid control equations are solved using the finite volume method to obtain the updated pore pressure field for the current iteration step. In constructing the mass accumulation term of the fluid equations, a dynamic porosity evolution model is coupled to reflect the feedback of solid deformation to the fluid domain. The updated pore pressure field, Young's modulus field characterizing the heterogeneity of the medium, and boundary load conditions are input into a pre-trained mechanical surrogate model to infer the predicted value of the volumetric strain field at the current moment, and the predicted value of the volumetric strain field is denormalized to restore its physical dimensions. The relative norm error between the pore pressure field obtained in the current iteration step and the previous iteration step is calculated and compared with a preset convergence tolerance. If the error is greater than the convergence tolerance, the predicted volumetric strain field obtained in the current iteration step is used to update the rock porosity through the dynamic porosity evolution model, and the solution of the fluid control equation and subsequent inference process are re-executed to start the next coupling iteration. If the error is less than or equal to the convergence tolerance, the fluid-structure interaction calculation in the current time step is determined to have converged, the physical field update of the time step is completed, and the process is advanced to the next time step. This process is repeated until all simulations are completed.

[0016] Furthermore, the dynamic porosity evolution model is expressed as: in, It is the volumetric strain predicted by the trained mechanical surrogate model at the current time step. It is the volumetric strain at the initial time step. It is the porosity at the current time step. It is the porosity at the initial time step. It is the Biot coefficient. C r It is the overall compression coefficient. It is the pore pressure at the current time step. P 0 It is the pore pressure at the initial time step; The method for performing inverse normalization on the predicted volumetric strain field values ​​to restore their physical dimensions is as follows: in, The output of the trained mechanical surrogate model, This is the result of dimensionless feature reconstruction of the output of the trained mechanical surrogate model. The standard deviation of the normalized Young's modulus. The mean of the normalized Young's modulus. E ref The baseline Young's modulus value is typically taken as 1e9. E ( x ) represents the Young's modulus field as input.

[0017] The iterative accelerated computation method for fluid-structure interaction in porous media based on Fourier neural operators provided by this invention has at least the following beneficial effects: This invention deeply integrates the Fourier Neural Operator (FNO) and the Finite Volume Method (FVM) based on a constant stress splitting strategy. This significantly simplifies the computational complexity of mechanical solutions while ensuring fluid mass conservation through physical equation constraints, thus significantly improving the computational efficiency and physical fidelity of the hybrid model. Furthermore, this invention designs a multi-channel tensor input mechanism capable of handling heterogeneous Young's modulus and varying boundary loads, establishing a rapid proxy model for fluid-structure interaction with strong generalization capabilities. This enables real-time prediction of reservoir stress-sensitive effects under different geological conditions and injection-production scenarios, helping oilfields to formulate development plans more efficiently, shortening the production optimization decision-making cycle, and providing strong technical support for the long-term stable production and efficient development of complex reservoirs. Attached Figure Description

[0018] Figure 1A flowchart of an iterative accelerated computation method for porous media fluid-structure interaction based on Fourier neural operators is provided for an embodiment of the present invention. Figure 2 The iterative coupled computational framework diagram of the surrogate model replacing the mechanical solver provided in this embodiment of the invention; Figure 3 A schematic diagram of the model configuration and boundary conditions for a two-dimensional heterogeneous mining example provided in an embodiment of the present invention; Figure 4 A schematic diagram of representative samples of a heterogeneous Young's modulus field provided in an embodiment of the present invention; wherein, (a) sample 1; (b) sample 2; (c) sample 3; (d) sample 4; Figure 5 The loss function convergence curve of the proxy model training process provided in this embodiment of the invention; Figure 6 A time evolution comparison diagram of pore pressure and volumetric strain at the central monitoring point provided in an embodiment of the present invention; wherein, (a) is a comparison curve of pore pressure at the central monitoring point; and (b) is a comparison curve of volumetric strain at the central monitoring point. Figure 7 The histogram of the statistical distribution of the relative L2 error of the test set provided in the embodiments of the present invention; Figure 8 The pore pressure field and volumetric strain field comparison cloud map at the 100th time step provided in the embodiment of the present invention; wherein, (a) pressure field - true value; (b) pressure field - prediction; (c) pressure field absolute error; (d) volumetric strain - true value; (e) volumetric strain - prediction; (f) volumetric strain absolute error.

[0019] The accompanying drawings have illustrated specific embodiments of the invention, which will be described in more detail below. These drawings and descriptions are not intended to limit the scope of the invention in any way, but rather to illustrate the concept of the invention to those skilled in the art through reference to specific embodiments. Detailed Implementation

[0020] Exemplary embodiments will now be described in detail, examples of which are illustrated in the accompanying drawings. When the following description relates to the drawings, unless otherwise indicated, the same numerals in different drawings denote the same or similar elements. The embodiments described in the following exemplary embodiments do not represent all embodiments consistent with the present invention. Rather, they are merely examples of apparatuses and methods consistent with some aspects of the invention as detailed in the appended claims.

[0021] It should be noted that in the embodiments of the present invention, certain software, components, models and other existing solutions in the industry may be mentioned. These should be regarded as exemplary and are only intended to illustrate the feasibility of implementing the technical solution of the present invention. However, they do not mean that the inventor has used or necessarily used the solution.

[0022] The technical solution of the present invention and how the technical solution of the present invention solves the above-mentioned technical problems are described in detail below with specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments. The embodiments of the present invention will now be described with reference to the accompanying drawings.

[0023] This invention provides an accelerated iterative computational method for fluid-structure interaction (FSI) in porous media based on Fourier neural operators. First, a sample dataset containing various heterogeneous geological conditions is constructed. A mechanical surrogate model based on Fourier neural operators (FNO) is established and trained offline, enabling it to learn the complex nonlinear mapping relationship between input features (including pore pressure field, heterogeneous Young's modulus field, Poisson's ratio, and load conditions at the upper and right boundaries) and output features (volume strain field). Next, a hybrid iterative computational framework for FSI based on a constant stress splitting strategy is constructed. During the time-step-progressive online simulation, the finite volume method (FVM) is used to solve the fluid equations to strictly ensure local mass conservation, and the trained FNO model directly replaces the traditional solid mechanics solution steps. Within each iteration step, a porosity evolution model serves as a bidirectional coupling link, feeding back the skeleton deformation to the fluid equations, driving the two modules to iteratively close the loop until the convergence criterion is met, thereby achieving efficient and high-fidelity simulation of the FSI process in heterogeneous porous media. This method introduces a Fourier neural operator surrogate model to replace the traditional finite element method and directly participates in the fluid-structure interaction iterative calculation process. While ensuring physical consistency and calculation accuracy, it significantly reduces the computational cost of solving large linear systems and solves the computational bottleneck problem of traditional fluid-structure interaction simulation in long-cycle, multi-query tasks. It is of great significance for efficient reservoir development and geomechanical risk assessment.

[0024] Specifically, such as Figure 1 As shown, the porous medium fluid-structure interaction iterative acceleration calculation method based on Fourier neural operators can be implemented through the following steps S10 to S30.

[0025] S10: Construct a multi-condition input sample containing heterogeneous parameters and boundary conditions based on random sampling, use the finite volume method to solve the fluid equation to obtain the time-varying pressure field, input the time-varying pressure field and the corresponding physical property boundary parameters into the mechanical solver to calculate the volumetric strain field as the true value label, and complete the data normalization and set partitioning to construct a sample dataset of the fluid-structure interaction physical field of heterogeneous porous media.

[0026] In some embodiments, step S10 is specifically implemented through the following steps S101-S104.

[0027] S101: The heterogeneity of underground reservoirs is simulated using the Gaussian random field generation algorithm to generate a heterogeneous Young's modulus field. At the same time, the Latin hypercube sampling (LHS) method is combined to sample key physical parameters such as boundary load, initial pore pressure, Poisson's ratio, and bottom hole flowing pressure to construct an input parameter set covering a variety of physical conditions.

[0028] In some embodiments, the Gaussian random field generation algorithm in step S101 is as follows: To simulate the spatial heterogeneity and continuity of the mechanical properties of underground reservoir rocks, this invention employs a random field generation algorithm based on Gaussian kernel convolution. This method has higher computational efficiency when processing high-resolution grids. The specific operation steps are as follows: First of all, Generate data on a two-dimensional grid that follows a standard normal distribution. random noise matrix As a spatially uncorrelated fundamental random field, then using a two-dimensional Gaussian kernel function For random noise matrix Perform convolution operations to obtain a smooth field with spatial correlation. : (1) in, For grid coordinate information, The standard deviation of the Gaussian kernel function is used to control the spatial correlation length, i.e., the spatial scale of the heterogeneous patch.

[0029] For smooth fields Z-score standardization is performed to make its mean 0 and standard deviation 1, resulting in the standardized field. : (2) Finally, based on the mean Young's modulus set for the current sample... and standard deviation The normalized field is mapped to the physical parameter space using a linear transformation, and physical upper and lower limits are applied for truncation. Ultimately, Young's modulus field was obtained. : (3) Clip is a truncation operation.

[0030] In some embodiments, the physical parameter sampling method in step S101 is as follows: This embodiment simulates the extraction process of an underground oil reservoir. Utilizing geometric and physical symmetry, it simulates only one quadrant of the reservoir (usually the upper right quadrant), thus reducing computational load. To ensure that the constructed dataset has both broad parameter coverage and strict adherence to rock mechanics and seepage physics laws (such as constant positive effective stress and no negative pressure), this embodiment constructs a Latin hypercube sampling space with six independent dimensions and calculates the final input parameters through physical constraint mapping. First, using LHS... D =6-dimensional unit hypercube [0,1] 6 Generated in Group Normalized Samples To ensure uniform distribution and space-filling of the sample within the parameter space, the six dimensions correspond to: initial pore pressure, etc. Vertical effective stress Lateral pressure coefficient Bottom hole pressure drop ratio and the mean Young's modulus Compared to Poisson The processing of the 6-dimensional parameters is as follows: Initial pore pressure Linearly map the first-dimensional sample to a preset range; upper boundary load To avoid the non-physical phenomenon of "negative effective stress" caused by random sampling, this method does not directly sample the load, but first maps the second-dimensional sample to vertical effective stress. Then, the upper boundary load is calculated back using the effective stress principle: (4) Right boundary load Introducing the lateral pressure coefficient from rock mechanics (The third-dimensional sample) is used to characterize the relationship between horizontal and vertical stress. First, the effective horizontal stress is calculated. Then calculate the right boundary load. This step ensures that the generated stress field conforms to the geological structure, covering a variety of conditions from extensional to compressional structures. (5) Bottom hole flowing pressure Utilizing the bottom hole pressure drop ratio (4th dimension sample) Determine the production well pressure to ensure a reasonable production pressure differential and prevent negative pressure: (6) Mean Young's modulus Map the 5th dimension sample to the mean of Young's modulus. , used for Gaussian random field sampling in step S101; Poisson's ratio Linearly map the 6th dimension sample to a preset range.

[0031] S102: The unsteady seepage simulation of the above-generated sample (input parameter set) is performed using FVM. Based on the logarithmic distribution strategy, key time nodes are selected, and the transient pore pressure field distribution covering the entire process from the initial violent fluctuations to the later pseudo-steady state is calculated and extracted.

[0032] In some embodiments, in step S102, the unsteady-state seepage simulation is based on the law of conservation of mass and Darcy's law to construct the pressure diffusion equation for a single-phase microcompressible fluid in a porous medium: (7) in, Pore ​​pressure, For rock porosity, For penetration rate, For fluid viscosity, The overall compression coefficient is... For injection and production well source confluence, For gradient operators, t For time.

[0033] For spatial discretization, a central difference scheme is used to discretize the continuous computational domain into a structured grid system; for temporal discretization, a fully implicit backward first-order difference scheme is used. (8) in, For time step, For fluid density, For mesh volume, subscript Indicates adjacent units and A suitable average value of physical quantities at the interface. A set of elements, with superscript For the previous moment , n +1 represents the current time. , The mobility coefficient, For unit and conductivity, For source and sink items, For unit At the present moment Pressure value, For unit At the present moment Pressure value, For unit At the present moment The source and sink items, For unit At the present moment fluid mass, For unit In the previous moment fluid mass, Adjacent units and The interface shows the mass flow rate parameters at the current moment. The fully implicit scheme has unconditional stability, ensuring that the solver can still converge stably without numerical oscillations when using a large time step with a logarithmic distribution for simulation.

[0034] In some embodiments, the logarithmic distribution sampling strategy for key time nodes in step S102 is as follows: Considering the nonlinear physical characteristics of the pressure diffusion process in porous media, which involves drastic initial changes followed by a gradual plateauing, this invention employs a logarithmic distribution strategy to select key time nodes. First, the simulation start time is set. With end time and the preset number of sampling frames. Time series are generated using geometric series: in, t l For the first l A time slice l It is a time series.

[0035] The logarithmic distribution sampling strategy at this key time point enables adaptive adjustment of the time step. In the initial stage of well operation (transient phase), the time step is extremely small, and high-frequency sampling accurately captures the steep gradient changes in wellbore pressure, providing rich high-frequency features for the neural network. In the later stage of production (quasi-steady-state phase), the time step increases exponentially to avoid data redundancy and cover the long-cycle pressure dissipation process. Finally, this... The pressure field snapshots corresponding to each key time step form a time-varying pressure dataset with high physical information density.

[0036] S103: The time-varying pore pressure field obtained in step S102 (i.e., the transient pore pressure field distribution covering the entire process from the violent fluctuations in the early stage of mining to the pseudo-steady state in the later stage) is mapped to the pore elastic force of the rock skeleton. Combined with the heterogeneous mechanical parameters and boundary loads in step S101, it is input into the solid mechanics solver for solving, and the corresponding volumetric strain field is obtained as the ground truth label for neural network training.

[0037] In some embodiments, in step S103, the solid mechanics solution is based on Biot's porous elastic theory: assuming the rock skeleton is a linear elastic isotropic medium with minimal deformation, its mechanical behavior follows quasi-static equilibrium equations. (10) Introducing the effective stress principle ,in For the total stress tensor, For constitutive equations, For the effective stress tensor, It is a fourth-order elastic stiffness tensor. For strain tensor, This represents the tensor contraction operation. For Biot coefficient, The pore pressure calculated in step 1.2, For unit tensors, Let be the force vector. The strain tensor can be obtained from the following equation: (11) in, Let be the displacement vector of all nodes in the field. This is the transpose of the gradient operator. Substituting the geometric equation (11) into the constitutive equation, and then combining it with the effective stress principle into the quasi-static equilibrium equation (10), the total stress is eliminated by eliminating intermediate variables. With strain It can be deduced that the displacement The fluid-structure interaction deformation control equation for porous media with the main variable as the primary variable is: (12) The first term of the above equation (12) represents the elastic internal force of the rock skeleton, the second term represents the pore elastic volume force transformed from the fluid pore pressure gradient, this term realizes the unidirectional coupling drive of the fluid field to the solid field, and the third term is the body force term represented by gravity.

[0038] In some embodiments, in step S103, the solution method is the finite element method. First, based on the heterogeneous Young's modulus field generated in step S101... The global stiffness matrix, which reflects the heterogeneity of the reservoir space, is obtained through unit integration assembly. Secondly, the nodal pore pressure field calculated in step S102 at the current time step is... It is mapped to an equivalent fluid nodal force vector through shape function interpolation. And the boundary load conditions (upper boundary load) sampled in step S101 and right boundary load Transform into an external load vector Finally, the following system of linear algebraic equations is established: (13) By solving the above system of linear equations, the total field nodal displacement vector at the current time step can be obtained. , T For matrix transpose, u x and u y For nodes in x direction and y Displacement in the direction. Obtaining the total field nodal displacement vector. Then, based on the small deformation assumption, the volumetric strain was calculated using the Cauchy strain geometry equation: (14) in, For volumetric strain, For spatial location ( x , y ),time t The volumetric strain below, for u x right x Partial derivatives of coordinates, for u y right y Partial derivatives of coordinates.

[0039] For each logarithmic time step of each sample, step S103 is executed to obtain a time-varying volumetric strain field that strictly corresponds to the input pressure field and the inhomogeneous parameter field. This time-varying volumetric strain field is finally extracted and stored as the standard ground truth label for subsequent training of the Fourier neural operator (FNO).

[0040] S104: To eliminate the influence of differences in the dimensions of different physical quantities on network training, the input data (pressure field, Young's modulus field, load conditions, Poisson's ratio) and output data (volume strain field) are normalized, and the preprocessed paired data are randomly divided into training set and test set.

[0041] In some embodiments, the normalization method used in step S104 is as follows: the physical quantities involved in the neural network input and output data span a huge difference in magnitude (e.g., Young's modulus is 10). 10 The level is 10, while the volumetric strain is 10. -4 If the input is directly fed into a neural network, the tiny fluctuations in the input will be masked by the huge cardinality, and the gradient of the output will vanish because the value is too small, causing the model to fail to converge. Therefore, it is necessary to eliminate the huge order-of-magnitude gap between the physical quantities in a fluid-structure interaction system. This invention uses a Gaussian normalization strategy based on statistical properties to preprocess the input and output data.

[0042] First, a logarithmic transformation is performed on all input variables. Then, based on the statistical data of the training set, Z-score standardization is performed on the transformed logarithmic features, uniformly mapping all input features to a dimensionless distribution with zero mean and unit variance. (15) in, For standardized data, The data is after logarithmic transformation. The mean of the data after logarithmic transformation. This represents the variance of the data after logarithmic transformation. For the output data, volumetric strain field... Because its value is extremely small and easily affected by numerical noise, a reference Young's modulus is introduced. A physical space transformation is performed to map the dimensionless volumetric strain field into an equivalent pseudo-stress field. : (16) Subsequently Z-score standardization is performed, and to prevent data leakage, all the statistical parameters required for the above normalization ( The parameters are calculated solely based on the training set and then frozen for use on the test set and in future inference processes.

[0043] S20: Construct a neural network architecture containing multiple Fourier transform layers as a mechanical surrogate model, using pore pressure, Young's modulus, Poisson's ratio and boundary load as inputs and volumetric strain as output, and train the mechanical surrogate model offline using the sample dataset until convergence.

[0044] In some embodiments, step S20 is specifically implemented through the following steps S201-S203.

[0045] S201: The normalized two-dimensional pore pressure field and Young's modulus field are extracted as the basic channels; using spatial broadcasting technology, scalar parameters such as Poisson's ratio, overburden pressure, and horizontal stress are expanded into a full-field feature map as constant-value channels; simultaneously, deterministic position encoding is introduced to generate normalized coordinate grid channels. All the above features are concatenated along the channel dimensions to form a 7-channel composite input tensor containing physical fields, boundary conditions, and spatial coordinates.

[0046] In some embodiments, in step S201, the spatial broadcasting technology specifically works as follows: Since the Fourier neural operator (FNO) typically processes spatially distributed data, it cannot directly receive scalar inputs, and for each sample's time-varying pore pressure field... With static inhomogeneous Young's modulus field All A 2D matrix. Therefore, for each scalar Build a size of The full-field constant matrix As input, it is necessary that the value of an element at any position in the matrix is ​​equal to the normalized value of the scalar.

[0047] In some embodiments, in step S201, the deterministic location encoding is specifically as follows: Considering that the rock mechanical response has a significant spatial location dependence, and that a simple physical field input lacks absolute coordinate information, deterministic location encoding is introduced to enhance the model's ability to perceive spatial distribution and non-periodic boundaries. Specifically, this involves generating two [grids] with the same dimensions as the physical grid. Consistent first coordinate grid matrix Second coordinate grid matrix First coordinate grid matrix The value increases linearly from 0 to 1 along the X-axis, while the value remains unchanged along the Y-axis. Second coordinate grid matrix The values ​​increase linearly from 0 to 1 along the Y-axis and remain unchanged along the X-axis. These two coordinate matrices constitute two additional geometric feature channels. By introducing coordinate channels, the translation invariance of traditional convolutional networks is explicitly broken, forcing the model to learn the differentiated mechanical response patterns at different spatial locations.

[0048] In some embodiments, in step S201, the 7-channel composite input tensor is specifically: (17) in, P For pore pressure, Concat is the feature splicing operation. E For Young's modulus field, This is the global constant matrix corresponding to Poisson's ratio. This is the global constant matrix corresponding to the overburden pressure or upper boundary load. This is the global constant matrix corresponding to the right boundary load.

[0049] The generated composite input tensor The dimension is (in (For batch size), this tensor fully contains four key pieces of information: fluid state, medium properties, boundary conditions, and geometric position, providing complete physical feature inputs for high-precision inference of subsequent Fourier neural operators.

[0050] S202: The input data is first mapped to a high-dimensional feature space through a fully connected layer; then it enters a multi-layer stacked Fourier layer, each layer adopting a dual-path architecture: one path extracts global low-frequency deformation features through two-dimensional fast Fourier transform (FFT) and modal truncation, and the other path preserves local details through spatial convolution. The two are added together and then activated by the GELU layer; finally, the high-dimensional features are regressed and mapped to the target volume strain field through the decoding projection layer.

[0051] In some embodiments, in step S202, the multi-layered stacked Fourier layers specifically refer to: in this invention, the number of Fourier layers is 4, and each layer adopts a dual-path architecture.

[0052] In some embodiments, in step S202, the dual-path architecture specifically comprises a global path and a local path. The global path performs spectral convolution operations, responsible for capturing long-range dependencies and overall trends; the local path is connected in parallel with a kernel of size [missing information]. The spatial convolutional layers directly combine features linearly in the spatial domain, responsible for transmitting local node information and improving gradient flow. The outputs of the two pathways are added element-wise, fusing global trends and local details, and then passed through the GELU nonlinear activation function. This design ensures that the model can simulate both macroscopic rock skeleton deformation and finely characterize the local response of heterogeneous regions.

[0053] In some embodiments, in step S202, the modal truncation and complex linear transformation in the frequency domain specifically involves: in the global path, firstly, a two-dimensional fast Fourier transform is performed to transform the feature map in the real domain to the frequency domain. Utilizing the conjugate symmetry of the real input, only the non-redundant half-spectrum is calculated and stored, thereby reducing the computational load; secondly, low-frequency modal truncation is performed. Based on the physical prior that rock deformation is usually dominated by overall continuous deformation with fewer local abrupt changes, a truncation wavenumber is set. In the frequency domain, only the low-frequency modes in the upper left corner of the spectrum (i.e., large-scale features with concentrated energy) are retained, while the coefficients of other modes representing high-frequency noise are forced to zero. This is equivalent to applying a learnable low-pass filter to the physical field. Then, complex weighted multiplication is performed to linearly transform the retained low-frequency modes, thereby achieving the mixing of global information and feature extraction. Finally, the processed frequency domain features are restored back to the spatial domain through inverse Fourier transform.

[0054] In some embodiments, in step S202, the decoding projection specifically involves processing the high-dimensional deep features extracted after passing through L layers (4 layers in this embodiment) of stacked Fourier layers, and then introducing these features into the decoding projection layer. This layer consists of two cascaded pointwise fully connected layers with a GELU activation function in between, the function of which is to channel the high-dimensional latent features. The dimensions are gradually reduced and compressed, and finally mapped back to the target output channel. (i.e., volumetric strain field).

[0055] S203: The mean squared error is used as the loss function, the Adam optimizer is used to update the parameters, and the ReduceLROnPlateau dynamic scheduling strategy is introduced to monitor the training set loss in real time, automatically decay the learning rate to promote model convergence, and finally save the model parameters with the lowest loss as the trained mechanical surrogate model.

[0056] In some embodiments, in step S203, the loss function is specifically: assuming a batch contains... For the nth sample, for the nth For each sample, the normalized prediction tensor output by the model is: The corresponding normalized truth label is The loss function is then defined as: (18) in, This represents the Euclidean norm. This loss function provides a smooth and continuous gradient guide for updating network weights during backpropagation, driving the model to approximate the true mechanical response in the normalized feature space.

[0057] In some embodiments, in step S203, the ReduceLROnPlateau dynamic scheduling strategy specifically involves: if the training loss does not decrease significantly within 50 consecutive training rounds (i.e., entering a convergence plateau), then a learning rate decay mechanism is triggered, halving the current learning rate (e.g., from...). Down to This coarse-to-fine search strategy effectively prevents the model from oscillating near the minimum due to excessively large step sizes, thus promoting further convergence of the model.

[0058] S30: Construct a fluid-structure interaction hybrid iterative acceleration computation framework based on constant stress splitting. The fluid-structure interaction hybrid iterative acceleration computation framework includes decoupling the fluid-structure interaction system into independent fluid solution and solid inference modules using a constant stress splitting strategy. The fluid module uses the finite volume method to solve the pore pressure field, and the solid module uses a trained mechanical surrogate model to infer the volumetric strain field. A porosity evolution model is established as a two-way coupling link to feed back the skeleton deformation to the fluid equation, driving the fluid solution and solid inference modules to alternately iterate in a closed loop until the physical field satisfies the convergence criterion.

[0059] In some embodiments, step S30 is specifically implemented through the following steps S301-S303.

[0060] S301: The implicit FVM discrete scheme constructed in step S102 is used. However, when constructing the mass accumulation term, porosity is no longer regarded as a simple function of pressure. Instead, a dynamic porosity evolution model is introduced, and the updated pore pressure field of the current iteration step is obtained by solving it.

[0061] In some embodiments, the dynamic porosity evolution model in step S301 is as follows: (19) in, It is the volumetric strain field predicted by a trained mechanical surrogate model. This represents the pore pressure at the current time step, while variables marked with a superscript of 0 represent values ​​from the initial time step. This equation explicitly incorporates the pore volume change caused by skeleton deformation into the fluid equation, achieving fluid-structure interaction.

[0062] S302: After obtaining the updated pore pressure field, input the pressure field, the heterogeneous Young's modulus field and the boundary load into the trained mechanical surrogate model FNO to quickly infer the volumetric strain field at the current moment, replacing the traditional finite element stiffness matrix solution process.

[0063] In some embodiments, the input data in step S302 must undergo strict consistency normalization. Using the mean and variance of the training set obtained in step S10, the pore pressure field and Young's modulus field of the current iteration step are logarithmically transformed and Z-score standardized. This ensures that the statistical distribution of the online inference data is consistent with that of the offline training domain, preventing inference failures caused by differences in numerical magnitude.

[0064] In some embodiments, the FNO output in step S302 needs to be denormalized. First, the dimensionless features of the network output are restored to the statistical space before normalization, and then the physical dimensions are restored by combining the current Young's modulus field. (20) (twenty one) in, The output of the trained mechanical surrogate model, This is the result of dimensionless feature reconstruction of the output of the trained mechanical surrogate model. The mean of the normalized Young's modulus. E ref The baseline Young's modulus value is typically taken as 1e9. E ( x ) represents the Young's modulus field as input.

[0065] S303: Calculate the relative norm error of the pore pressure field between the current iteration step and the previous iteration step. If the error is greater than the preset convergence tolerance, the rock porosity is corrected using the latest volumetric strain field predicted by FNO, and the process returns to step S301 to continue the next coupled iteration. If the error meets the convergence requirements, the current time step is determined to be completed, the physical state of the entire field (including pressure, strain and porosity) is updated, and the process is advanced to the next time step until the entire process simulation is completed.

[0066] In some embodiments, the convergence determination criterion in step S303 is as follows: calculating the current iteration step. k +1 and the previous iteration step k pore pressure p k+1 and p k The absolute maximum value of the difference is normalized using the characteristic amplitude of the current pressure field: (twenty two) in, This represents the component with the largest absolute value in the vector. To preset the convergence tolerance (e.g.) ), To prevent extremely small stable terms with a denominator of zero, this criterion ensures that pressure fluctuations at any grid node within the computational domain have decayed to a negligible level, guaranteeing the pointwise convergence of the numerical solution.

[0067] In some embodiments, such as Figure 2 The diagram shows the iterative coupled computational framework of the surrogate model replacing the mechanical solver. The specific data flow of this iterative coupled computational framework is as follows: In the iteration process of each time step, it starts with the flow module. The fluid solution adopts the traditional finite volume method. The solver receives the porosity and other parameters determined in the previous iteration step, solves the discretized pressure diffusion equation, and obtains the pore pressure field distribution at the current moment; subsequently, this pressure field... P Compared with static geomechanical parameters (heterogeneous Young's modulus) E Poisson's ratio v and boundary loads ( q top , q right Channel concatenation is performed to dynamically construct the multi-channel feature tensor of the input FNO. Then, a forward mapping is performed using the trained FNO model to directly output the volumetric strain field under the current working condition. (output space) This process achieves instantaneous response to mechanical problems; finally, it updates reservoir porosity based on volumetric strain, feeding back the effects of mechanical deformation to the fluid system; during iteration, if the residual is less than the preset tolerance, the current time step is considered converged, and the process proceeds to the next time step; otherwise, the updated physical property parameters are used to return to the fluid solution module. The lower right corner of the diagram shows a network schematic, in which... Indicates Fourier transform, It is a complex matrix. Indicates the inverse Fourier transform. W It is a linear matrix. This is the activation function.

[0068] In some embodiments, after step S30, the method further includes step S40, which is used to verify the accuracy of the hybrid iterative acceleration algorithm and evaluate its generalization ability under multiple working conditions. Based on the results of traditional full numerical simulation, the relative error between the pressure field and the volumetric strain field obtained by the hybrid iterative calculation is calculated to evaluate the algorithm's accuracy. The trained mechanical surrogate model is applied to non-homogeneous parameter distributions and varying boundary load conditions that were not trained to verify the model's extrapolation and generalization ability. Simultaneously, the number of iterations and computation time are statistically analyzed to quantitatively evaluate the acceleration benefits of the proposed method.

[0069] The specific process of step S40 is as follows: S401: For the physical conditions in the test set, a full numerical coupled simulation is performed using the traditional finite element method (FEM) mechanical solver and the FVM fluid solver. Under the same time step and mesh generation conditions, the pore pressure field and volumetric strain field of the whole cycle are calculated and used as the "true value" to evaluate the accuracy of the hybrid iterative algorithm. Among them, the traditional finite element method (FEM) mechanical solver and the FVM fluid solver are the methods described in step S10, and the coupling strategy between the fluid and solid fields is the method described in step S30.

[0070] S402: Compare the predicted field obtained by the hybrid iteration calculation in step S30 with the reference field in step S401 point by point, calculate the relative error norm and point-by-point absolute error distribution of the pore pressure field and the volumetric strain field at different times, and verify whether the physical fidelity of the fluid-structure interaction calculation can still be maintained after replacing the mechanical solver with FNO.

[0071] The relative error norm is as follows: for any time... t The physical field vectors predicted by the hybrid algorithm are denoted as... The full numerical reference field is denoted as relative error norm The calculation formula is: (twenty three) This index can comprehensively reflect the average deviation of the pore pressure field and the volumetric strain field throughout the entire computational domain, effectively avoiding evaluation bias caused by local extreme points.

[0072] Among them, point-by-point absolute error E abs ( x , y The details are as follows: (twenty four) in, For in position ( x , y The physical field vector predicted by the hybrid algorithm at position ) For in position ( x , y The full numerical reference field at ().

[0073] This error can verify whether the hybrid algorithm has non-physical numerical artifacts or smoothing effects, ensuring that the FNO surrogate model not only fits the overall trend, but also accurately captures the local high gradient features in fluid-structure interaction.

[0074] S403: Select entirely new samples outside the training set distribution for testing, and evaluate the extrapolation and generalization ability of the mechanical surrogate model by comparing the prediction errors under different working conditions; S404: Record and compare the total computation time of traditional full numerical simulation and FVM-FNO hybrid iterative algorithm when processing the same physical time span, count the average number of nonlinear iterations of the two in each time step, calculate the speedup ratio of the hybrid algorithm compared with the traditional method, and thus quantitatively evaluate the efficiency advantage of the proposed method in engineering applications while ensuring computational accuracy.

[0075] To demonstrate the feasibility and superiority of the present invention, the following embodiments are provided.

[0076] This embodiment uses a two-dimensional heterogeneous mining scenario as a benchmark example to evaluate the performance of a surrogate model replacing the mechanical solver to accelerate fluid-structure interaction iterations. The model configuration is as follows: Figure 3 As shown: The simulation area is 150m × 150m in size, discretized into a 30 × 30 grid; a constant-pressure production well is set in the lower left corner of the area, with a production pressure of [value missing]. p wf A constant compressive load is applied to the top and right boundaries, while rolling supports (normal displacement constraints) are applied to the left and bottom boundaries. The simulation process evolves over 100 time steps, with each step lasting 50 seconds. Figure 4 Several representative implementations of nonhomogeneous Young's modulus fields in the test set are presented.

[0077] Figure 5The loss curves during training are shown, indicating that the model converges and stabilizes during training. To further verify the coupling performance, Figure 6 The evolution of pore pressure and volumetric strain at the monitoring point at the center of the region over time was compared. The solid and dashed lines in the figure represent the iterative coupled calculation results using the standard finite element method (FEM) solver and the FNO surrogate model as the mechanical solver, respectively. After 100 time steps, the curves show a high degree of agreement, confirming that the FNO surrogate model can effectively replace the traditional mechanical solver, thus verifying the reliability of the proposed framework.

[0078] Figure 7 The relative L2 error distribution across the entire test dataset was statistically analyzed. Histograms show that the relative error for the vast majority of samples is below 1%. Specifically, a mean error of 0.3% reflects the model's high mean accuracy, confirming the surrogate model's effective generalization ability to unseen data. To further evaluate the field distribution accuracy under coupled conditions, Figure 8 This paper presents a detailed comparison of the pore pressure field and volumetric strain field obtained by standard FEM coupling and surrogate-assisted coupling at the 100th time step (using...). Figure 4 (The first non-homogeneous Young's modulus field sample in the model). The results show that even after 100 time-step iterative coupling, the calculation results based on the surrogate model still agree very well with the finite element solution. This result proves that even with repeated calls to the surrogate model, the solution remains stable and accurate, without significant cumulative error.

[0079] In terms of computational efficiency, the proposed method demonstrates significant advantages. The total simulation runtime is reduced from 63.12 seconds for the standard FEM to 37.76 seconds for the surrogate model, a reduction of 40.2%. Notably, the surrogate model accelerates specific mechanics solution stages by approximately 24.7 times (single-step time reduced from 0.0857 seconds to 0.0035 seconds). These results validate the effectiveness of the proposed framework as an efficient alternative to long-term coupled simulations. Furthermore, given the rapidly increasing computational burden of solving linear systems with the standard FEM, this speedup advantage is expected to be even more pronounced in larger-scale or more complex models.

[0080] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for accelerated iterative computation of fluid-structure interaction in porous media based on Fourier neural operators, characterized in that, include: Based on random sampling, a multi-condition input sample containing heterogeneous parameters and boundary conditions is constructed. The time-varying pressure field is obtained by solving the fluid equation using the finite volume method. The time-varying pressure field and the corresponding physical property boundary parameters are input into the mechanical solver to calculate the volumetric strain field as the true value label. The data is then normalized and partitioned to construct a sample dataset of the fluid-structure interaction physical field of heterogeneous porous media. A neural network architecture containing multiple Fourier transform layers was constructed as a mechanical surrogate model. Pore pressure, Young's modulus, Poisson's ratio and boundary load were used as inputs and volumetric strain was used as output. The mechanical surrogate model was trained offline using the sample dataset until it converged. A fluid-structure interaction (FSI) hybrid iterative acceleration computation framework based on constant stress splitting is constructed. This framework includes a constant stress splitting strategy to decouple the FSI system into independent fluid solution and solid inference modules. The fluid module uses the finite volume method to solve the pore pressure field, while the solid module uses a trained mechanical surrogate model to infer the volumetric strain field. A porosity evolution model is established as a two-way coupling link, feeding back the skeleton deformation to the fluid equations to drive the fluid solution and solid inference modules in alternating iterative closed-loop cycles until the physical field satisfies the convergence criterion. A multi-condition input sample containing heterogeneous parameters and boundary conditions is constructed based on random sampling. The time-varying pressure field is obtained by solving the fluid equation using the finite volume method. This time-varying pressure field and the corresponding physical property boundary parameters are input into a mechanical solver to calculate the volumetric strain field as a truth label. Data normalization and set partitioning are then performed to construct a sample dataset of the fluid-structure interaction physics field in heterogeneous porous media, including: The heterogeneity of underground reservoirs is simulated using a Gaussian random field generation algorithm to generate a heterogeneous Young's modulus field. At the same time, the key physical parameters are sampled using the Latin hypercube sampling method to construct an input parameter set covering a variety of physical conditions. The key physical parameters include boundary load, initial pore pressure, Poisson's ratio, and bottom hole flowing pressure. Unsteady seepage simulation was performed on the input parameter set using the finite volume method. Key time nodes were selected based on the logarithmic distribution strategy. The transient pore pressure field distribution covering the entire process from the violent fluctuations in the early stage of mining to the pseudo-steady state in the later stage was calculated and extracted. The transient pore pressure field distribution is mapped to the pore elastic force of the rock skeleton. Based on the heterogeneous mechanical parameters and boundary loads, the data is input into a solid mechanics solver to obtain the corresponding volumetric strain field as the ground truth label for neural network training. The input and output data are normalized, and the preprocessed pairs of data are randomly divided into training and test sets.

2. The iterative accelerated computation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 1, characterized in that, The heterogeneity of subsurface reservoirs is simulated using a Gaussian random field generation algorithm, generating a heterogeneous Young's modulus field, including: A random noise matrix following a standard normal distribution is generated on a two-dimensional grid as a spatially uncorrelated basic random field. The random noise matrix is ​​then convolved using a two-dimensional Gaussian kernel function to obtain a smooth field with spatial correlation. The smoothed field is Z-score normalized to obtain the normalized field; Based on the mean and standard deviation of Young's modulus set for the current sample, the standardized field is mapped to the physical parameter space using a linear transformation, and physical upper and lower limits are applied to truncate the field, thus obtaining the Young's modulus field. Simultaneously, by combining the Latin hypercube sampling method, key physical parameters are sampled to construct an input parameter set covering various physical conditions, including: Multiple sets of normalized samples are generated in a six-dimensional unit hypercube, where the six dimensions are the initial pore pressure. Vertical effective stress Lateral pressure coefficient Bottom hole pressure drop ratio and the mean Young's modulus Compared to Poisson ;in: The first-dimensional sample is linearly mapped to a preset range as the initial pore pressure. ; Map the second-dimensional sample to vertical effective stress. The upper boundary load is calculated using the effective stress principle. : Based on vertical effective stress The right boundary load is calculated using the following formula. : In the formula, Indicates the effective horizontal stress; Based on bottom hole pressure drop ratio The bottom hole flowing pressure is calculated using the following formula. : Map the 5th dimension sample to the mean of Young's modulus. , used for Gaussian random field sampling; The 6th dimension sample is linearly mapped to a preset range as the Poisson ratio. .

3. The iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 1, characterized in that, Unsteady-state seepage simulation was performed on the input parameter set using the finite volume method. Based on a logarithmic distribution strategy, key time points were selected, and the transient pore pressure field distribution covering the entire process from the initial violent fluctuations in mining to the later quasi-steady state was calculated and extracted, including: Based on the law of conservation of mass and Darcy's law, an unsteady seepage simulation was conducted to construct the pressure diffusion equation for a single-phase microcompressible fluid in a porous medium: in, Pore ​​pressure, For rock porosity, For penetration rate, For fluid viscosity, The overall compression coefficient is... For source and sink items, For gradient operators, t For time; For spatial discretization, a central difference scheme is used to discretize the continuous computational domain into a structured grid system; for temporal discretization, a fully implicit backward first-order difference scheme is used. in, For time step, For fluid density, For the mesh volume, A set of units, The mobility coefficient, For unit and conductivity, For source and sink items, For unit The pressure value at the current moment, For unit The pressure value at the current moment, For unit At the current moment, the source and sink items, For unit The fluid mass at the current moment, For unit The fluid mass at the previous moment, Adjacent units and The interface shows the current mass flow rate parameters. Set simulation start time With end time and the preset number of sampling frames. Time series are generated using geometric series: in, t l For the first l A time slice l It is a time series.

4. The iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 1, characterized in that, The transient pore pressure field distribution is mapped to the pore elastic force of the rock skeleton. Based on the heterogeneous mechanical parameters and boundary loads, the data is input into a solid mechanics solver for solving, and the corresponding volumetric strain field is obtained as the ground truth label for neural network training, including: If the rock skeleton is a linear elastic isotropic medium and the deformation is less than a set threshold, then the mechanical behavior of the rock skeleton is determined to follow the following quasi-static equilibrium equation: in, For gradient operators, For the total stress tensor, , For the effective stress tensor, , It is a fourth-order elastic stiffness tensor. For strain tensor, This represents the tensor contraction operation. For Biot coefficient, Pore ​​pressure, For unit tensors, This is the force vector; strain tensor It is calculated using the following formula: in, Let be the displacement vector of all nodes in the field. This is the transpose of the gradient operator; Substituting the formula for calculating the strain tensor into the constitutive equation, and then combining it with the effective stress principle into the quasi-static equilibrium equation, the intermediate variable total stress is eliminated. With strain The displacement is derived. The fluid-structure interaction deformation control equation for porous media with the main variable as the primary variable is: Based on the heterogeneous Young's modulus field The global stiffness matrix, which reflects the heterogeneity of the reservoir space, is obtained through unit integration assembly. The nodal pore pressure field calculated at the current time step It is mapped to an equivalent fluid nodal force vector through shape function interpolation. The boundary load conditions are then transformed into external load vectors. Establish the following system of linear algebraic equations: By solving the above system of linear algebraic equations, the total field nodal displacement vector at the current time step is obtained. , T For matrix transpose, u x and u y For nodes in x direction and y Displacement components in the direction; Obtaining the total field nodal displacement vector Then, the volumetric strain was calculated using Cauchy's strain geometry equations: in, For spatial location ( x , y ),time t The volumetric strain below, for u x right x Partial derivatives of coordinates, for u y right y Partial derivatives of coordinates; For each logarithmic time step of each sample, the above steps are performed to obtain the time-varying volumetric strain field corresponding to the input pressure field and the inhomogeneous parameter field, and the time-varying volumetric strain field is used as the standard truth label.

5. The iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 1, characterized in that, A neural network architecture containing multiple Fourier transform layers is constructed as a mechanical surrogate model. Pore pressure, Young's modulus, Poisson's ratio, and boundary loads are used as inputs, and volumetric strain is used as the output. The mechanical surrogate model is trained offline using the aforementioned sample dataset until convergence, including: The normalized two-dimensional pore pressure field and Young's modulus field are extracted as the basic channels; Poisson's ratio, overburden pressure and horizontal stress are extended into full-field feature maps as constant value channels using spatial broadcasting technology; deterministic position coding is introduced to generate normalized coordinate grid channels; all features are spliced ​​in the channel dimension to form a 7-channel composite input tensor. The input data is mapped to a high-dimensional feature space through a fully connected layer, and then enters a multi-layer stacked Fourier layer. Each layer adopts a dual-path architecture, and the two outputs are added together and then activated by a GELU layer. Finally, the high-dimensional features are regressed and mapped to the target volume strain field through a decoding projection layer. We use mean squared error as the loss function, use the Adam optimizer to update parameters, introduce a dynamic scheduling strategy, monitor the training set loss in real time, automatically decay the learning rate, and save the model parameters with the lowest loss as the mechanical surrogate model.

6. The iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 5, characterized in that, The space broadcasting technology includes: for each scalar Build a size of The full-field constant matrix As input; where Poisson's ratio, For the upper boundary load; For right boundary load; Deterministic position encoding is introduced to generate normalized coordinate grid channels; all features are concatenated along the channel dimension to form a 7-channel composite input tensor, including: Generate two first coordinate grid matrices with the same size as the physical grid. Second coordinate grid matrix Wherein, the first coordinate grid matrix The value increases linearly from 0 to 1 along the X-axis, while remaining constant along the Y-axis; this is the second coordinate grid matrix. The value increases linearly from 0 to 1 along the Y-axis, while the value remains unchanged along the X-axis. Based on the first coordinate grid matrix Second coordinate grid matrix A 7-channel composite input tensor is generated using the following formula: in, P For pore pressure, Concat is the feature splicing operation. E For Young's modulus field, This is the global constant matrix corresponding to Poisson's ratio. This is the global constant matrix corresponding to the overburden pressure or upper boundary load. This is the global constant matrix corresponding to the right boundary load.

7. The iterative accelerated computation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 5, characterized in that, The dual-path architecture includes a global path and a local path, wherein the global path performs spectral convolution operations; The local path is connected in parallel with a spatial convolutional layer, which directly performs linear combination of features in the spatial domain; The dynamic scheduling strategy includes: if the training loss does not decrease significantly within a set number of consecutive training rounds, a learning rate decay mechanism is triggered to halve the current learning rate.

8. The iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 1, characterized in that, A fluid-structure interaction (FSI) hybrid iterative acceleration computation framework based on constant stress splitting is constructed. This framework includes a constant stress splitting strategy to decouple the FSI system into independent fluid solution and solid inference modules. The fluid module uses the finite volume method to solve the pore pressure field, while the solid module uses a trained mechanical surrogate model to infer the volumetric strain field. A porosity evolution model is established as a two-way coupling link, feeding back the skeleton deformation to the fluid equations to drive the fluid solution and solid inference modules in alternating iterative closed-loop cycles until the physical field satisfies the convergence criterion, including: The fluid control equations are solved using the finite volume method to obtain the updated pore pressure field for the current iteration step. In constructing the mass accumulation term of the fluid equations, a dynamic porosity evolution model is coupled to reflect the feedback of solid deformation to the fluid domain. The updated pore pressure field, Young's modulus field characterizing the heterogeneity of the medium, and boundary load conditions are input into a pre-trained mechanical surrogate model to infer the predicted value of the volumetric strain field at the current moment, and the predicted value of the volumetric strain field is denormalized to restore its physical dimensions. The relative norm error between the pore pressure field obtained in the current iteration step and the previous iteration step is calculated and compared with a preset convergence tolerance. If the error is greater than the convergence tolerance, the predicted volumetric strain field obtained in the current iteration step is used to update the rock porosity through the dynamic porosity evolution model, and the solution of the fluid control equation and subsequent inference process are re-executed to start the next coupling iteration. If the error is less than or equal to the convergence tolerance, the fluid-structure interaction calculation in the current time step is determined to have converged, the physical field update of the time step is completed, and the process is advanced to the next time step. This process is repeated until all simulations are completed.

9. The iterative accelerated calculation method for porous media fluid-structure interaction based on Fourier neural operators according to claim 8, characterized in that, The dynamic porosity evolution model is expressed as follows: in, It is the volumetric strain predicted by the trained mechanical surrogate model at the current time step. It is the volumetric strain at the initial time step. It is the porosity at the current time step. It is the porosity at the initial time step. It is the Biot coefficient. C r It is the overall compression coefficient. It is the pore pressure at the current time step. P 0 It is the pore pressure at the initial time step; The method for performing inverse normalization on the predicted volumetric strain field values ​​to restore their physical dimensions is as follows: in, The output of the trained mechanical surrogate model, This is the result of dimensionless feature reconstruction of the output of the trained mechanical surrogate model. The standard deviation of the normalized Young's modulus. The mean of the normalized Young's modulus. E ref The baseline Young's modulus value, E ( x ) represents the Young's modulus field as input.

Images