A cavern group-aquifer non-matching grid coupling modeling method

Abstract

The application discloses a cavern group-aquifer non-matching grid coupling modeling method and belongs to the technical field of three-dimensional geological modeling. The method comprises the following steps: collecting multi-source detection data to construct a unified three-dimensional geological framework; generating a fine non-structured tetrahedral grid on the side of a cavern and a coarse structured hexahedral grid on the side of an aquifer to form a non-matching double-grid system; creating a Mortar intermediate interface on a hydraulic contact boundary; constructing two-stage projection operators to realize the conservative transfer of physical quantities on the non-matching interface; establishing a seepage-stress coupling saddle point system with the Mortar multiplier as a constraint and synchronously solving; and realizing model dynamic updating and uncertainty quantification based on the Bayes theory and ensemble Kalman filtering. The application solves the problems of non-conservative transfer of physical quantities between non-matching grids, inability of strong coupling and difficulty in dynamic updating of the model, and significantly improves the precision and efficiency of seepage-stress coupling simulation of large underground cavern groups.

Description

Technical Field

[0001] This invention relates to the field of three-dimensional geological modeling, and in particular to a method for modeling cavern groups and aquifers with mismatched meshes. Background Technology

[0002] During the construction and operation phases of underground cavern complex projects, the groundwater seepage field, stress field, and deformation field between the surrounding rock and aquifer are highly coupled. Accurately constructing an integrated multi-field coupled model of the cavern complex and aquifer is crucial for early warning of water inrush risks, seepage prevention design, and surrounding rock stability analysis. Currently, the related technologies for multi-field coupled modeling of underground cavern surrounding rock and aquifer have significant limitations, specifically as follows: Chinese patent "Quantitative Calculation Method and System for Underground Surrounding Rock Damage Degree Coupled by Multiphysics Fields" (Publication No.: CN121723936A) proposes a quantitative analysis method for underground surrounding rock damage through bidirectional coupling of stress field and seepage field. By constructing a multiphysics field coupled numerical model, it achieves joint solution of stress and seepage, and completes quantitative calculation of surrounding rock damage and output of cloud maps based on damage potential and damage degree functions. This method can consider the influence of multi-field coupling on surrounding rock stability, improve the accuracy of damage evaluation, and solve the problem of insufficient consideration by traditional single-physics field analysis mechanisms. However, this method uses a globally unified mesh model, which cannot simultaneously consider the computational efficiency of the fine local structure of the cavern and the large area of ​​the aquifer; interface data transfer relies on node matching and does not support non-matched mesh coupling; it does not set conservation constraints for physical quantities such as flow rate and pressure, which easily leads to seepage field distortion; and the model only supports static analysis and cannot achieve dynamic updates and adaptive corrections based on construction exposure and monitoring data, making it difficult to meet the dynamic coupling requirements of cavern groups and aquifers.

[0003] Chinese patent "A Fluid-Structure Coupling Simulation Method and System for Water Immersion Weakening of Tunnel Surrounding Rock" (Application No.: CN202510657478.X) considers the water immersion weakening effect, establishes an equivalent variation evolution equation for water immersion, and constructs a seepage stress-damage coupled elastoplastic constitutive model, which can reflect the deteriorating influence of groundwater weakening on the mechanical and seepage parameters of the surrounding rock. This method can improve the rationality of the stability analysis of water-rich tunnel surrounding rock and is more in line with the actual engineering situation under long-term seepage. However, it uses a traditional finite element unified mesh, which does not support non-matched mesh coupling; interface data is only transmitted through interpolation, without conservation constraints; the model focuses on constitutive and damage, and does not achieve cross-scale integrated modeling of the aquifer in the tunnel, and cannot meet the needs of local fine detail and efficient global calculation.

[0004] Chinese patent "A Dynamic Prediction Method and System for Tunnel Lining Stress Based on GMSMODFLOW and ABAQUS" (Publication No.: CN119962329A) uses groundwater simulation software and finite element software to solve the problem, mapping the seepage field calculation results to a mechanical model to achieve dynamic prediction of tunnel lining stress. This method can achieve indirect coupling between the seepage field and the stress field, improving the rationality of the lining stress analysis. However, it uses unidirectional data transfer between software, resulting in loose coupling rather than strong coupling, and lacks real-time bidirectional feedback between physical fields; the interface data is mapped only through simple interpolation, failing to satisfy flow conservation and pressure balance; the model uses a uniform mesh, which cannot handle the problem of mismatched mesh coupling between the cavern and the aquifer; and it lacks adaptive dynamic updating capabilities, making it difficult to adapt to real-time model correction and dynamic risk assessment during the construction of underground cavern groups.

[0005] In summary, existing multi-field coupling modeling techniques generally suffer from drawbacks such as non-conservation of interface physical quantities, inability to couple with mismatched meshes, and difficulty in dynamically updating the model during construction and monitoring. Therefore, there is an urgent need for an integrated modeling method for cavern groups and aquifers that can achieve high precision, high efficiency, strong conservation coupling under mismatched mesh conditions, and support dynamic evolution. Summary of the Invention

[0006] The purpose of this invention is to provide a non-matching mesh coupling modeling method for cavern groups and aquifers, addressing the shortcomings of existing technologies such as the inability to couple non-matching interfaces, non-conservation of interface flow, discontinuous pressure, difficulty in balancing computational efficiency and accuracy, and difficulty in dynamically updating the model.

[0007] To achieve the above objectives, the present invention adopts the following technical solution: A method for modeling cavern group-aquifer mismatched mesh coupling includes the following steps: S1. Collect multi-source detection data, perform coordinate unification, spatiotemporal registration, standardization and noise removal processing on the multi-source detection data, and construct a unified three-dimensional geological framework model including strata, cavern groups, aquifers, water-rich zones and water-conducting channels; S2. Within the unified three-dimensional geological framework model, unstructured tetrahedral fine meshes are generated for the cavern group and the area near the surrounding rock, and structured hexahedral coarse meshes are generated for the far-field aquifer area, thus constructing a non-matching dual-mesh system where nodes do not coincide and elements do not correspond at the contact interface. S3. At the hydraulic contact boundary of the mismatched dual-grid system, an independent Mortar intermediate interface is created, which serves as a virtual transition layer for the conservation and transfer of physical quantities between the mismatched grids. S4, Construct the first stage L 2An orthogonal projection operator projects the physical quantities in the fine grid of the cavern onto the Mortar intermediate interface. Then, a second-level dual basis function projection operator is constructed to distribute the conservation information on the Mortar intermediate interface to the coarse grid of the aquifer, thereby realizing the conservation and transfer of physical quantities at the mismatched interface. S5. Using Mortar multipliers as interface constraints, a hybrid variational equation is established to form a seepage-stress coupled saddle point system. The groundwater seepage control equation, the effective stress balance equation of the surrounding rock, and the interface flux continuity equation are solved simultaneously to realize real-time bidirectional exchange and global coupling of the physical field between the cavern and the aquifer. S6. When construction reveals, supplements exploration, or monitoring data is updated, the surrounding rock parameters of the tunnel side are updated based on Bayesian theory, and the first-level L is corrected. 2 The orthogonal projection operator uses ensemble Kalman filtering to invert the aquifer parameter field and simultaneously optimizes the second-level dual basis function projection operator to achieve dynamic updating of all elements of the model and quantification of uncertainty.

[0008] A storage medium that stores instructions and data for implementing a cavern group-aquifer mismatched mesh coupling modeling method.

[0009] A cavern group-aquifer mismatched mesh coupling modeling device includes: a processor and a storage medium; the processor loads and executes instructions and data in the storage medium to implement a cavern group-aquifer mismatched mesh coupling modeling method.

[0010] The beneficial effects of this invention are as follows: 1. Seamless strong coupling of non-matching meshes is achieved: by constructing a Mortar intermediate interface and a two-level projection operator (L... 2 (Orthogonal projection + dual basis function projection) For the first time, the transfer of physical quantities that satisfy normal flow conservation, pressure continuity and energy conservation was realized between the fine grid of the cavern group and the coarse grid of the aquifer, overcoming the problem of lack of conservation caused by the traditional method relying on node matching and only being able to perform interpolation mapping; 2. Balancing computational accuracy and efficiency: A partitioned mesh discretization strategy is adopted. Fine tetrahedral meshes of 0.5-2m are used in the cavern and near-surrounding rock areas to ensure geometric and physical field resolution, while coarse hexahedral meshes of 5-50m are used in the far-field aquifer area to reduce computational load. In typical engineering cases, this can reduce the total degrees of freedom by more than 60%, while the reconstruction error of the projection operator constant field is less than 10. -12 ; 3. Supports real-time bidirectional coupling of seepage and stress: The established hybrid variational equation and saddle point system simultaneously solve the groundwater seepage control equation, the effective stress balance equation of the surrounding rock, and the interface flux continuity equation. The Uzawa iteration that satisfies the LBB condition is adopted to achieve weak continuity and global conservation of interface pressure and flow rate, avoiding physical field distortion caused by loose coupling or unidirectional mapping. 4. Possesses dynamic model updating and uncertainty quantification capabilities: Based on Bayesian theory, the cavern side parameters are updated and the projection operator is corrected. An ensemble Kalman filter is used to invert the aquifer permeability coefficient field and pore water pressure field, realizing the synchronous evolution of all elements of the grid, Mortar interface, projection operator and coupling field, and outputting parameter confidence intervals and risk levels, filling the gap in existing methods that cannot adaptively correct based on construction exposure data. 5. Strong engineering applicability: It can be widely used in the dynamic prediction of water inflow, surrounding rock stability analysis and seepage prevention design during the construction period of deep underground cavern groups such as water conservancy and hydropower, transportation tunnels, and mining, providing an integrated numerical simulation tool for early warning of water inrush risk. Attached Figure Description

[0011] Figure 1 This is a schematic diagram of the method flow of the present invention; Figure 2 This is a schematic diagram of a non-matching mesh system and the Mortar interface; Figure 3 This is a schematic diagram of the principle of the dual-projection operator; Figure 4 This is a schematic diagram of the hardware device operation according to an embodiment of the present invention. Detailed Implementation

[0012] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific examples described herein are merely illustrative and not intended to limit the scope of the invention.

[0013] Before formally describing the present invention, a general description of the solution of the present invention will be given first to facilitate understanding.

[0014] Example 1: Please refer to Figure 1 The present invention provides a method for modeling a cavern group-aquifer mismatched mesh coupling, comprising the following steps: S1. Collect multi-source detection data, perform coordinate unification, spatiotemporal registration, standardization and noise removal processing on the multi-source detection data, and construct a unified three-dimensional geological framework model including strata, cavern groups, aquifers, water-rich zones and water-conducting channels; It should be noted that the multi-source detection data collected in step S1 includes ground-penetrating radar data, borehole test data, hydrological test data, and groundwater monitoring data.

[0015] Specifically, this invention collects multi-source data such as ground-penetrating radar, borehole testing, hydrological experiments, and groundwater monitoring, and completes coordinate unification, spatiotemporal registration, data standardization, and noise removal. It constructs a unified three-dimensional geological framework model that includes strata, cavern groups, aquifers, water-rich zones, and water-conducting channels, providing a unified geometric benchmark and parameter basis for subsequent generation of partitioned grids, and ensuring spatial consistency and semantic integrity of multi-source information.

[0016] S2. Within the unified three-dimensional geological framework model, unstructured tetrahedral fine meshes are generated for the cavern group and the area near the surrounding rock, and structured hexahedral coarse meshes are generated for the far-field aquifer area, thus constructing a non-matching dual-mesh system where nodes do not coincide and elements do not correspond at the contact interface. It should be noted that in step S2, the unit size of the unstructured tetrahedral fine mesh at the cavern outline is controlled within a first preset range, and the unit size of the structured hexahedral coarse mesh is controlled within a second preset range. The two types of meshes form a natural mismatched interface at the hydraulic contact interface.

[0017] Preferably, in this invention, the first preset range is 0.5-2m; the second preset range is 5-50m.

[0018] In other words, within a unified three-dimensional geological framework, partitioned mesh discretization is performed. Unstructured tetrahedral fine meshes are generated for the cavern complex and near-surrounding rock areas, with element sizes controlled between 0.5 and 2 m, accurately reproducing the cavern outline, lining, fissures, and local aquifer structures. Structured hexahedral coarse meshes are generated for the far-field aquifer areas, with element sizes between 5 and 50 m, ensuring computational efficiency over a large area. The nodes of the two types of meshes do not coincide and the elements do not correspond at the hydraulic contact interface, forming a natural mismatched interface.

[0019] S3. At the hydraulic contact boundary of the mismatched dual-grid system, an independent Mortar intermediate interface is created, which serves as a virtual transition layer for the conservation and transfer of physical quantities between the mismatched grids. Please refer to Figure 2 It should be noted that step S3, creating the Mortar intermediate interface, further includes: S31. Extract the intersection of the cavern group region boundary and the aquifer region boundary from the non-matching dual grid system, mark the cavern group side interface as the non-Mortar side, and mark the aquifer side interface as the Mortar side. Specifically, the spatial interface between the cavern complex and the aquifer is precisely extracted from the grids of both regions, and Boolean operations are used to calculate the intersection of the two boundaries. ,in This marks the boundary of the cavern complex area. The aquifer boundary was defined, and the cavern group side interface was then marked as the non-Mortar side, while the aquifer side interface was marked as the Mortar side.

[0020] S32. Map the intersection in three-dimensional space to two-dimensional parameter space, calculate the curvature of the intersection interface, and perform local encryption identification in the high curvature region. Specifically, the interface in three-dimensional space Mapping to a two-dimensional parameter space (u, v) facilitates subsequent mesh generation. For each interface surface patch, a local coordinate system is established, and a harmonic mapping is used to flatten the surface to the two-dimensional parameter domain, ensuring the bijectivity of the mapping and avoiding overlap and flipping. Then, the Gaussian curvature and average curvature of the interface are calculated, and high curvature regions (such as cavern corners and densely fractured zones) are extracted.

[0021] The high curvature region is defined as having a Gaussian curvature absolute value greater than 1. Or the average rate of change of curvature is greater than The area. The cell size after local refinement is taken as the original Mortar interface mesh size. of to Specifically, when the curvature exceeds twice the threshold, take... When it exceeds 1 times the threshold, take The remaining areas retain their original size.

[0022] S33. Generate a mesh in the two-dimensional parameter space. The cell size is calculated as the geometric mean of the mesh sizes on both sides, and Delaunay triangulation is used to generate the surface mesh. Specifically, a structured or semi-structured mesh is generated in the parameter domain (u, v), and the cell size is calculated using geometric mean. ,in For the size of the cavern grid, The surface mesh is generated using Delaunay triangulation, which is the mesh size for the aquifer.

[0023] S34. Establish the mathematical mapping relationship between the Mortar intermediate interface and the grids on both sides, define the basis functions of the cavern group side, the Mortar interface side and the aquifer side respectively, and calculate the Mortar interface quality matrix. Specifically, to establish the mathematical mapping relationship between the Mortar interface and the meshes on both sides, the basis functions for each region are first defined, with linear triangular element basis functions used for the cavern group survey: , specifically, ,in , , , All coefficients are determined by the node coordinates, through... Confirmed; the Mortar interface uses piecewise linear basis functions: , where this function satisfies To ensure the conservation of projection; bilinear quadrilateral basis functions are used on the aquifer side: , specifically, In the parameter range [-1, 1] 2 Define the above; then calculate the Mortar interface quality matrix: , ,in This represents the number of nodes in the Mortar interface.

[0024] S35. Establish the adjacency relationship between the Mortar middle interface and the meshes on both sides, adopt the region decomposition strategy to support parallel computing, and establish a version control mechanism to automatically trigger interface updates when the mesh is reconstructed. Specifically, an adjacency relationship is established between the Mortar interface and the meshes on both sides to support subsequent parallel computation and dynamic updates. First, for each Mortar interface cell, a list of adjacent cavern group side cell numbers and aquifer side cell numbers is recorded. Then, a domain decomposition strategy is adopted to assign the Mortar interface to processors in adjacent subdomains. In the MPI parallel environment, Mortar interface data is redundantly stored between adjacent processors to reduce communication overhead. Finally, a version control mechanism is established to record the generation of Mortar meshes and projection matrices. When the cavern group side mesh is reconstructed, the updates of Mortar interfaces and projection operators are automatically triggered.

[0025] S36. Perform a closure check and external normal direction unification on the Mortar intermediate interface, verify that the interface area matches the interface area of ​​the meshes on both sides, and achieve geometric consistency and area conservation.

[0026] Specifically, to ensure the geometric consistency of the Mortar interface, a closure check is performed on the Mortar interface to ensure that it forms a closed surface without cracks or overlaps. Simultaneously, the outward normal direction of the Mortar interface (pointing towards the aquifer side) is standardized and verified. The area is matched with the interface area of ​​the grids on both sides, thus achieving area conservation.

[0027] S4, Construct the first stage L 2 An orthogonal projection operator projects the physical quantities in the fine grid of the cavern onto the Mortar intermediate interface. Then, a second-level dual basis function projection operator is constructed to distribute the conservation information on the Mortar intermediate interface to the coarse grid of the aquifer, thereby realizing the conservation and transfer of physical quantities at the mismatched interface. Please refer to Figure 3 It should be noted that constructing the dual projection operator in step S4 further includes: S41. Verify the spatial overlap between the Mortar intermediate interface and the boundary grids on both sides, establish a local coordinate system and calculate the Jacobian determinant; Specifically, based on the Mortar interface mesh generated in step S3, its spatial overlap with the boundary meshes of the cavern group side and the aquifer side is verified. A local coordinate system is established and the Jacobian determinant is calculated to ensure the accuracy of the area element transformation.

[0028] S42. Define the linear triangular basis functions of the cavern group side boundary element, the piecewise linear basis functions of the Mortar interface element, and the bilinear quadrilateral basis functions of the aquifer side boundary element, and configure Gauss-Legend integration points on the Mortar interface element. Specifically, based on step S34, linear triangular basis functions are defined for the side boundary elements of the cavern group. Define piecewise linear basis functions for Mortar interface elements. Define bilinear quadrilateral basis functions for the aquifer side boundary elements. Configure Gauss-Legend integration points on the Mortar interface unit to establish a mapping relationship between physical coordinates and integration weights.

[0029] S43. Calculate the first projection operator from the cavern group side to the Mortar interface, construct the dual basis function that satisfies the biorthogonality condition, and then calculate the second projection operator from the Mortar interface to the aquifer side. Specifically, the first projection operator (cavity group → Mortar interface) is calculated: , , ,in The number of side boundary nodes of the cavern group, in passing Obtain Mortar interface coefficients Then, dual basis functions satisfying the bioorthogonality condition are constructed. To ensure the conservation of the second projection matrix, specifically, When the bioorthogonality condition is satisfied: ; Finally, the second projection operator (Mortar interface → aquifer) is calculated: First, the original projection matrix is ​​calculated: ,in The number of boundary nodes on the aquifer side is used to obtain the second projection matrix by correcting it using the dual basis function. Finally, through the projection formula The pressure values ​​at the aquifer side nodes were obtained.

[0030] S44. Perform constant field accuracy and global conservation checks on the first and second projection operators, check the condition number of the mass matrix, and output the projection operator that satisfies the conservation constraints.

[0031] Specifically, the projection operator is validated by performing constant field accuracy and global conservation checks to verify its ability to accurately reconstruct the constant field and the conservation of the total integral. The condition number of the mass matrix is ​​checked to ensure numerical stability, and the final output is the projection of the side pressures of the cavern group. Aquifer lateral pressure projection .

[0032] S5. Using Mortar multipliers as interface constraints, a hybrid variational equation is established to form a seepage-stress coupled saddle point system. The groundwater seepage control equation, the effective stress balance equation of the surrounding rock, and the interface flux continuity equation are solved simultaneously to realize real-time bidirectional exchange and global coupling of the physical field between the cavern and the aquifer. It should be noted that step S5, which involves constructing and solving the multi-field coupled system, further includes: S51. Construct the system matrices for each physical field, including the side stiffness matrix of the cavern group, the side stiffness matrix of the aquifer, the Biot coupling matrix, the side seepage system matrix of the cavern group, and the side seepage system matrix of the aquifer. Specifically, the system matrices for each physical field are constructed, including the stress field matrix, Biot coupling matrix, and seepage field matrix. lateral stiffness matrix of cavern group lateral stiffness matrix of aquifer Where B is a geometric matrix, satisfying D is the elasticity matrix; The Biot coupling matrix Q achieves the coupling of seepage to stress, converting pore pressure into equivalent nodal forces. Specifically, the Biot coupling matrix on the cavern group side... Biot coupling matrix on the aquifer side ,in Here, m is the Biot coefficient, and m is the unit normal vector matrix. The pressure shape function matrix is ​​used to interpolate the pressure at each node to the integration point; The seepage field matrix includes: the seepage system matrix of the cavern group. Aquifer lateral seepage system matrix Where H is the penetration matrix (transmission matrix), satisfying R is the water storage matrix (capacity matrix), which satisfies , Permeability coefficient (m / s) represents the ability of a rock mass to allow water to flow through it; The specific water storage coefficient (m) -1This represents the amount of water released due to a unit change in head. Additionally, through... To achieve the consolidation effect, among which , This represents the volumetric strain rate within the current time step.

[0033] Among the parameters mentioned above, the geometric matrix Satisfy the strain-displacement relationship Elasticity matrix By Young's modulus Compared to Poisson Determined; Biot coefficient The range of values ​​is Complete rock samples Fractured rock mass Specific water storage coefficient Units are The value range is generally 100%. It can be obtained by inversion from the pressure water test.

[0034] S52. Based on the system matrix, a complete seepage-stress coupled saddle point system is assembled, wherein the traction force Mortar multiplier and the flow rate Mortar multiplier are introduced as Lagrange multipliers for displacement constraint and pressure constraint, respectively. Specifically, a complete saddle point system is assembled based on the matrix formed in step S51:

[0035] Specifically, the mechanical equilibrium equations for the cavern group are:

[0036] Specifically, this can be expressed as: internal stress balance of the cavern group + volume force generated by pore pressure + traction force at the Mortar interface = external load; Aquifer mechanical equilibrium equations:

[0037] Specifically, this can be expressed as: internal stress balance of the aquifer + volume force generated by pore pressure - traction force at the Mortar interface = external load; The conservation equation for seepage flow in a group of caverns is as follows:

[0038] Specifically, it can be expressed as: volumetric strain change (consolidation) + pressure diffusion + Mortar interface inflow = consolidation source term (historical volumetric strain change). Aquifer seepage conservation equation:

[0039] in The negative sign is because, from the Mortar interface, the flow rate into the cavern group is equal to the flow rate out of the aquifer, thus satisfying the continuity requirement. Mortar displacement constraint equations:

[0040] The displacements on both sides should be equal after calculation by the projection operator. At this point, there is no relative displacement at the interface, which satisfies the displacement compatibility condition.

[0041] Mortar pressure constraint equation:

[0042] in, For the stiffness of the cavern group, For aquifer stiffness, This is the transpose of the Biot coupling matrix on the side of the cavern group. This is the transpose of the Biot coupling matrix on the aquifer side. This is the transpose of the side displacement projection operator for the cavern group. This is the transpose of the displacement projection operator on the water-bearing side. This is the transpose of the side pressure projection operator for the cavern group. This is the transpose of the projection operator for the aquifer side pressure. The matrix of the side seepage system of the cavern group. The matrix represents the seepage system on the aquifer side. Let be the displacement vector of the cavern group nodes. Let be the displacement vector of the water-bearing side node. For the water pressure at the nodes of the cavern group, The water pressure at the aquifer node. For traction force Mortar multipliers, For flow Mortar multipliers.

[0043] It should be noted that the displacement projection operator and the pressure projection operator are constructed in a completely similar manner; the only difference lies in the physical quantity being projected. Specifically: First-level displacement projection operator (cavity fine mesh → Mortar interface): in: , The displacement shape function of the cavern side (same as the pressure shape function, using linear triangular elements) is projected to obtain the displacement field on the Mortar interface.

[0044] Second-level displacement projection operator (Mortar interface → aquifer rough mesh): in: , The aquifer lateral displacement shape function (bilinear quadrilateral element) is used to achieve weak continuity constraint on the interface displacement. .

[0045] S53. Based on the LBB condition, verify the stability of the Mortar interface constraint matrix. Use Uzawa iteration for global coupling calculation. Use the Mortar multiplier update equation to make the residual difference of physical quantities projected on both sides of the interface approach zero. It should be noted that the convergence condition of the Uzawa iteration in step S53 is: the displacement difference and pressure difference after projection on both sides of the interface are used as residuals, and convergence is determined when the residuals approach zero; the Mortar multiplier update equation is: , in and These are the Mortar multipliers for traction force and the Mortar multiplier for flow rate, respectively. The Mortar interface quality matrix. and These are the displacement residual and pressure residual after the m-th iteration, respectively. and is the iterative relaxation factor.

[0046] Specifically, based on the saddle point system in step S52, the Mortar interface constraint matrix is ​​verified according to the LBB condition requirements. That is, there exists a constant that is independent of the mesh size. Such that for all non-zero Mortar multipliers For each, a corresponding physical field variable (u, p) can be found that satisfies and In actual calculations, this is equivalent to requiring the projection operator constructed in step 3. , , , The row rank is equal to the number of Mortar nodes, ensuring that the interface constraint equations are not over-constrained or under-constrained; then, Uzawa iteration is used to calculate the residuals (the difference in physical quantities after projection on both sides). After the equation approaches zero, check for convergence and then perform the Mortar multiplier update equation: , Ultimately, this achieves weak continuity and global conservation of interface pressure and traffic.

[0047] S54. Extract the water pressure, flow rate and traction force distribution on the Mortar interface, calculate engineering indicators, and realize real-time bidirectional exchange and global coupling of the physical field between the cavern and the aquifer.

[0048] Finally, the water pressure, flow rate, and traction force distribution on the Mortar interface are extracted, and engineering indicators such as the total inflow of the cavern group, the effective stress around the cavern, and the thickness of the water seal layer are calculated to realize real-time bidirectional exchange and global coupling of water pressure, seepage flow, and effective stress between the cavern and the aquifer.

[0049] S6. When construction reveals, supplements exploration, or monitoring data is updated, the surrounding rock parameters of the tunnel side are updated based on Bayesian theory, and the first-level L is corrected. 2 The orthogonal projection operator uses ensemble Kalman filtering to invert the aquifer parameter field and simultaneously optimizes the second-level dual basis function projection operator to achieve dynamic updating of all elements of the model and quantification of uncertainty.

[0050] It should be noted that in step S6, the dynamic update of all elements of the model includes: updating the surrounding rock parameters of the cavern side and correcting the first projection operator, inverting the aquifer permeability coefficient field and pore water pressure field and simultaneously optimizing the second projection operator, and finally realizing the synchronous update of the grid, Mortar interface, projection operator and coupling field, and outputting the parameter confidence interval, field distribution probability and risk level.

[0051] The specific method of updating the surrounding rock parameters of the tunnel based on Bayesian theory is as follows: using the fracture opening and filling status revealed during construction as observation data. The permeability coefficient of the cavern side As the parameter to be estimated, the prior distribution is taken as a log-normal distribution. Likelihood function Based on the local analytical solution of seepage, the posterior distribution is obtained by Markov chain Monte Carlo (MCMC) sampling. The posterior mean is taken as the updated parameter value, and the first projection operator is recalculated accordingly. .

[0052] The specific method for retrieving aquifer parameter fields using ensemble Kalman filtering involves constructing a system containing... aquifer permeability field of each set member and pore water pressure field The prediction step is obtained by advancing one time step based on the seepage equation. Analysis steps utilize monitoring water levels Update each member: in Let be the covariance matrix of permeability coefficient and water pressure. The covariance matrix for water pressure prediction. To determine the observation error covariance, the ensemble mean is taken as the inversion result, and the second projection operator is optimized accordingly. .

[0053] Example 2: This embodiment uses the dynamic prediction of water inflow and the analysis of surrounding rock stability during the construction of the underground powerhouse cavern group of a deep-buried hydropower station.

[0054] This embodiment focuses on a deep-buried hydroelectric power station's underground powerhouse cavern complex. The project comprises three main caverns: the main powerhouse, the main transformer room, and the tailrace surge chamber, along with connecting tunnels. The maximum burial depth is approximately 600 m. The aquifer in the area is fractured rock mass, with the initial groundwater level approximately 200 m above the cavern arches. Construction employs layered excavation. This embodiment applies the method of this invention to dynamically predict the cavern water inflow, pore water pressure distribution in the surrounding rock, and effective stress evolution throughout the construction process. The model is dynamically updated based on statistical data of the fractures revealed during excavation and borehole water level monitoring data.

[0055] Data from borehole columnar sections, borehole pressure test data, ground-penetrating radar profiling, and long-term water level observation wells arranged during construction were collected from the engineering geological survey report. The coordinates of these data were unified to the engineering-independent coordinate system, and the time was registered to the construction progress timeline. Abnormal pressure values ​​from borehole tests were removed. A three-dimensional geological framework model was constructed, including the precise outlines of the main powerhouse, main transformer chamber, and tailrace chamber, three main fracture zones (F1, F2, F3), and the regional aquifer boundary. The model's extent was set at five times the borehole diameter (approximately 1500 m × 1200 m × 800 m) extending outwards from the cavern group. Unstructured tetrahedral fine meshes were then generated for the cavern group and its surrounding rock area (extending 30 m outwards from the cave wall). The mesh size was 0.8 m at the cavern outline and a maximum size of 2 m near the surrounding rock. Structured hexahedral coarse meshes were generated for the far-field aquifer area, with unit sizes of 15 m horizontally and 8 m vertically. The nodes of the two sets of grids do not overlap at the outer boundary of the cavern group (i.e., the hydraulic contact interface), forming a typical non-matching interface.

[0056] The Mortar interface is constructed using step S3. The intersection surface between the boundary of the cavern group region and the boundary of the aquifer region is extracted. After mapping to the two-dimensional parameter domain, Delaunay triangulation is performed, and the interface mesh cell size is taken as... Due to the significant water pressure gradient at the intersection of the F1 and F2 fracture zones and the cavern, the Mortar interface in this area was locally densified to 1.5 m. After completing the sealing check and external normal unification, the Mortar interface quality matrix was stored. and adjacency table.

[0057] Construct the first projection operator according to step S4. (Cave group → Mortar interface) and second projection operator (Mortar interface → aquifer). Number of side boundary nodes of the cavern group. Number of aquifer lateral boundary nodes Mortar interface node count The projection matrix elements are calculated using a fourth-order Gauss-Legend integral. Frequent field verification shows that the reconstruction error of the projection operator for the uniform pressure field is less than [value missing]. The global flow conservation error is less than .

[0058] A seepage-stress coupled saddle point system was constructed. The constitutive model of the surrounding rock of the cavern complex adopted the Hoek-Brown equivalent Mohr-Coulomb parameters, and the aquifer rock mass was modeled using a linear elastic model. Initial permeability coefficient values ​​were derived from pressure water tests, and the time step was taken as the excavation cycle of each layer during construction, simulating a total of 12 excavation steps. Uzawa iterative solutions were used, with approximately 18-25 iterations converging within each time step. Finally, the total inflow of water into the caverns, the pore water pressure distribution around the cavern, and the volume of the plastic zone in the surrounding rock were calculated for each excavation step.

[0059] When construction reveals new fracture information or monitoring water level changes, the permeability parameters of key fracture zones on the tunnel side are updated using Bayesian theory, and the aquifer permeability coefficient field is inverted using ensemble Kalman filtering, while simultaneously optimizing the projection operator and Mortar interface. The model is then automatically updated, the coupled solution is re-executed, and the corrected inflow prediction curve and surrounding rock stability indices are output. This achieves integrated modeling throughout the entire process, from data fusion, mismatched mesh coupling, transfer of conserved physical quantities, strong coupling solution to dynamic adaptive updates.

[0060] Example 3: Please see Figure 4 , Figure 4 This is a schematic diagram of the hardware device in operation according to an embodiment of the present invention. The hardware device specifically includes: a cavern group-aquifer non-matching mesh coupling modeling device 401, a processor 402, and a storage medium 403.

[0061] A cavern group-aquifer non-matching mesh coupling modeling device 401: The cavern group-aquifer non-matching mesh coupling modeling device 401 implements the cavern group-aquifer non-matching mesh coupling modeling method.

[0062] Processor 402: The processor 402 loads and executes the instructions and data in the storage medium 403 to implement the cavern group-aquifer non-matching mesh coupling modeling method.

[0063] Storage medium 403: The storage medium 403 stores instructions and data; the storage medium 403 is used to implement the cavern group-aquifer non-matching mesh coupling modeling method.

[0064] This invention is not limited to the specific embodiments described above. Those skilled in the art can implement this invention using various other specific embodiments based on the content disclosed herein. Therefore, any design that adopts the design structure and concept of this invention and makes some simple changes or modifications falls within the scope of protection of this invention.

Claims

1. A method for modeling a cavern group-aquifer mismatched mesh coupling, characterized in that: Includes the following steps: S1. Collect multi-source detection data, perform coordinate unification, spatiotemporal registration, standardization and noise removal processing on the multi-source detection data, and construct a unified three-dimensional geological framework model including strata, cavern groups, aquifers, water-rich zones and water-conducting channels; S2. Within the unified three-dimensional geological framework model, unstructured tetrahedral fine meshes are generated for the cavern group and the area near the surrounding rock, and structured hexahedral coarse meshes are generated for the far-field aquifer area, thus constructing a non-matching dual-mesh system where nodes do not coincide and elements do not correspond at the contact interface. S3. At the hydraulic contact boundary of the mismatched dual-grid system, an independent Mortar intermediate interface is created, which serves as a virtual transition layer for the conservation and transfer of physical quantities between the mismatched grids. S4, Construct the first stage L 2 An orthogonal projection operator projects the physical quantities in the fine grid of the cavern onto the Mortar intermediate interface. Then, a second-level dual basis function projection operator is constructed to distribute the conservation information on the Mortar intermediate interface to the coarse grid of the aquifer, thereby realizing the conservation and transfer of physical quantities at the mismatched interface. S5. Using Mortar multipliers as interface constraints, a hybrid variational equation is established to form a seepage-stress coupled saddle point system. The groundwater seepage control equation, the effective stress balance equation of the surrounding rock, and the interface flux continuity equation are solved simultaneously to realize real-time bidirectional exchange and global coupling of the physical field between the cavern and the aquifer. S6. When construction reveals, supplements exploration, or monitoring data is updated, the surrounding rock parameters of the tunnel side are updated based on Bayesian theory, and the first-level L is corrected. 2 The orthogonal projection operator uses ensemble Kalman filtering to invert the aquifer parameter field and simultaneously optimizes the second-level dual basis function projection operator to achieve dynamic updating of all elements of the model and quantification of uncertainty.

2. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 1, characterized in that, Step S3, creating the Mortar intermediate interface, further includes: S31. Extract the intersection of the cavern group region boundary and the aquifer region boundary from the non-matching dual grid system, mark the cavern group side interface as the non-Mortar side, and mark the aquifer side interface as the Mortar side. S32. Map the intersection in three-dimensional space to two-dimensional parameter space, calculate the curvature of the intersection interface, and perform local encryption identification in the high curvature region. S33. Generate a mesh in the two-dimensional parameter space. The cell size is calculated as the geometric mean of the mesh sizes on both sides, and Delaunay triangulation is used to generate the surface mesh. S34. Establish the mathematical mapping relationship between the Mortar intermediate interface and the grids on both sides, define the basis functions of the cavern group side, the Mortar interface side and the aquifer side respectively, and calculate the Mortar interface quality matrix. S35. Establish the adjacency relationship between the Mortar middle interface and the meshes on both sides, adopt the region decomposition strategy to support parallel computing, and establish a version control mechanism to automatically trigger interface updates when the mesh is reconstructed. S36. Perform a closure check and external normal direction unification on the Mortar intermediate interface, verify that the interface area matches the interface area of ​​the meshes on both sides, and achieve geometric consistency and area conservation.

3. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 1, characterized in that: Step S4, constructing the dual projection operator, further includes: S41. Verify the spatial overlap between the Mortar intermediate interface and the boundary grids on both sides, establish a local coordinate system and calculate the Jacobian determinant; S42. Define the linear triangular basis functions of the cavern group side boundary element, the piecewise linear basis functions of the Mortar interface element, and the bilinear quadrilateral basis functions of the aquifer side boundary element, and configure Gauss-Legend integration points on the Mortar interface element. S43. Calculate the first projection operator from the cavern group side to the Mortar interface, construct the dual basis function that satisfies the biorthogonality condition, and then calculate the second projection operator from the Mortar interface to the aquifer side. S44. Perform constant field accuracy and global conservation checks on the first and second projection operators, check the condition number of the mass matrix, and output the projection operator that satisfies the conservation constraints.

4. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 1, characterized in that: Step S5, which involves constructing and solving the multi-field coupled system, further includes: S51. Construct the system matrices for each physical field, including the side stiffness matrix of the cavern group, the side stiffness matrix of the aquifer, the Biot coupling matrix, the side seepage system matrix of the cavern group, and the side seepage system matrix of the aquifer. S52. Based on the system matrix, a complete seepage-stress coupled saddle point system is assembled, wherein the traction force Mortar multiplier and the flow rate Mortar multiplier are introduced as Lagrange multipliers for displacement constraint and pressure constraint, respectively. S53. Based on the LBB condition, verify the stability of the Mortar interface constraint matrix. Use Uzawa iteration for global coupling calculation. Use the Mortar multiplier update equation to make the residual difference of physical quantities projected on both sides of the interface approach zero. S54. Extract the water pressure, flow rate and traction force distribution on the Mortar interface, calculate engineering indicators, and realize real-time bidirectional exchange and global coupling of the physical field between the cavern and the aquifer.

5. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 1, characterized in that, The multi-source detection data collected in step S1 includes ground-penetrating radar data, borehole test data, hydrological test data, and groundwater monitoring data.

6. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 3, characterized in that, In step S6, the dynamic update of all elements of the model includes: updating the surrounding rock parameters of the cavern side and correcting the first projection operator, inverting the aquifer permeability coefficient field and pore water pressure field and simultaneously optimizing the second projection operator, and finally realizing the synchronous update of the grid, Mortar interface, projection operator and coupling field, and outputting the parameter confidence interval, field distribution probability and risk level.

7. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 4, characterized in that: The convergence condition for the Uzawa iteration in step S53 is as follows: the displacement difference and pressure difference projected onto both sides of the interface are used as residuals, and convergence is determined when the residuals approach zero; the Mortar multiplier update equation is: , in and These are the Mortar multipliers for traction force and the Mortar multiplier for flow rate, respectively. The Mortar interface quality matrix. and These are the displacement residual and pressure residual after the m-th iteration, respectively. and This is the iterative relaxation factor.

8. The method for modeling a cavern group-aquifer mismatched mesh coupling according to claim 1, characterized in that: In step S2, the unit size of the unstructured tetrahedral fine mesh at the cavern outline is controlled within a first preset range, and the unit size of the structured hexahedral coarse mesh is controlled within a second preset range. The two types of meshes form a natural mismatched interface at the hydraulic contact interface.

9. A storage medium, characterized in that: The storage medium stores instructions and data to implement the cavern group-aquifer non-matching mesh coupling modeling method according to any one of claims 1 to 8.

10. A device for modeling mismatched meshes between cavern groups and aquifers, characterized in that: include: A processor and a storage medium; the processor loads and executes instructions and data in the storage medium to implement the cavern group-aquifer non-matching mesh coupling modeling method according to any one of claims 1 to 8.

Images