Fluid-solid coupling seepage evolution simulation method for pore-fracture type oil and gas reservoir

Abstract

The present application belongs to oil and gas exploitation and numerical simulation technology, and relates to a pore-fracture type oil and gas reservoir fluid-solid coupling seepage evolution simulation method. Reservoir parameters and initial stress field are obtained, and a pore and fracture grid mapping relationship is established; the fracture normal displacement is calculated according to the stress field, the opening and permeability are updated, and a dynamic topology graph is generated; the fluid-solid coupling partial differential equation is discretized, the node features of the dynamic topology graph are used to represent the pore pressure and skeleton displacement, and the edge features are used to represent the fracture permeability and normal stiffness; a graph neural network is constructed, the dynamic topology graph is input, the balance residual of fluid pressure and skeleton deformation is used as a loss function, and the pore pressure field and stress field distribution are solved end to end and the results are output. The present application eliminates the time lag and error accumulation of the physical field alternating transmission, overcomes the oscillation caused by the fracture parameter mutation, solves the strong fluid-solid coupling alternating iteration divergence problem, and improves the numerical stability and calculation efficiency.

Description

Technical Field

[0001] This invention belongs to the field of oil and gas extraction and numerical simulation technology, and relates to a method for simulating the fluid-solid coupling seepage evolution of pore-fracture type oil and gas reservoirs. Background Technology

[0002] Oil and gas reservoirs contain both pore and fractured media, making fluid-structure interaction (FSI) seepage evolution simulation crucial for assessing production dynamics. Conventional techniques often employ a dual-media model, equating the fracture network to a continuous medium connecting pore blocks. In implementation, the reservoir geometry is meshed using the finite element method (FEM) or finite volume method, assigning porosity and permeability attributes to the pore matrix and fracture network, respectively. During computation, a crossflow coefficient is introduced between the pore and fracture systems. A system of partial differential equations is established based on mass and momentum conservation, and the pressure and stress fields are solved using an alternating iterative approach to characterize the coupling process between fluid flow and solid deformation. This method reduces modeling complexity through macroscopic averaging and is suitable for scenarios with relatively uniform fracture distribution where the mechanical behavior of individual fractures is not a concern.

[0003] When the fracture distribution exhibits strong heterogeneity and dominant connectivity, conventional techniques employ a discrete fracture network model. This approach explicitly characterizes fractures as low-dimensional geometric surfaces embedded within a porous matrix mesh. During the fluid-structure interaction (FSI) solution, matrix deformation control equations and fracture seepage control equations are established separately. The stress field is solved using the finite element method to obtain nodal displacements and fracture normal stresses, and empirical formulas are used to correct fracture aperture and permeability. The seepage field is solved by calculating fluid pressure distribution based on the corrected permeability, and then feeding pore pressure back as a volume load to the stress field calculation. The two physical fields transfer variables to each other through an alternating iterative mechanism until the change in physical quantities in adjacent iteration steps is less than a set threshold, completing one time step of the coupled calculation. This iterative process is repeated to obtain the seepage evolution results.

[0004] The aforementioned discrete fracture network model, combined with the conventional alternating iteration approach, is highly susceptible to numerical oscillations and computational divergence under strong fluid-structure interaction (FSI) conditions. Under stress, fractures undergo normal closure or shear slip, leading to abrupt changes in aperture and permeability across orders of magnitude, resulting in dynamic evolution of the topology. When dealing with such strong nonlinearity and topological abrupt changes, the alternating iteration mechanism suffers from a time lag effect in the transfer of physical quantities between the seepage field and the stress field. The computational error in a particular field is amplified during iterations, violating conservation laws and preventing the entire numerical calculation process from converging to a stable state. Therefore, the core technical problem of existing techniques lies in the ease with which alternating iteration solutions diverge under strong FSI conditions due to the dynamic topological evolution of fractures. Summary of the Invention

[0005] The purpose of this invention is to provide a fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs, which can solve the problems in the background art mentioned above.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows: A fluid-structure interaction simulation method for seepage evolution in pore-fracture type oil and gas reservoirs includes: obtaining reservoir geometric parameters and initial stress field, and establishing the mapping relationship between pore continuous medium grid and fracture discrete medium grid; The normal displacement of the fracture is calculated based on the stress field distribution, the fracture aperture and permeability are updated, and a dynamic topology map characterizing the connectivity between pores and fractures is generated. The mass and momentum conservation equations in the fluid-structure interaction partial differential equations are discretized, and the pore pressure and skeleton displacement are characterized by the node features of the dynamic topology graph, while the fracture permeability and normal stiffness are characterized by the edge features. A graph neural network is constructed, and the dynamic topology graph is input into the graph neural network. The balance residual between fluid pressure and skeleton deformation is used as the loss function to solve the pore pressure field and stress field distribution end-to-end, and the seepage evolution results are output.

[0007] Preferably, establishing the mapping relationship between the porous continuous medium grid and the fractured discrete medium grid includes: calculating the shortest spatial distance from the center point of the porous continuous medium grid to the surface of the fractured discrete medium grid, and constructing a spatial proximity matrix; Extract the set of boundary nodes of the porous continuous medium mesh, and project the set of boundary nodes onto the normal direction of the fractured discrete medium mesh to obtain the set of projection intersection points; Based on the spatial proximity matrix and the set of projection intersections, equivalent weighting coefficients for cross-scale mass and momentum transfer are constructed, and a bidirectional mapping relationship between the porous continuous medium grid and the fractured discrete medium grid is established according to the equivalent weighting coefficients.

[0008] Preferably, the calculation of the fracture normal displacement based on the stress field distribution, updating the fracture aperture and permeability, and generating a dynamic topology map characterizing the connectivity between pores and fractures includes: extracting the effective normal stress and tangential shear stress of the fracture discrete medium mesh based on the stress field distribution, and calculating the normal closure amount and shear slip amount of the fracture discrete medium mesh. Based on the normal closure amount and the shear slip amount, the Mohr-Coulomb criterion is introduced to determine the contact state of the fractured discrete medium mesh, and the contact state includes fully open, slip, and locked. The crack aperture and permeability are dynamically corrected based on the contact state, and the corrected crack aperture and permeability are assigned to the edge attributes of the topology graph to generate the dynamic topology graph.

[0009] Preferably, the discretization of the mass conservation and momentum conservation equations in the fluid-structure interaction partial differential equations, and the characterization of pore pressure and skeleton displacement by the node features of the dynamic topology graph, and the characterization of fracture permeability and normal stiffness by the edge features, includes: using the control volume finite element method to spatially discretize the mass conservation equations and the momentum conservation equations to obtain the pore node control volume and the fracture edge control volume. Asynchronous time steps are set for the pore node control body and the crack edge control body respectively, and time difference is performed on the discretized equations based on the asynchronous time steps; The differential pore pressure and skeleton displacement variables are mapped to the node feature vectors of the dynamic topology graph, and the differential fracture permeability and normal stiffness variables are mapped to the edge feature vectors of the dynamic topology graph.

[0010] Preferably, the construction of the graph neural network involves inputting the dynamic topology graph into the graph neural network, using the balance residual between fluid pressure and skeleton deformation as the loss function, solving the pore pressure field and stress field distribution end-to-end, and outputting the seepage evolution results. This includes: in the graph neural network, calculating the message passing weights between adjacent nodes based on the edge feature vectors of the dynamic topology graph, performing weighted aggregation and nonlinear transformation on the node feature vectors, and obtaining the hidden layer features of the nodes. The hidden layer features of the nodes are residually connected with the original node feature vectors to output the predicted pore pressure and skeleton displacement values. The equilibrium residuals of the predicted pore pressure and the predicted skeleton displacement with respect to the mass conservation equation and the momentum conservation equation are calculated. A loss function is constructed using the norm of the equilibrium residuals. The parameters of the graph neural network are updated by backpropagation until the loss function converges, and the seepage evolution results are output.

[0011] Preferably, projecting the set of boundary nodes onto the normal direction of the fractured discrete medium mesh to obtain a set of projection intersection points includes: calculating the projection line of each boundary node in the set of boundary nodes along the normal direction of the fractured discrete medium mesh; Determine the intersection state between the projected line and the fractured discrete medium mesh surface. If there are multiple intersection points, select the intersection point closest to the boundary node as the valid intersection point. If the projected line does not intersect the fractured discrete medium mesh surface, then the nearest face point of the boundary node on the fractured discrete medium mesh surface is searched as a substitute intersection point. Aggregate all the valid intersection points and the alternative intersection points to construct the projected intersection point set.

[0012] Preferably, the step of dynamically correcting the fissure aperture and permeability based on the contact state includes: when the contact state is fully open, setting the permeability to a first correlation value that is proportional to the cube of the fissure aperture according to the cubic law; When the contact state is slippage, the shear dilatation increment is calculated based on the shear dilatation angle, and the permeability is set to a second correlation value that is proportional to the cube of the sum of the fracture aperture and the shear dilatation increment. When the contact state is locked, a residual opening threshold is introduced, and the permeability is set to a third correlation value that is proportional to the cube of the residual opening threshold. The fracture opening and permeability are dynamically corrected using the first correlation value, the second correlation value, and the third correlation value.

[0013] Preferably, asynchronous time steps are set for the pore node control body and the fracture edge control body respectively, and time difference is performed on the discretized equations according to the asynchronous time steps, including: adaptively determining the first pore time step of the pore node control body and the first fracture time step of the fracture edge control body according to the rate of change of fluid pressure and skeleton displacement. When the time step of the first pore is not synchronized with the time step of the first crack, an interpolation method is used on the control volume with a small time scale to obtain the physical quantity of the control volume with a large time scale at the boundary moment. Based on the physical quantities obtained by the interpolation, flux compensation calculations are performed on the control volume with a large time scale so that the pore node control volume and the crack edge control volume satisfy local conservation within their respective time steps.

[0014] Preferably, before acquiring the reservoir geometric parameters and initial stress field, the method further includes: acquiring well logging data and seismic anisotropy data of the target reservoir, and extracting the shear wave time difference from the well logging data and the azimuth anisotropy amplitude from the seismic anisotropy data. A shear wave splitting model is constructed, and the shear wave time difference and the azimuth anisotropic amplitude are input into the shear wave splitting model to invert and obtain the initial maximum horizontal principal stress direction and the initial minimum horizontal principal stress direction of the target reservoir. Boundary loads are applied to the geomechanical model based on the initial maximum horizontal principal stress direction and the initial minimum horizontal principal stress direction, and the initial stress field is obtained by calibration.

[0015] Preferably, after outputting the seepage evolution results, the method further includes: extracting the pore pressure field distribution and fracture velocity field distribution at the target time from the seepage evolution results; Microscale fluid tracer testing was performed on the physical model of the target reservoir to obtain tracer breakthrough curves; The fracture velocity field distribution is converted into a simulated breakthrough curve, and the dynamic error between the simulated breakthrough curve and the tracer breakthrough curve is calculated. If the dynamic error exceeds a preset threshold, the dynamic error is added as a penalty to the loss function of the graph neural network, and the parameters of the graph neural network are updated again until the dynamic error is less than the preset threshold.

[0016] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. The mapping relationship between the porous continuous medium mesh and the fractured discrete medium mesh is transformed into a dynamic topological graph structure. A graph neural network is then used to solve the pore pressure and stress field distributions end-to-end, replacing the traditional alternating iterative solution process for fluid-structure interaction (FSI). By discretizing the mass and momentum conservation equations in the FSI partial differential equations, the node features of the dynamic topological graph represent pore pressure and skeleton displacement, while the edge features represent fracture permeability and normal stiffness, thus constructing an isomorphic mapping between physical field variables and graph structure data. In the graph neural network computation, the dynamic topological graph is input into the network, and a loss function constructed from the equilibrium residuals of fluid pressure and skeleton deformation is used as a constraint, forcing the network to satisfy the physical conservation laws of FSI during parameter updates. This end-to-end solution mechanism eliminates the time lag and error accumulation caused by the alternating transmission of the seepage field and stress field, overcomes the physical quantity oscillation caused by the sudden change of crack aperture and permeability across orders of magnitude, solves the technical problem that alternating iterative solution is prone to divergence under strong fluid-structure interaction, and improves the numerical stability and computational efficiency of fluid-structure interaction seepage evolution simulation of discontinuous media.

[0017] 2. When updating the fracture aperture and permeability to generate a dynamic topology map, the effective normal stress and tangential shear stress of the fracture discrete medium mesh are extracted based on the stress field distribution. The normal closure amount and shear slip amount are calculated, and the Mohr-Coulomb criterion is introduced to determine the contact state of the fracture being fully open, slipping, and locked. Based on different contact states, the aperture and permeability are dynamically corrected and assigned to the edge attributes of the topology map, realistically reproducing the mechanical closure and shear dilatation behavior of the fracture under different stress paths. During the equation discretization process, asynchronous time steps are set for the pore node control volume and the fracture edge control volume, respectively. When the time scales are mismatched, interpolation is used to obtain the physical quantities at the boundary moment for flux compensation calculation, adapting to the objective difference in seepage rate between the pore matrix and the fracture network, reducing unnecessary small time step calculation overhead while maintaining local mass conservation. Inside the graph neural network, the message passing weights between adjacent nodes are calculated based on edge features, and weighted aggregation and nonlinear transformation are performed. The hidden features of the nodes are then residually connected with the original node feature vectors, which avoids gradient vanishing in deep network training and ensures the physical fidelity of solving fluid-structure interaction partial differential equations. Attached Figure Description

[0018] Figure 1 The overall main flow chart for simulating the fluid-structure interaction seepage evolution of pore-fracture type oil and gas reservoirs; Figure 2 A flowchart for constructing a two-way mapping relationship between porous continuous medium meshes and fractured discrete medium meshes; Figure 3 A flowchart for discretization of fluid-structure interaction equations, graph eigenvector mapping, and asynchronous time step processing; Figure 4 This is a flowchart of the end-to-end solution and parameter iterative update process for a graph neural network. Detailed Implementation

[0019] refer to Figures 1 to 4 In one embodiment, reservoir geometric parameters and an initial stress field are obtained, and a mapping relationship is established between the pore continuous medium grid and the fracture discrete medium grid. The reservoir geometric parameters include the three-dimensional spatial extent of the reservoir, the porosity distribution of the pore matrix, and the geometric morphology parameters and spatial distribution coordinates of the fractures. The initial stress field includes the spatial distribution of the vertical principal stress, the maximum horizontal principal stress, and the minimum horizontal principal stress. The pore continuous medium grid uses tetrahedral or hexahedral elements to partition the reservoir matrix, with each grid element corresponding to a pore control volume. The fracture discrete medium grid uses triangular or quadrilateral facet elements to partition a single fracture, with each facet element corresponding to a fracture control volume. When establishing the mapping relationship, the spatial association between each pore continuous medium grid element and all fracture discrete medium grid facet elements is calculated to determine the transmission path and intensity of fluid and stress between the two media.

[0020] The normal displacement of fractures is calculated based on the stress field distribution, and the fracture aperture and permeability are updated to generate a dynamic topology map characterizing the connectivity between pores and fractures. The stress field distribution is formed by superimposing the initial stress field and the subsequently calculated skeleton deformation field. The normal and tangential stress components on each surface element of the fracture discrete medium mesh are extracted, and the normal closing displacement and tangential slip displacement of the fracture surface are calculated according to the rock mechanics constitutive relation. The fracture aperture is equal to the initial aperture minus the normal closing displacement. When shear slip occurs, an additional aperture increment is added according to the dilatation effect. The permeability is calculated using the cubic law based on the corrected fracture aperture. The nodes of the dynamic topology map correspond to the nodes of the pore continuous medium mesh and the nodes of the fracture discrete medium mesh, and the edges of the dynamic topology map correspond to the element edges of the pore continuous medium mesh and the surface element edges of the fracture discrete medium mesh. The existence of an edge indicates that there is a fluid connection channel between two nodes, and the attribute values ​​of the edge represent the permeability and mechanical stiffness of the connection channel.

[0021] The mass and momentum conservation equations in the fluid-structure interaction partial differential equations (FSIs) are discretized. The pore pressure and skeleton displacement are characterized by node features in the dynamic topology graph, while the fracture permeability and normal stiffness are characterized by edge features. The FSIs describe the flow patterns of fluids in pores and fractures, as well as the deformation of the solid skeleton under fluid pressure and geostress. The mass conservation equation ensures that the mass difference between the fluid flowing into and out of the control volume equals the rate of change of fluid mass within the control volume. The momentum conservation equation ensures that the resultant force of all forces acting on the solid skeleton equals the inertial force of the skeleton. A numerical discretization method is used to transform the continuous FSIs into a discrete system of algebraic equations, converting continuous physical field variables into physical quantities at discrete grid nodes. The discretized pore pressure and skeleton displacement variables are combined into node eigenvectors, and the discretized fracture permeability and normal stiffness variables are combined into edge eigenvectors, achieving a correspondence between the physical field variables and the graph structure data.

[0022] A graph neural network (Graph Neural Network) is constructed by inputting a dynamic topological graph into the network. The residuals of the equilibrium between fluid pressure and skeleton deformation are used as the loss function to solve for the pore pressure and stress field distributions end-to-end, outputting the seepage evolution results. The Graph Neural Network consists of multiple stacked graph convolutional layers. Each convolutional layer receives node and edge features from the previous layer and updates the node features through a message passing mechanism. The message passing mechanism calculates the information transfer weights between adjacent nodes based on the edge features, weights and aggregates the features of adjacent nodes, and then obtains the updated node features through a nonlinear activation function. The loss function is a weighted sum of the residuals of the mass conservation equation and the momentum conservation equation, where the residuals represent the deviation between the predicted physical quantity values ​​and the theoretical values ​​required by the conservation equations. The gradient of the loss function with respect to the Graph Neural Network parameters is calculated using the backpropagation algorithm, and the network parameters are updated using gradient descent until the loss function converges to a preset threshold. The converged network output represents the pore pressure and skeleton displacement field distributions at the current moment. This is used as the initial condition for the next moment, and the above process is repeated to obtain the seepage evolution results over the entire time range.

[0023] In this embodiment, reservoir geometric parameters are extracted from 3D seismic and well logging data using geological modeling software. The initial stress field is obtained by combining geostress test data with geomechanical simulation. The mesh density of the pore continuum grid is determined based on the degree of reservoir heterogeneity; a denser mesh is used in areas with drastic changes in reservoir properties, while a sparser mesh is used in areas with relatively uniform properties. The mesh density of the fracture discrete medium grid is determined based on the geometric complexity of a single fracture; fewer surface elements are used for long fractures with regular shapes, while more surface elements are used for short fractures with complex shapes. The number of nodes in the dynamic topology graph is equal to the sum of the number of nodes in the pore continuum grid and the number of nodes in the fracture discrete medium grid, and the number of edges in the dynamic topology graph is equal to the sum of the number of cell edges in the pore continuum grid and the number of surface cell edges in the fracture discrete medium grid. The number of layers in the graph neural network is set to 3 to 5, and the number of hidden units in each layer is set to 64 to 256. The convergence threshold of the loss function is set to 10. -6 Up to 10 -4 .

[0024] This embodiment transforms the coupling relationship between a porous continuous medium and a fractured discrete medium into a dynamic topological graph structure. It utilizes the message passing mechanism of a graph neural network to simulate fluid flow between pores and fractures, as well as stress transmission within the solid framework. The residuals of the fluid-structure interaction conservation equations are used as a loss function to constrain network training, ensuring that the network output automatically satisfies physical conservation laws. The end-to-end solution method avoids the lag and error accumulation in the transmission of physical quantities between the seepage field and stress field in traditional alternating iterative methods, thus improving the stability of numerical calculations under strong fluid-structure interaction conditions.

[0025] In a preferred embodiment, when establishing the mapping relationship between the porous continuous medium mesh and the fractured discrete medium mesh, the shortest spatial distance from the center point of the porous continuous medium mesh to the surface of the fractured discrete medium mesh is calculated, and a spatial proximity matrix is ​​constructed. The center point of the porous continuous medium mesh is the geometric center coordinate of each tetrahedral or hexahedral element. The surface of the fractured discrete medium mesh is the plane containing each triangular or quadrilateral surface element. The shortest spatial distance is the perpendicular distance from the mesh center point to the plane. When the foot of the perpendicular falls inside the surface element, the perpendicular distance is the shortest distance; when the foot of the perpendicular falls outside the surface element, the shortest distance is the shortest distance from the mesh center point to the boundary of the surface element. The rows of the spatial proximity matrix correspond to the porous continuous medium mesh elements, and the columns correspond to the fractured discrete medium mesh surface elements. The values ​​of the matrix elements are the reciprocals of the corresponding shortest distances. The closer the distance, the larger the element value, indicating a stronger spatial association between the two elements.

[0026] The boundary node set of the porous continuous medium mesh is extracted, and then projected onto the normal direction of the fractured discrete medium mesh to obtain the projection intersection set. The boundary nodes of the porous continuous medium mesh are the pore mesh nodes located at the reservoir geometric boundary or intersecting with fractures. The normal direction of the fractured discrete medium mesh is the unit normal vector of each surface element. The projection process involves drawing a straight line from the boundary node along the fracture normal direction and calculating the intersection points of this line with the fracture surface elements. If multiple intersection points exist, the intersection point closest to the boundary node is selected as the valid intersection point; if no intersection point exists, the nearest point of the boundary node on the fracture surface element is searched as the alternative intersection point. All valid intersection points and alternative intersection points are collected to form the projection intersection set.

[0027] Based on the spatial proximity matrix and the set of projection intersection points, equivalent weighting coefficients for cross-scale mass and momentum transfer are constructed. A bidirectional mapping relationship is established between the porous continuous medium mesh and the fractured discrete medium mesh based on these equivalent weighting coefficients. The equivalent weighting coefficients are obtained by weighted summation of the spatial proximity weights and the projection intersection point weights. The spatial proximity weights are obtained by normalizing the elements of the spatial proximity matrix, reflecting the natural transfer capability between the two elements due to spatial proximity. The projection intersection point weights are jointly determined by the position of the projection intersection point on the fracture surface element and the position of the boundary node on the porous mesh element, reflecting the forced transfer capability between the two elements due to geometric intersection. The bidirectional mapping relationship includes mappings from porous media to fractured media and mappings from fractured media to porous media. The mapping from porous media to fractured media transfers the physical quantity values ​​at the porous mesh nodes to the fractured mesh nodes through the equivalent weighting coefficients; conversely, the mapping from fractured media to porous media transfers the physical quantity values ​​at the fractured mesh nodes to the porous mesh nodes through the equivalent weighting coefficients.

[0028] In this embodiment, the calculation of the shortest spatial distance is achieved using vector dot product and cross product operations. For a triangular facet element, the coordinate vectors of the three vertices are first calculated, and then the normal vector of the facet element and the plane equation are calculated. For a given spatial point, the perpendicular distance from the point to the plane is calculated, and it is determined whether the foot of the perpendicular falls inside the triangle. If the foot of the perpendicular falls inside the triangle, the perpendicular distance is the shortest distance; if the foot of the perpendicular falls outside the triangle, the shortest distances from the point to the three sides of the triangle are calculated separately, and the minimum value is taken as the final shortest distance. For a quadrilateral facet element, it is decomposed into two triangular facet elements, and the shortest distances from the spatial point to the two triangles are calculated separately, and the minimum value is taken as the shortest distance to the quadrilateral facet element.

[0029] The spatial proximity matrix is ​​normalized using row-wise normalization. For each row, all elements are divided by the sum of all elements in that row, ensuring the sum of all elements in that row equals 1. This ensures that the total weight transferred from a single pore grid cell to all fractured grid cells is 1, satisfying the conservation of mass and momentum. The weight of the projected intersection is calculated using an area-weighted average method. For each projected intersection, the area of ​​the fractured grid cell containing that intersection is calculated as the ratio of the area of ​​that intersection to the total area of ​​all fractured grid cells containing the intersection. This ratio is then used as the weight of that intersection.

[0030] The formula for calculating the equivalent weighting coefficient is: ; in, For the first The first pore grid cell and the first The equivalent weighting coefficients between the surface elements of the fractured mesh This is a balance coefficient between the spatial proximity weight and the projection intersection weight, with a value ranging from 0 to 1. For the first From the center point of the first pore grid cell to the... The shortest spatial distance of each fractured mesh surface element. For the first From the center point of the first pore grid cell to the... The shortest spatial distance of each fractured mesh surface element. The total number of surface elements in the fractured mesh. For the first The area of ​​each fractured mesh element. For the first The area of ​​each fractured mesh element. For including the first A collection of cracked mesh surface elements at the intersection of the projected boundary nodes of a porous mesh element.

[0031] This embodiment constructs a bidirectional mapping relationship by combining a spatial proximity matrix and a set of projected intersection points. This considers both the spatial positional relationship between porous and fractured media, as well as their geometric intersection relationships. The introduction of equivalent weighting coefficients quantifies the transfer intensity between different units, making cross-scale mass and momentum transfer more accurate. The establishment of the bidirectional mapping relationship achieves seamless coupling between the porous continuous medium and the fractured discrete medium, laying the foundation for the subsequent generation of dynamic topology graphs and the solution of graph neural networks.

[0032] In a preferred embodiment, when calculating the fracture normal displacement based on the stress field distribution, updating the fracture aperture and permeability, and generating a dynamic topology map characterizing the connectivity between pores and fractures, the effective normal stress and tangential shear stress of the fracture discrete medium mesh are extracted based on the stress field distribution, and the normal closure and shear slip of the fracture discrete medium mesh are calculated. The effective normal stress is the total normal stress acting on the fracture surface minus the pore fluid pressure, reflecting the actual contact stress between the solid skeletons. The tangential shear stress is the stress component parallel to the fracture surface acting on it. The normal closure is the relative displacement of the fracture surface under the action of the effective normal stress, and it has a nonlinear relationship with the effective normal stress. The shear slip is the relative tangential displacement of the fracture surface under the action of the tangential shear stress; shear slip occurs when the tangential shear stress exceeds the shear strength of the fracture.

[0033] Based on normal closure and shear slip, the Mohr-Coulomb criterion is introduced to determine the contact state of a fractured discrete medium mesh. The contact states include fully open, slip, and locked. The Mohr-Coulomb criterion describes the shear failure condition of rock materials, and its expression is: ; in, The shear strength of the crack. The cohesive force within the fracture. For the normal effective stress, The internal friction angle of the crack is denoted as . When the effective normal stress is less than or equal to zero, the crack is fully open, with no contact between the crack surfaces and zero shear strength. When the effective normal stress is greater than zero and the tangential shear stress is less than or equal to the shear strength, the crack is closed, with close contact between the crack surfaces and no relative slippage. When the effective normal stress is greater than zero and the tangential shear stress is greater than the shear strength, the crack is in a slipping state, with relative shear slippage occurring between the crack surfaces.

[0034] The fracture aperture and permeability are dynamically adjusted based on the contact state, and the adjusted aperture and permeability are assigned to the edge attributes of the topology graph to generate a dynamic topology graph. When the fracture is fully open, the fracture aperture equals the initial aperture, and the permeability is calculated according to the cubic law. When the fracture is in a slipping state, the fracture aperture equals the initial aperture minus the normal closure amount plus the shear dilatation aperture increment, where the shear dilatation aperture increment is proportional to the shear slip and the shear dilatation angle. When the fracture is in a closed state, the fracture aperture equals the residual aperture, which is the minimum aperture that the fracture surface cannot completely close even under high pressure, and the permeability is calculated based on the residual aperture. The adjusted fracture aperture and permeability are used as attribute values ​​for the corresponding fracture edges in the dynamic topology graph, and the existence state of the edges is updated simultaneously. When the fracture aperture is less than or equal to zero, the corresponding edge is deleted, indicating that the connection channel is closed; when the fracture aperture is greater than zero, the corresponding edge is retained, indicating that the connection channel is open.

[0035] In this embodiment, the calculation of the normal closure amount adopts a hyperbolic constitutive model, the expression of which is: ; in, For normal closed-loop quantity, The initial normal stiffness of the crack is... The maximum normal closure of the crack is equal to the initial opening. This model describes the nonlinear characteristic of the crack normal stiffness increasing with the increase of the effective normal stress. The shear slip is calculated using an elastoplastic constitutive model. When the tangential shear stress is less than the shear strength, the shear deformation is elastic, and the shear slip is proportional to the tangential shear stress. When the tangential shear stress is greater than the shear strength, the shear deformation is plastic, and the shear slip continuously increases with the increase of the tangential shear stress.

[0036] The formula for calculating the shear dilatation increment is: ; in, This is the increment of shear dilatation opening. This is the shear slip. This is the shear dilatation angle of the fracture. The shear dilatation angle reflects the volume expansion effect of the fracture surface due to its roughness during shearing.

[0037] This embodiment introduces the Mohr-Coulomb criterion to accurately determine the three contact states of the fracture. Different aperture and permeability correction formulas are used for different contact states, realistically reproducing the mechanical behavior of the fracture under complex stress conditions. The edge attributes and existence states of the dynamic topology graph are updated in real time with the changes in the fracture contact state, accurately characterizing the dynamic evolution process of the connectivity between pores and fractures. This dynamic topology provides an accurate geometric and physical basis for graph neural networks to solve fluid-structure interaction problems.

[0038] In a preferred embodiment, the mass and momentum conservation equations in the fluid-structure interaction partial differential equations are discretized. The pore pressure and skeleton displacement are characterized by node features of the dynamic topology graph, and the fracture permeability and normal stiffness are characterized by edge features. The control volume finite element method is then used to spatially discretize the mass and momentum conservation equations, obtaining the pore node control volume and the fracture edge control volume. The control volume finite element method combines the geometric adaptability of the finite element method with the conservation properties of the finite volume method. For porous continuous media, the region surrounding each mesh node is divided into a node control volume, and the boundary of the control volume is formed by connecting the midpoint of the edge of the adjacent element to the center of the element. For fracture discrete media, each fracture edge element is divided into an edge control volume, and the boundary of the control volume is formed by connecting the two endpoints of the fracture edge element to the midpoint of the adjacent element.

[0039] Asynchronous time steps are set for the pore node control volume and the fracture edge control volume, respectively, and time difference is performed on the discretized equations based on the asynchronous time steps. The fluid seepage velocity in the porous matrix is ​​relatively slow, and the changes in physical quantities are relatively gradual, allowing for a larger time step. The fluid seepage velocity in the fracture network is relatively fast, and the changes in physical quantities are more drastic, requiring a smaller time step. The asynchronous time step allows the pore control volume and the fracture control volume to be updated at different time points, accommodating the difference in seepage rates between the two media. The time difference employs an implicit difference scheme, improving the stability of the numerical calculation.

[0040] The differentially derived pore pressure and skeleton displacement variables are mapped to node feature vectors of the dynamic topology graph, and the differentially derived fracture permeability and normal stiffness variables are mapped to edge feature vectors of the dynamic topology graph. The node feature vector is a two-dimensional vector, with the first element being the pore pressure value and the second element being the skeleton displacement value. For pore nodes, the skeleton displacement value is the modulus of the displacement components of the node in the three coordinate axes; for fracture nodes, the skeleton displacement value is the displacement component of the fracture surface in the normal direction. The edge feature vector is a two-dimensional vector, with the first element being the fracture permeability value and the second element being the fracture normal stiffness value. For pore edges, the permeability value is the permeability of the pore matrix, and the normal stiffness value is the bulk modulus of the pore matrix; for fracture edges, the permeability value is the corrected fracture permeability, and the normal stiffness value is the corrected fracture normal stiffness.

[0041] In this embodiment, the spatial discretization of the mass conservation equation is performed using the control volume integral method. The law of mass conservation is applied to each control volume to obtain the discretized algebraic equation: ; in, For the boundary of the control volume, For fluid density, For fluid at the boundary surface Normal flow velocity on the surface For boundary surface The unit normal vector, For boundary surface area, Porosity To control the volume of the body.

[0042] The spatial discretization of the momentum conservation equation is performed using the finite element method. By applying the principle of virtual work to each element, the discretized algebraic equation is obtained: ; in, The element stiffness matrix, Let be the nodal displacement vector. For external load vector, This is the equivalent nodal force vector generated by pore pressure.

[0043] The time difference is performed using a backward Euler scheme. For the pore pressure variable, the discretized equation is: ; in, for The pore pressure vector at time t. for The pore pressure vector at time t. For time step, The discrete Laplace operator matrix, This is the source term vector.

[0044] This embodiment employs the controlled volume finite element method for spatial discretization, ensuring the conservation of numerical computation and geometric adaptability. The introduction of asynchronous time steps accommodates the difference in seepage rates between the porous matrix and the fracture network, improving computational efficiency while maintaining accuracy. Mapping the discretized physical quantities to the node and edge features of a graph achieves an organic combination of traditional numerical discretization methods and graph neural network methods, providing a unified data representation for end-to-end solutions to fluid-structure interaction problems.

[0045] In a preferred embodiment, a graph neural network is constructed. A dynamic topology graph is input into the graph neural network, and the balance residual between fluid pressure and skeleton deformation is used as the loss function. The pore pressure and stress field distributions are solved end-to-end. When outputting the seepage evolution results, the message passing weights between adjacent nodes are calculated based on the edge feature vectors of the dynamic topology graph in the graph neural network. Weighted aggregation and nonlinear transformation are performed on the node feature vectors to obtain the hidden layer features of the nodes. The message passing weights are calculated by passing the edge feature vectors through a two-layer fully connected neural network, and the output is a scalar weight value. The weighted aggregation process involves multiplying the feature vectors of all adjacent nodes of each node by their corresponding message passing weights and then summing the results to obtain the aggregated feature vector of that node. The nonlinear transformation process involves inputting the aggregated feature vector into a fully connected neural network with an activation function to obtain the hidden layer feature vector of that node.

[0046] The hidden layer features of a node are residually concatenated with the original node feature vector to output predicted pore pressure and skeleton displacement values. The residual concatenation involves element-wise addition of the input and output feature vectors from the previous layer before inputting them into the next layer. Residual concatenation can alleviate the vanishing gradient problem in deep neural network training, improving the network's training stability and convergence speed. After processing through multiple graph convolutional layers, the final layer outputs the predicted feature vector of the node, with the first element being the predicted pore pressure value and the second element being the predicted skeleton displacement value.

[0047] The equilibrium residuals of predicted pore pressure and skeleton displacement values ​​relative to the mass and momentum conservation equations are calculated. A loss function is constructed using the norm of the equilibrium residuals. Backpropagation is used to update the parameters of the graph neural network until the loss function converges, at which point the seepage evolution results are output. The mass conservation residual is the residual obtained by substituting the predicted pore pressure value into the discretized mass conservation equation. The momentum conservation residual is the residual obtained by substituting the predicted skeleton displacement value into the discretized momentum conservation equation. The loss function is the weighted sum of the L2 norms of the mass and momentum conservation residuals. The gradient of the loss function with respect to all trainable parameters of the network is calculated using automatic differentiation, and the network parameters are updated using the Adam optimizer. When the value of the loss function is less than a preset convergence threshold, training stops, and the current pore pressure and skeleton displacement field distributions are output.

[0048] In this embodiment, the expression for the message passing process of the graph convolutional layer is: ; in, For the first Layer The feature vector of each node For the first Layer The feature vector of each node For the first The set of neighboring nodes of a node. For the first Layer The node and the first Message passing weights between nodes For the first The trainable weight matrix of the layer, For the first The trainable bias vectors of the layer. The activation function is non-linear, using the ReLU function.

[0049] The formula for calculating message passing weight is: ; in, To connect the first The node and the first The feature vectors of the edges of each node. For the first The two-layer fully connected neural network has 32 hidden units in the first layer and 1 output unit in the second layer. The activation function is the Sigmoid function, which restricts the output value to between 0 and 1.

[0050] The formula for calculating the loss function is: ; in, For loss function, and represents the weighting coefficients for the mass conservation residual and the momentum conservation residual. Let the mass-conserved residual vector be... The momentum-conserving residual vector, This is the pore pressure prediction vector. This is the skeleton displacement prediction vector.

[0051] The graph neural network constructed in this embodiment simulates the propagation and interaction of physical fields in space through a message-passing mechanism. The introduction of residual connections improves the network's depth and expressive power. Using the residuals of the physical conservation equations as the loss function allows the network to automatically learn solutions that satisfy physical laws during training, without requiring a large amount of labeled data. The end-to-end solution approach unifies spatial discretization, temporal discretization, and equation solving within a single neural network framework, improving computational efficiency and numerical stability.

[0052] In a preferred embodiment, the set of boundary nodes is projected onto the normal direction of the fractured discrete medium mesh. When obtaining the set of projection intersection points, the projection lines of each boundary node in the set of boundary nodes along the normal direction of the fractured discrete medium mesh are calculated. The coordinates of the boundary nodes are three-dimensional spatial coordinates. The unit normal vector of the fracture surface element is Then the parametric equation of the projected line is: ; in, For scalar parameters, when When, the corresponding boundary node position; when When, it extends along the positive direction of the normal; when At that time, it extends in the negative direction along the normal.

[0053] Determine the intersection state between the projected line and the fractured discrete medium mesh surface. If multiple intersection points exist, select the intersection point closest to the boundary node as the valid intersection point. For triangular surface elements, substitute the parametric equations of the projected line into the equations of the plane containing the triangle to obtain the solution. The value of the intersection point is used to obtain the coordinates of the intersection point. Then, it is determined whether the intersection point lies inside the triangle. If it does, it is a valid intersection point; if it falls outside, it is not. For quadrilateral surface elements, the parametric equations of the projected line are substituted into the equations of the plane containing the quadrilateral to obtain the coordinates of the intersection point. Then, it is determined whether the intersection point lies inside the quadrilateral. If the projected line intersects multiple fracture surface elements, the distance from each valid intersection point to the boundary node is calculated, and the intersection point with the smallest distance is selected as the final valid intersection point.

[0054] If the projected line does not intersect the fractured discrete medium mesh surface, the nearest facet point of the boundary node on the fractured discrete medium mesh surface is used as the alternative intersection point. The nearest facet point is the point on the fractured surface element that is closest to the boundary node, and its calculation method is the same as the method for calculating the shortest distance in space. For triangular surface elements, the nearest facet point may be a vertex, a point on an edge, or an interior point of the triangle. For quadrilateral surface elements, the nearest facet point may be a vertex, a point on an edge, or an interior point of the quadrilateral. All valid intersection points and alternative intersection points are collected to form the projected intersection point set.

[0055] In this embodiment, the centroid coordinate method is used to determine whether a point lies inside a triangle. For triangle ABC and point P, the vector is calculated. , and Then calculate the three centroid coordinate parameters: like , and If point P is inside triangle ABC, then point P is outside triangle ABC; otherwise, point P is outside triangle ABC.

[0056] To determine if a point lies inside a quadrilateral, use the ray method. Starting from point P, draw a ray in any direction and count the number of intersections between the ray and each of the four sides of the quadrilateral. If the number of intersections is odd, point P lies inside the quadrilateral; if the number of intersections is even, point P lies outside the quadrilateral. It is important to handle the special case where the ray intersects with a vertex of the quadrilateral to avoid double counting.

[0057] This embodiment details the specific implementation method of projecting boundary nodes onto the fracture surface, including the parametric equations of the projected lines, the solution and validity judgment of intersection points, and the search method for the nearest face point when there are no intersection points. These methods ensure the accuracy and completeness of the set of projected intersection points, providing a reliable geometric foundation for the subsequent construction of cross-scale equivalent weighting coefficients.

[0058] In a preferred embodiment, when dynamically adjusting the fracture aperture and permeability based on the contact state, when the contact state is fully open, the permeability is set to a first correlation value proportional to the cube of the fracture aperture, according to the cubic law. The cubic law describes the seepage law of fluid in a parallel plate fracture, and its expression is: ; in, For fracture permeability, Let f be the fracture aperture. This law assumes that the fracture surfaces are smooth parallel plates and that the fluid flows in a laminar manner. For a fully open fracture, there is no contact between the fracture surfaces, and the fluid can flow freely; therefore, the cubic law is used to calculate permeability.

[0059] When the contact state is slip, the shear dilatation aperture increment is calculated based on the dilatation angle. The permeability is set to a second correlation value that is proportional to the cube of the sum of the fracture aperture and the shear dilatation aperture increment. During shear slip, due to the rough and uneven surface, the two fracture surfaces will climb each other, leading to an increase in fracture aperture. This phenomenon is called the shear dilatation effect. The shear dilatation aperture increment is proportional to the shear slip amount and the dilatation angle. The corrected fracture aperture is the initial aperture minus the normal closure amount plus the shear dilatation aperture increment, i.e.: ; in, The initial fracture aperture is given. Then, based on the corrected fracture aperture, the cubic law is used to calculate the fracture permeability under slip conditions.

[0060] When the contact state is closed, a residual aperture threshold is introduced, and the permeability is set to a third correlation value that is proportional to the cube of the residual aperture threshold. Even under very high normal effective stress, due to the roughness of the fracture surface and the interlocking of mineral grains, the fracture cannot be completely closed, and a small residual aperture will always exist. The magnitude of the residual aperture is related to the roughness of the fracture surface, the rock strength, and the normal effective stress. For closed fractures, the fracture surfaces are in close contact, and fluid can only seep through the residual aperture; therefore, the residual aperture is used to calculate permeability.

[0061] In this embodiment, the correction rules for crack opening and permeability under different contact states are shown in Table 1.

[0062] Table 1 Correction rules for fracture aperture and permeability under different contact states ; To address tangential shear stress, this embodiment establishes corresponding aperture and permeability correction rules for three different contact states of the fracture, fully considering the influence of normal closure, shear slip, and dilatation effects on the fracture seepage characteristics. Table 1 clearly illustrates the calculation methods under different conditions, facilitating understanding and implementation by those skilled in the art. This dynamic correction mechanism can accurately simulate the seepage evolution process of fractures under complex stress paths.

[0063] In a preferred embodiment, asynchronous time steps are set for the pore node control body and the fracture edge control body, respectively. When performing time difference on the discretized equations according to the asynchronous time steps, the first pore time step for the pore node control body and the first fracture time step for the fracture edge control body are adaptively determined based on the rate of change of fluid pressure and skeleton displacement. The rate of change of fluid pressure is the absolute value of the difference between the pore pressure at the current time and the previous time divided by the time step at the previous time. The rate of change of skeleton displacement is the absolute value of the difference between the skeleton displacement at the current time and the previous time divided by the time step at the previous time. For each pore node control body, the maximum value of its fluid pressure change rate and skeleton displacement change rate is calculated, and the next time step is determined based on this maximum value. The larger the rate of change, the smaller the time step; the smaller the rate of change, the larger the time step. The same method is used to determine the next time step for each fracture edge control body.

[0064] When the time step of the first pore is not synchronized with the time step of the first crack, an interpolation method is used on the control volume with a smaller time scale to obtain the physical quantities of the control volume with a larger time scale at the boundary moment. Assume the pore time step is... The time step of the crack is ,and Within one pore time step, the fracture control volume will undergo... This is the second update. During each fracture update, the pore pressure and skeleton displacement values ​​of the pore control volume at that moment need to be obtained. Since the pore control volume is only updated at the beginning and end of the pore time step, the physical quantity values ​​at intermediate moments need to be obtained through interpolation. Linear interpolation or quadratic interpolation is used, predicting the physical quantity values ​​at the current moment based on the physical quantity values ​​from the previous few moments.

[0065] Flux compensation calculations are performed on the control volumes with large time scales based on the physical quantities obtained through interpolation to ensure that the pore node control volume and the fracture edge control volume satisfy local conservation within their respective time steps. During each update of the fracture control volume, the fluid flux and momentum flux between pores and fractures are calculated. When the pore control volume is updated, the fluxes generated by all fracture updates within a pore time step are accumulated as the total flux of the pore control volume. This process ensures that the total flux flowing into and out of the pore control volume equals the change in physical quantities within the control volume over the entire pore time step, satisfying local mass and momentum conservation.

[0066] In this embodiment, the formula for calculating the adaptive time step is: ; in, For the next time step, For the current time step, This represents the difference in pore pressure between the current moment and the previous moment. This represents the difference in skeleton displacement between the current moment and the previous moment. and This represents the maximum permissible variation in pore pressure and skeleton displacement, with a value range of 10. -3 Up to 10 -1 .

[0067] The formula for linear interpolation is: ; in, For a moment The pore pressure value, For the first The starting point of a time step. For the first The endpoint of a time step. for The pore pressure value at time t. for The pore pressure value at a given time.

[0068] This embodiment employs an adaptive asynchronous time-step method, dynamically adjusting the time step based on the rate of change of physical quantities. Smaller time steps are used in regions of rapid physical quantity changes to ensure computational accuracy, while larger time steps are used in regions of gradual changes to improve computational efficiency. Interpolation and flux compensation methods guarantee local conservation under asynchronous time steps, avoiding non-physical losses of mass and momentum. This method effectively handles situations where there are significant differences in seepage rates between the pore matrix and the fracture network, improving the overall performance of the numerical simulation.

[0069] In a preferred embodiment, before acquiring reservoir geometry parameters and initial stress field, well logging data and seismic anisotropy data of the target reservoir are acquired, and shear wave transit time from the well logging data and azimuth anisotropy amplitude from the seismic anisotropy data are extracted. Well logging data is obtained during drilling using wireline logging tools, including sonic logging, density logging, and resistivity logging. Shear wave transit time is the time it takes for a shear wave to travel a unit distance through rock, reflecting the elastic properties of the rock. Seismic anisotropy data is obtained through 3D seismic exploration. When seismic waves pass through rock with directional fractures, shear wave splitting occurs, meaning the shear wave splits into two shear waves with mutually perpendicular polarization directions, and their propagation speeds and amplitudes differ. Azimuth anisotropy amplitude is the difference in seismic wave amplitude at different azimuth angles, reflecting the orientation and development degree of the fractures.

[0070] A shear wave splitting model is constructed, and the shear wave time difference and azimuth anisotropic amplitude are input into the model to invert and obtain the initial maximum and minimum horizontal principal stress directions of the target reservoir. The shear wave splitting model describes the propagation law of shear waves in anisotropic media, and its parameters include the fast shear wave direction, the slow shear wave direction, and the shear wave splitting time difference. The fast shear wave direction is parallel to the direction of the maximum horizontal principal stress, and the slow shear wave direction is parallel to the direction of the minimum horizontal principal stress. By fitting the shear wave time difference from well logging data and the azimuth anisotropic amplitude from seismic data, the parameters of the shear wave splitting model can be inverted, and thus the direction of the initial geostress can be obtained.

[0071] Boundary loads were applied to the geomechanical model based on the initial maximum and minimum horizontal principal stress directions to calibrate and obtain the initial stress field. The geomechanical model is a three-dimensional geological model including the reservoir and surrounding rocks. Based on the inverted stress directions, corresponding horizontal and vertical loads were applied to the model's boundaries. The vertical load was calculated from the weight of the overlying strata. The magnitude of the horizontal load was calibrated based on the regional tectonic stress background and in-situ stress test data. The geomechanical model was solved using the finite element method to obtain the initial stress field distribution throughout the entire reservoir.

[0072] In this embodiment, the theoretical model for shear wave splitting is based on Hudson's fractured medium theory. This theory treats rocks containing directional fractures as equivalent to transversely isotropic media with an elastic stiffness matrix... It can be represented as: ; in, The elastic stiffness matrix of the background rock. This represents the increment in the elastic stiffness matrix caused by the crack. The elements involved are related to the density, aperture, aspect ratio, and direction of the fracture. By measuring the transverse wave velocity in different polarization directions, the directional parameters of the fracture can be deduced, and thus the direction of the geostress can be obtained.

[0073] For directional, coin-shaped cracks, an explicit quantitative formula for calculating the elastic stiffness increment is constructed using Hudson's first-order approximation theory. Crack density is defined. : ; In the formula, per unit volume The number of internal cracks Define the fracture radius; define the fracture aspect ratio. : ; In the formula, let θ be the azimuth angle of the crack normal and φ be the dip angle. The tensor of the elastic stiffness increment caused by the crack. The explicit expression is: ; In the formula, , The values ​​are the Lamé constant and shear modulus of the background rock, respectively. , , , For the four Kronecker functions corresponding to different dimension indices in the fourth-order elastic stiffness tensor, ... The normalized crack factor is related to the crack aspect ratio; for cracks filled with fluid, the crack factor calculation formula is: ; By transforming the stiffness increment in the local coordinate system of the fracture to the global coordinate system using a coordinate transformation matrix, and substituting the fracture azimuth angle θ and dip angle φ, the values ​​of each element of the elastic stiffness increment matrix corresponding to any directional fracture can be obtained.

[0074] This embodiment constructs a complete forward simulation-inversion fitting closed-loop process, with the following specific steps: Step 1: Forward simulation calculation: Input the set of fracture parameters to be inverted, including fracture density, fracture aperture, fracture aspect ratio, fracture azimuth and dip angle. Calculate the equivalent elastic stiffness matrix of the rock according to the above elastic stiffness increment formula. Further solve for the fast shear wave velocity, slow shear wave velocity, shear wave splitting time difference corresponding to different polarization directions, and seismic wave amplitude response values ​​at different azimuth angles. Step 2: Constructing the objective function for fitting: Using the measured shear wave time difference from well logging and the measured azimuth anisotropic amplitude from seismic measurements as the fitting benchmarks, a weighted least squares method is used to construct the objective function for fitting. The expression is: ; In the formula, The time difference of shear wave splitting obtained from forward simulation. This refers to the time difference of shear wave splitting measured in well logging. The azimuth anisotropy amplitude difference is obtained from the forward simulation. This represents the azimuth anisotropy amplitude difference measured during earthquakes. , These are the fitting weighting coefficients for well logging data and seismic data, respectively. Step 3: Iterative Inversion Solution: With minimizing the objective function F as the optimization objective, a nonlinear optimization algorithm is used to iteratively update the fracture parameter set. When the objective function... When the value is less than the preset convergence threshold, the optimal set of fracture parameters is output. Based on the optimal fracture azimuth and dip angle, the directions of the initial maximum and minimum horizontal principal stresses are determined, and the geostress direction inversion is completed.

[0075] The initial stress field calibration employs a combination of trial and error and optimization algorithms. First, initial boundary load values ​​are set based on regional geological data, and the initial stress field is calculated. Then, the calculated stress values ​​are compared with measured values ​​from field stress test points, and the error between the two is calculated. By adjusting the boundary load values, the error is continuously reduced until it falls below a preset threshold. The optimization algorithm uses either a genetic algorithm or a particle swarm optimization algorithm to improve calibration efficiency and accuracy.

[0076] Taking the genetic algorithm as an example, the core parameter settings for initial stress field calibration are as follows: Inversion variables and their range: Inversion variables include the maximum horizontal principal stress value. Minimum horizontal principal stress value azimuth of maximum horizontal principal stress The value range is set based on the regional tectonic stress background. The value ranges from 0.8γH to 2.5γH. The value range is 0.6γ Up to 1.5γH, β ranges from 0° to 180°, where This represents the average unit weight of the overlying rock strata. The target reservoir depth; Constraints: The following basic mechanical constraints must be satisfied: ; The stress is the vertical principal stress, and the stress gradient between adjacent depths does not exceed 0.03 MPa / m, ensuring that the stress distribution conforms to the basic laws of rock mechanics. Fitness function: The fitness function is the reciprocal of the total objective error, which is the weighted sum of the simulated stress value and the error of the multi-source measured data; Multi-source data weighting rules: When fusing measured geostress data, well logging data, and seismic data, the weighting coefficients are set as follows: =0.5、 =0.3、 =0.2; where measured in-situ stress data, being direct data from field tests, is assigned the highest weight; well logging data, being continuous data at the single-well scale, is assigned the second highest weight; and seismic data, being indirect data at the regional scale, is assigned a lower weight; total error The expression is: ; In the formula This represents the relative error between measured and simulated values ​​of ground stress. The relative error of stress inverted from shear wave velocity in well logging. This represents the relative error in seismic anisotropic inversion stress.

[0077] This embodiment uses well logging data and seismic anisotropy data to invert the initial geostress direction, and then calibrates the initial stress field through geomechanical simulation, thereby improving the accuracy and reliability of the initial stress field. An accurate initial stress field is a prerequisite for fluid-structure interaction seepage evolution simulation and directly affects the calculation results of subsequent fracture mechanical behavior and seepage characteristics.

[0078] In a preferred embodiment, after outputting the seepage evolution results, the pore pressure field distribution and fracture velocity field distribution at the target time are extracted from the seepage evolution results. The target time is any time specified by the user, which can be the early, middle, or late stage of mining. The pore pressure field distribution consists of the pore pressure values ​​at all grid nodes throughout the entire reservoir. The fracture velocity field distribution consists of the fluid normal velocity values ​​at all fracture edge cells. The velocity values ​​are calculated according to Darcy's law and are proportional to the fracture permeability and pressure gradient.

[0079] Microscale fluid tracer testing was conducted in a physical model of the target reservoir to obtain tracer breakthrough curves. The physical model was either a core model or a large-scale physical simulation model created based on the geological characteristics of the target reservoir. Tracers were injected at the injection end of the model, and the tracer concentration over time was monitored at the production end to obtain the tracer breakthrough curves. The tracer breakthrough curves reflect the flow path and migration time of the fluid in the reservoir, containing information about reservoir heterogeneity and connectivity.

[0080] The fracture velocity field distribution is converted into a simulated breakthrough curve, and the dynamic error between the simulated breakthrough curve and the tracer breakthrough curve is calculated. Based on the fracture velocity field distribution, the tracer transport time in the fracture network is calculated. The tracer transport time from the injection end to the production end is equal to the path length divided by the velocity. By tracing all possible flow paths, the concentration of the tracer reaching the production end at different times is calculated, yielding the simulated breakthrough curve. The dynamic error is the sum of the squares of the concentration differences between the simulated and measured breakthrough curves at each time point.

[0081] If the dynamic error exceeds a preset threshold, it is added as a penalty term to the loss function of the graph neural network, and the parameters of the graph neural network are updated again until the dynamic error is less than the preset threshold. The weight coefficients of the penalty term are adjusted according to the magnitude of the error. The new loss function containing the penalty term is then substituted into the training process of the graph neural network, and the network parameters are updated through backpropagation. This process ensures that the network output not only satisfies the physical conservation laws but also matches the actual tracer test data, thus improving the accuracy of the simulation results.

[0082] In this embodiment, the tracer transport time is calculated using the particle tracking method. A large number of tracer particles are released at the injection end, and the velocity of each particle is equal to the fluid velocity at its location. Based on the velocity field distribution, the displacement of each particle within each time step is calculated. When a particle reaches the output end, its arrival time is recorded. By statistically analyzing the arrival times of all particles, the tracer concentration change curve over time is obtained, i.e., the simulated breakthrough curve.

[0083] The formula for calculating dynamic error is: ; in, For dynamic error, This represents the total number of time steps. for Simulated tracer concentration at time, for The measured tracer concentration at any given time.

[0084] The loss function including the penalty term is: ; in, For the new loss function, This represents the weighting coefficient of the penalty term.

[0085] This embodiment calibrates the numerical simulation results using microscale fluid tracer test data and incorporates dynamic error as a penalty term into the loss function, achieving a solution method that combines data-driven and physics-driven approaches. This method effectively corrects simulation errors caused by model parameter uncertainties and the complexity of geological conditions, improving the accuracy and reliability of seepage evolution simulation results.

[0086] In a preferred embodiment, the update frequency of the dynamic topology graph is consistent with the time step. After each time step, the normal and tangential displacements of the fractures are recalculated based on the current skeleton displacement field distribution, the fracture aperture and permeability are updated, and then the edge attributes and existence state of the dynamic topology graph are updated. The updated dynamic topology graph serves as the input to the graph neural network for the next time step. This real-time update mechanism ensures that the dynamic topology graph accurately reflects the connectivity state of pores and fractures at the current moment.

[0087] The parameters of a graph neural network are updated once at each time step. In the first time step, the network is trained using randomly initialized parameters. In subsequent time steps, the network parameters converged in the previous time step are used as the initial parameters. This parameter inheritance mechanism accelerates the network's convergence speed and improves computational efficiency. For time steps where the physical field changes relatively smoothly, the network only needs a few iterations to converge; for time steps where the physical field changes drastically, the network requires more iterations to converge.

[0088] In this embodiment, the number of network training iterations for different time steps is shown in Table 2.

[0089] Table 2. Number of network training iterations under different degrees of change in physical fields. ; This embodiment establishes a collaborative mechanism between dynamic topology graph updating and graph neural network training, ensuring the accuracy and efficiency of numerical simulation. Table 2 shows the number of network training iterations under different physical field variations, providing a reference for parameter settings in practical applications. This collaborative mechanism can adapt to changes in fluid-structure interaction intensity at different mining stages, achieving efficient and stable numerical simulation.

[0090] In a preferred embodiment, the meshing of the porous continuous medium and the fractured discrete medium employs an adaptive mesh refinement technique. During the simulation, the mesh density is dynamically adjusted based on the magnitude of the physical quantity gradient. In regions with large physical quantity gradients, such as pressure abrupt change zones, stress concentration zones, and fracture tip regions, mesh refinement is performed to improve computational accuracy; in regions with small physical quantity gradients, mesh coarsening is performed to improve computational efficiency. The adaptive mesh refinement technique can minimize the number of meshes while ensuring computational accuracy, thereby improving the overall performance of the numerical simulation.

[0091] The number of nodes and edges in a dynamic topology graph changes dynamically with the adaptive refinement and coarsening of the mesh. When the mesh is refined, corresponding nodes and edges are added; when the mesh is coarsened, corresponding nodes and edges are merged. Graph neural networks can automatically adapt to dynamically changing graph structures without requiring the reconstruction of the network model. This adaptive graph structure processing capability is one of the key advantages of graph neural networks compared to traditional numerical methods.

[0092] This embodiment employs a combination of adaptive mesh refinement and dynamic graph structure processing to further improve the accuracy and efficiency of numerical simulation. Adaptive mesh refinement effectively captures the local abrupt changes in the physical field, while dynamic graph structure processing ensures that the graph neural network can adapt to constantly changing mesh structures. This approach is particularly suitable for simulating physical processes with strong local nonlinear characteristics, such as crack propagation and connectivity.

Claims

1. A fluid-structure interaction simulation method for seepage evolution in pore-fracture type oil and gas reservoirs, characterized in that, include: Obtain reservoir geometric parameters and initial stress field, and establish the mapping relationship between pore continuous medium grid and fracture discrete medium grid; The normal displacement of the fracture is calculated based on the stress field distribution, the fracture aperture and permeability are updated, and a dynamic topology map characterizing the connectivity between pores and fractures is generated. The mass and momentum conservation equations in the fluid-structure interaction partial differential equations are discretized, and the pore pressure and skeleton displacement are characterized by the node features of the dynamic topology graph, while the fracture permeability and normal stiffness are characterized by the edge features. A graph neural network is constructed, and the dynamic topology graph is input into the graph neural network. The balance residual between fluid pressure and skeleton deformation is used as the loss function to solve the pore pressure field and stress field distribution end-to-end, and the seepage evolution results are output.

2. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 1, characterized in that, The process of establishing the mapping relationship between the porous continuous medium mesh and the fractured discrete medium mesh includes: calculating the shortest spatial distance from the center point of the porous continuous medium mesh to the surface of the fractured discrete medium mesh, and constructing a spatial proximity matrix; Extract the set of boundary nodes of the porous continuous medium mesh, and project the set of boundary nodes onto the normal direction of the fractured discrete medium mesh to obtain the set of projection intersection points; Based on the spatial proximity matrix and the set of projection intersections, equivalent weighting coefficients for cross-scale mass and momentum transfer are constructed, and a bidirectional mapping relationship between the porous continuous medium grid and the fractured discrete medium grid is established according to the equivalent weighting coefficients.

3. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 1, characterized in that, The normal displacement of the fracture is calculated based on the stress field distribution, the fracture aperture and permeability are updated, and a dynamic topology map characterizing the connectivity between pores and fractures is generated. This includes: extracting the effective normal stress and tangential shear stress of the fracture discrete medium mesh based on the stress field distribution, and calculating the normal closure and shear slip of the fracture discrete medium mesh. Based on the normal closure amount and the shear slip amount, the Mohr-Coulomb criterion is introduced to determine the contact state of the fractured discrete medium mesh, and the contact state includes fully open, slip, and locked. The crack aperture and permeability are dynamically corrected based on the contact state, and the corrected crack aperture and permeability are assigned to the edge attributes of the topology graph to generate the dynamic topology graph.

4. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 1, characterized in that, The discretization of the mass conservation and momentum conservation equations in the fluid-structure interaction partial differential equations, and the characterization of pore pressure and skeleton displacement by the node features of the dynamic topology graph, and the characterization of fracture permeability and normal stiffness by the edge features, includes: using the control volume finite element method to spatially discretize the mass conservation equations and the momentum conservation equations to obtain the pore node control volume and the fracture edge control volume. Asynchronous time steps are set for the pore node control body and the crack edge control body respectively, and time difference is performed on the discretized equations based on the asynchronous time steps; The differential pore pressure and skeleton displacement variables are mapped to the node feature vectors of the dynamic topology graph, and the differential fracture permeability and normal stiffness variables are mapped to the edge feature vectors of the dynamic topology graph.

5. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 4, characterized in that, The construction of the graph neural network involves inputting the dynamic topology graph into the graph neural network, using the balance residual between fluid pressure and skeleton deformation as the loss function, solving the pore pressure field and stress field distribution end-to-end, and outputting the seepage evolution results. This includes: in the graph neural network, calculating the message passing weight between adjacent nodes based on the edge feature vectors of the dynamic topology graph, performing weighted aggregation and nonlinear transformation on the node feature vectors, and obtaining the hidden layer features of the nodes. The hidden layer features of the nodes are residually connected with the original node feature vectors to output the predicted pore pressure and skeleton displacement values. The equilibrium residuals of the predicted pore pressure and the predicted skeleton displacement with respect to the mass conservation equation and the momentum conservation equation are calculated. A loss function is constructed using the norm of the equilibrium residuals. The parameters of the graph neural network are updated by backpropagation until the loss function converges, and the seepage evolution results are output.

6. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 2, characterized in that, Projecting the set of boundary nodes onto the normal direction of the fractured discrete medium mesh to obtain the set of projection intersection points includes: calculating the projection line of each boundary node in the set of boundary nodes along the normal direction of the fractured discrete medium mesh; Determine the intersection state between the projected line and the fractured discrete medium mesh surface. If there are multiple intersection points, select the intersection point closest to the boundary node as the valid intersection point. If the projected line does not intersect the fractured discrete medium mesh surface, then the nearest face point of the boundary node on the fractured discrete medium mesh surface is searched as a substitute intersection point. Aggregate all the valid intersection points and the alternative intersection points to construct the projected intersection point set.

7. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 3, characterized in that, The method of dynamically correcting the fissure aperture and permeability based on the contact state includes: when the contact state is fully open, setting the permeability to a first correlation value that is proportional to the cube of the fissure aperture according to the cubic law. When the contact state is slippage, the shear dilatation increment is calculated based on the shear dilatation angle, and the permeability is set to a second correlation value that is proportional to the cube of the sum of the fracture aperture and the shear dilatation increment. When the contact state is locked, a residual opening threshold is introduced, and the permeability is set to a third correlation value that is proportional to the cube of the residual opening threshold. The fracture opening and permeability are dynamically corrected using the first correlation value, the second correlation value, and the third correlation value.

8. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 4, characterized in that, Asynchronous time steps are set for the pore node control body and the fracture edge control body respectively. Time difference is performed on the discretized equations based on the asynchronous time steps, including: adaptively determining the first pore time step of the pore node control body and the first fracture time step of the fracture edge control body based on the rate of change of fluid pressure and skeleton displacement. When the time step of the first pore is not synchronized with the time step of the first crack, an interpolation method is used on the control volume with a small time scale to obtain the physical quantity of the control volume with a large time scale at the boundary moment. Based on the physical quantities obtained by the interpolation, flux compensation calculations are performed on the control volume with a large time scale so that the pore node control volume and the crack edge control volume satisfy local conservation within their respective time steps.

9. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 1, characterized in that, Before acquiring the reservoir geometric parameters and initial stress field, the method further includes: acquiring well logging data and seismic anisotropy data of the target reservoir, and extracting the shear wave time difference from the well logging data and the azimuth anisotropy amplitude from the seismic anisotropy data. A shear wave splitting model is constructed, and the shear wave time difference and the azimuth anisotropic amplitude are input into the shear wave splitting model to invert and obtain the initial maximum horizontal principal stress direction and the initial minimum horizontal principal stress direction of the target reservoir. Boundary loads are applied to the geomechanical model based on the initial maximum horizontal principal stress direction and the initial minimum horizontal principal stress direction, and the initial stress field is obtained by calibration.

10. The fluid-structure interaction seepage evolution simulation method for pore-fracture type oil and gas reservoirs according to claim 9, characterized in that, After outputting the seepage evolution results, the method further includes: extracting the pore pressure field distribution and fracture velocity field distribution at the target time from the seepage evolution results; Microscale fluid tracer testing was performed on the physical model of the target reservoir to obtain tracer breakthrough curves; The fracture velocity field distribution is converted into a simulated breakthrough curve, and the dynamic error between the simulated breakthrough curve and the tracer breakthrough curve is calculated. If the dynamic error exceeds a preset threshold, the dynamic error is added as a penalty to the loss function of the graph neural network, and the parameters of the graph neural network are updated again until the dynamic error is less than the preset threshold.

Images