A calculation method for non-coordinated deformation of cushion material and face plate of high face rockfill dam in alpine region

Abstract

The application discloses a kind of high-cold area high faceplate rockfill dam cushion material-faceplate inharmonious deformation calculation method.The method obtains the structure parameters and operating boundary conditions of high faceplate rockfill dam, and establishes multi-physical field domain model with non-matching grid interface.A unified projection base operator based on contact mechanics is constructed at the non-matching grid interface, which is used to map the physical parameters across the interface to obtain the generalized jump and dynamic contact gap of the interface.According to the dynamic contact gap, the evolution state of the interface is determined, and then the equivalent thermal impedance and hydraulic conductivity of the interface are obtained through the gate feedback calculation to construct the bidirectional feedback interface term.Finally, the augmented control equation is constructed by combining the preset contact constraint condition, and the nonlinear simultaneous iteration is solved to output the multi-physical field distribution and inharmonious deformation state of the dam.The application realizes high-fidelity simulation of water and heat blockage and dredging effect caused by panel disengagement under strong coupling conditions, and improves the safety evaluation accuracy of cold region dams.

Description

Technical Field

[0001] This invention relates to the field of multiphysics numerical simulation technology for water conservancy projects, and in particular to a method for calculating the incoordination deformation of the cushion material and the face panel of a high-face rockfill dam in a cold region. Background Technology

[0002] Due to their advantages such as strong adaptability to terrain, fast construction speed, and relatively low project cost, rockfill dams with concrete face panels have been widely used in mountainous water conservancy and hydropower projects. Especially in high-altitude and cold regions, the dam body is subjected to low temperatures and frost heave for a long time. The temperature field, moisture field, and displacement field between the foundation material and the concrete face panel interact with each other, which can easily lead to problems such as frost heave deformation, contact cracking, and seepage evolution, thereby affecting the safety and stability of the dam body.

[0003] The face panel and the subbase material of a rockfill dam have a complex contact relationship. Due to the significant difference in material stiffness between the face panel and the subbase material, incompatible deformation problems such as separation, local voids, and tangential slippage may occur at the contact interface. Existing engineering analyses often use simplified methods such as node binding, penalty function contact, or rigid connections, which can easily lead to problems such as local penetration and discontinuous contact stress distribution under mismatched mesh conditions, making it difficult to accurately describe the actual contact force transmission process between the face panel and the subbase material. Especially under high-altitude and cold conditions, the contact state evolves continuously with temperature changes, moisture migration, and frost heave deformation, and traditional methods often struggle to balance computational accuracy and solution stability.

[0004] Therefore, there is an urgent need to propose a water-thermal-mechanical coupling analysis method suitable for high-face rockfill dams in cold regions. This method should rationally describe the incoordination deformation and contact response caused by the difference in material stiffness between the face material and the cushion material, taking into account freeze-thaw phase changes, unfrozen water migration, and frost heave deformation. This would improve the accuracy and stability of the analysis of the dam's temperature, moisture, and displacement fields, providing technical support for the design analysis, operational safety evaluation, and related disease research of high-face rockfill dams in cold regions. Summary of the Invention

[0005] Purpose of the invention: To provide a calculation method for the mismatch deformation between the subbase material and the panel of a high-altitude rockfill dam in cold regions, in order to solve the above-mentioned problems in the prior art.

[0006] Technical solution: A calculation method for the mismatch deformation of the subbase material and the face panel of a high-altitude rockfill dam in cold regions, including: Obtain the structural parameters and operational boundary conditions of a high-face rockfill dam; A multiphysics domain model with a non-matching mesh interface is established based on structural parameters. The multiphysics domain model includes a moisture field equation with a phase change source, a temperature field equation with latent heat, and a displacement field equation based on effective stress. Construct a unified projection basis operator based on contact mechanics at the interface of the mismatched mesh; By using the unified projection basis operator to perform cross-interface mapping of the physical parameters of displacement field, temperature field and moisture field, the generalized jump quantity and dynamic contact gap of the interface are obtained. The evolution state of the interface is determined based on the dynamic contact gap; Gated feedback calculations are performed based on the evolutionary state and generalized jump amount to obtain the interface equivalent thermal impedance and interface hydraulic conductivity, so as to construct a bidirectional feedback interface term. Augmented control equations are constructed based on a multiphysics domain model, bidirectional feedback interface terms, and preset contact constraints. The augmented control equations are solved nonlinearly and iteratively to output the multiphysics distribution and incoordination deformation state of the high panel rockfill dam.

[0007] Beneficial effects: By introducing a unified projection basis operator based on Mortar contact mechanics, high-precision cross-interface mapping of physical parameters of displacement, temperature, and moisture under non-matched mesh conditions is achieved; by constructing a gated feedback mechanism, the interface thermal impedance and hydraulic conductivity dynamically evolve with the contact state, overcoming the shortcomings of the traditional fixed boundary assumption that cannot reflect the sudden changes in physical flux caused by delamination; by constructing a globally consistent solution framework using nonlinear complementary functions and the semi-smooth Newton method, numerical oscillations caused by contact state switching are avoided, and the solution stability of strongly coupled nonlinear problems is improved. Attached Figure Description

[0008] Figure 1 This is a flowchart of the calculation process of the present invention.

[0009] Figure 2 This is a geometric diagram of a panel rockfill dam.

[0010] Figure 3 This is a grid diagram of a panel rockfill dam.

[0011] Figure 4 This is a water distribution cloud map of a panel rockfill dam.

[0012] Figure 5 This is a temperature field cloud map of a panel rockfill dam.

[0013] Figure 6 This is a displacement field cloud map of a panel rockfill dam.

[0014] Figure 7 It is a displacement curve graph.

[0015] Figure 2 In the attached diagram, the labels are: toe slab 1, concrete surface layer 2, cushion layer 3, transition layer 4, main rockfill layer 5, secondary rockfill layer 6, dam crest 7. TP1 to TP8 are measuring points 1 to 8. Figure 4 In the middle, 8 represents the wetting line. Detailed Implementation

[0016] To address the problems existing in the prior art, the main technical concept, main formulas, and model of this application are described as follows: A water-thermal-mechanical coupling calculation model for high-face rockfill dams in cold regions was established, considering the effects of frost heave of the subbase material and reservoir water level. Contact mechanics calculation methods were introduced into the water-thermal-mechanical coupling calculation, incorporating the contact interaction between the face panel and the subbase material into the water-thermal-mechanical coupling calculation framework. This model can more realistically reflect the redistribution of contact pressure between the face panel and the subbase material, separation of the face panel and the subbase material, voids, and uncoordinated deformation problems such as frost heave failure of the face panel under the influence of frost heave of the subbase material and reservoir water level.

[0017] The frost heave effect of the subbase material is calculated by obtaining the volumetric ice content at each unit through water-thermal coupling calculation; the volumetric ice content is used as a key state variable and substituted into the thermal expansion effect to calculate the frost heave strain; the frost heave strain is introduced as an additional strain term into the equilibrium equation of the displacement field and participates in the nonlinear iterative solution.

[0018] The reservoir water level effect is determined based on the actual operating water level of the high-face rockfill dam in the high-altitude cold region. The head conditions at the upstream and downstream boundaries are determined and used as variable inputs. The matrix potential-volume water content relationship function and the volumetric water content-permeability coefficient relationship function are established. The above relationship functions are applied as boundary head conditions to the saturated-unsaturated seepage control equation for solution. The volumetric water content and permeability coefficient distribution at each location of the dam body are obtained through iterative calculation.

[0019] Specifically, the head boundary conditions are determined based on the reservoir water level time series in the operational boundary conditions, and the equivalent relative permeability coefficient and the head boundary conditions are substituted into the water field equation containing the phase change source to perform saturated-unsaturated seepage calculations.

[0020] The periodic changes in reservoir water level are described using the Van Genuchten model. This model dynamically updates the volumetric water content, relative permeability coefficient, and hydraulic parameters based on the pressure head, overcoming the limitation of fixed seepage parameters in reflecting the changes in the unsaturated seepage state of the dam body under conditions of reservoir water level fluctuations and frost heave of the foundation material. This improves the calculation accuracy of moisture migration in the unsaturated zone of the dam body and its water-thermal-mechanical coupling process. Furthermore, this method can reasonably reflect the repeated freeze-thaw effects of periodic reservoir water level changes on the foundation material of high-face rockfill dams in cold regions under conditions of panel defects, and accurately describe the process of frost heave deformation of the foundation material continuously intensifying and the accelerated development of uncoordinated deformation between the panel and the foundation material under repeated freeze-thaw cycles.

[0021] Contact mechanics calculations are performed using either the multiplier transformation method or the augmented Lagrange method. Taking the multiplier transformation method as an example, its basic idea is as follows: Based on the principle of minimum potential energy, a variational equation is established; contact constraints are applied, introducing constraints for the normal non-penetration condition, the tangential bonding point set, and the sliding friction point set; Lagrange multipliers are condensed through local basis transformations, defining generalized displacements for each node in the active contact state; further, a global transformation matrix P is constructed to transform the multiplier degrees of freedom introduced by the contact constraints into equivalent constraint terms, thereby obtaining linearized governing equations suitable for iterative solutions using the Newton-Raphson method.

[0022] The reservoir water level affects the pore water pressure and pore ice pressure of the subbase material, while the frost heave of the subbase material affects the deformation of the framework, thus influencing the distribution and magnitude of the total stress. These two factors jointly determine the effective stress state of the subbase material through the principle of effective stress. Effective stress is transferred through the contact between the subbase material and the panel, forming the effective contact stress at the interface. When the effective contact stress at the interface is high, the panel and subbase material maintain continuous contact and the load is transferred normally. When the effective contact stress at the interface decreases to near zero, the panel separates from the subbase material, resulting in a void. The panel loses local support, and the surrounding stress redistributes, causing the void to expand in adjacent areas. The voided panel, under its own weight, water pressure, and the frost heave of the subbase material, experiences excessive bending tensile stress, easily inducing cracking, arching, and joint misalignment.

[0023] Furthermore, the contact constraint conditions are as follows: 1) Impenetrable: ; 2) Tangential bonding: ; 3) Friction and sliding: ; In the formula: n is a point on the surface The unit outward normal vector at that location, For points on the surface, To indicate from point The projection along the n-direction onto the principal plane. For elements that are not in the Mortar matrix D, For the elements of the Mortar matrix M, A This is the set of nodes that actually come into contact in the current iteration step; To be from the surface point The unit tangent vector at that point, Г This is the set of nodes that have not yet been visited. The coefficient of friction is given by the Lagrange multiplier, which represents the contact force at a point. For tangential contact force, For normal contact force, For generalized displacement increment, Let be the set of friction sliding points.

[0024] The local basis transformation is applied to each active sub-node. j Generalized displacement is defined : (1); In the formula: Let the generalized displacement be from node j on the surface. Let be the original nodal displacement vector from node j on the face. The elements are the inverse of the non-Mortar matrix D. For the elements of the Mortar matrix M, j , l It is a natural number.

[0025] Define the global transformation matrix P: (2); In the formula: I represents the identity matrix. N It is a set of non-contact nodes. M The set of main face nodes, Г The set of nodes that have not yet been visited. A This is the set of nodes that actually come into contact in the current iteration step; For the reason A matrix composed of elements.

[0026] The linearized governing equations obtained using the Newton-Raphson method are: (3); In the formula: , Finite element stiffness matrix For generalized displacement based on the local coordinate system; , This represents the displacement increment between different iteration steps; , The residual vector is represented by the equivalent nodal forces. It is a local Lagrange multiplier.

[0027] The method overcomes the problems of increased degrees of freedom, compromised positive definiteness of coefficient matrix, and reduced solution efficiency caused by Lagrange multipliers being directly involved in the solution of global equations as independent unknowns. It achieves point-by-point decoupling of contact constraints and improves the stability and efficiency of contact calculation.

[0028] Furthermore, the water-thermal-mechanical coupling calculation method: Based on Fourier's law and the law of conservation of energy, the temperature field governing equations are constructed, and water migration and ice-water phase transition are introduced into the heat transfer term; the water field governing equations are constructed based on the two-dimensional Richards saturated-unsaturated equations; to achieve coupling between the temperature field and the water field, an empirical relationship for the unfrozen water content is introduced, and the solid-liquid ratio is used. B I The ratio of pore ice to unfrozen water is described; the displacement field is described by the equilibrium equation, geometric equation and constitutive equation respectively.

[0029] The temperature field control equation is expressed as: (4); In the formula: The volumetric heat capacity of the soil, kJ / (m³) 3 ·℃); T is the transient temperature of the soil, ℃; t is time, s; is the thermal conductivity of the soil, W / (m·℃); V The seepage velocity is in m / s; For the Hamiltonian operator, for the two-dimensional problem, take [ / x, / y).

[0030] The governing equation for the moisture field is expressed as: (5) In the formula: H For the total water flow, S f For water storage rate, V The velocity of water migration in the frozen zone of the dam body is expressed in m / s.

[0031] The hydrothermal coupling equation is expressed as: (6) In the formula: This represents the initial volumetric water content of the soil. negative temperature T At that time, the volumetric water content of unfrozen water; T f The temperature at which water turns into ice during soil freezing, in °C; B It is a constant and depends on the soil type and salt content.

[0032] At the same time, the solid-liquid ratio is introduced. B I , expressed as ice content by volume θ I With volumetric moisture content θ u The ratio.

[0033] (7); In the formula: θ I Ice content by volume θ u Volumetric water content The control equations for the displacement field are expressed as follows: (8); In the formula: For the total stress, External force; In response, This represents the total displacement. For stress increment, The temperature-dependent elastic matrix of frozen soil; These represent the strain increments.

[0034] In summary, the main idea is as follows: obtain the geometric parameters and operational boundary conditions of the high-face rockfill dam, establish a multiphysics domain model with a mismatched mesh interface between the face area and the cushion material area, and the multiphysics domain model includes a moisture field containing a phase change source, a temperature field containing latent heat, and a displacement field based on effective stress. At the non-matching mesh interface, a unified projection basis operator based on the Mortar contact mechanics mechanism is constructed to map the physical quantities of displacement field, temperature field and moisture field through the unified projection basis operator to obtain the generalized jump amount across the interface and the dynamic contact gap of the interface. Based on the dynamic contact gap, the contact or de-void evolution state of the interface is determined; and based on the gating mechanism, the contact or de-void evolution state is fed back to the calculation of the interface equivalent thermal impedance and the interface leakage hydraulic conduction capacity, forming a two-way feedback between interface contact response and heat flow and moisture migration. An augmented control equation consistent with the multiphysics field variables and the contact state set is constructed. The volume control model, the bidirectional feedback interface term and the contact constraint conditions are solved by a fully coupled nonlinear simultaneous iterative solution. The multiphysics field distribution law of the dam body and the uncoordinated deformation state between the panel and the cushion material are output.

[0035] The unified projection basis operator refers to the integral projection matrix constructed based on the biorthogonal dual function in Mortar contact mechanics. Its function is to establish a high-precision mapping channel that satisfies energy conservation between the two nodes of the non-matched mesh interface, so that the physical parameters of the mechanical field, temperature field and moisture field can be transferred across the interface through the same set of matrix operators.

[0036] The pre-defined contact constraints in the augmented governing equations are implemented using the multiplier transformation method, specifically including: Define the non-penetration condition in the interface normal and the frictional sliding condition in the tangential direction; Based on local basis transformation, a generalized displacement is defined for interface nodes in the active contact state; Construct a global transformation matrix based on generalized displacement; By using the global transformation matrix, the Lagrange multiplier degrees of freedom introduced by solving the contact constraints are transformed into equivalent constraint terms. The Lagrange multipliers are eliminated, and the equivalent control equations after eliminating the multiplier degrees of freedom are obtained. This completes the point-by-point decoupling of the augmented control equations and satisfies the non-penetration condition and the friction sliding condition.

[0037] The pre-defined contact constraints in the augmented governing equations are implemented using the augmented Lagrangian method, specifically including: Define the master face and slave face on both sides of the non-matching mesh interface to construct interface contact pairs; Set the initial values ​​of the pre-configured penalty factor, Lagrange multiplier, and allowable penetration; In the nonlinear iterative solution of the augmented control equation, based on the interface contact pairs, the interface is made to approximate the non-penetration condition and friction constraint condition by alternating update calculation of the penalty factor and Lagrange multiplier, thereby obtaining the interface mechanical response that satisfies the contact constraint and realizing the control of the contact constraint.

[0038] The main details of this framework are described in detail below through examples.

[0039] like Figure 1 As shown in Example 1: This example provides a calculation method for considering the uncoordinated deformation of high-face rockfill dams in cold regions. The method includes the following steps: Step 101: Obtain the structural parameters and operational boundary conditions of the high-face rockfill dam; Specifically, structural parameters are physical parameters that characterize the geometric topology and material partitioning of the dam entity.

[0040] In the preliminary modeling stage of actual engineering, the structural parameters obtained may include the maximum dam height, dam crest length, and the spatial coordinate range of the dam body's panel area, cushion layer area, transition zone, and main rockfill area. Running boundary conditions are used to constrain the initial state of the multiphysics field's temporal evolution and external driving sources.

[0041] The acquired operational boundary conditions can specifically include time series of reservoir water level elevations that change dynamically over time, and time series of air temperatures monitored by meteorological stations. Obtaining these parameters provides a fundamental data source for subsequently constructing a numerical model that conforms to the scale of the actual project.

[0042] Step 102: Establish a multiphysics domain model with a non-matching mesh interface based on structural parameters and operational boundary conditions. The multiphysics domain model includes a moisture field equation containing a phase change source, a temperature field equation containing latent heat, and a displacement field equation based on effective stress. Specifically, based on the aforementioned structural characteristics, mesh discretization is performed in the spatial domain. Since the concrete panel material is dense and has a large deformation gradient, while the subbase material exhibits granular properties, a non-matching locally refined mesh is used for spatial discretization at the interface between the two materials to adapt to the solution accuracy requirements of different materials.

[0043] For the multiphysics domain model, a theoretical framework is constructed using a set of partial differential equations conventional in this field. Specifically, the water field equation containing a phase change source is realized using the governing equations that include the phase change source term and the mass conservation constitutive model; the temperature field equation containing latent heat is realized using the energy conservation equation based on the equivalent apparent heat capacity method; and the displacement field equation based on effective stress is realized using the modified effective stress formula that includes pore water pressure and pore ice pressure, as well as the momentum balance equation.

[0044] The differential assembly and solution rules for the relevant governing equations adopt the existing general standards of computational mechanics, and will not be described in detail here.

[0045] Step 103: Construct a unified projection basis operator based on contact mechanics at the interface of the mismatched mesh; In multiphysics computations, direct data transfer can lead to non-conservation of physical flux due to non-coincident grid nodes on both sides. To address this, a contact mechanics dual operator is extracted across the non-matching grid interface and used as a unified projection basis operator.

[0046] This operator is used to construct a high-precision integral mapping channel between master and slave nodes in mathematical space, ensuring the conservation and transfer of energy and mass at mesh discontinuities.

[0047] Step 104: Use the unified projection basis operator to perform cross-interface mapping on the physical parameters of displacement field, temperature field and moisture field in the multiphysics domain model to obtain the generalized jump quantity and the dynamic contact gap of the interface. In the specific implementation, matrix multiplication is performed by a unified projection basis operator to map the node parameters of the panel side and the subbase material side to the same integral surface for subtraction.

[0048] This allows for the calculation of a generalized jump quantity, including displacement jump, interface temperature jump, and water potential jump.

[0049] Furthermore, the component of the displacement jump in the interface normal direction is extracted and used as the dynamic contact gap to quantitatively characterize the physical distance of separation between the panel and the padding layer.

[0050] Step 105: Determine the evolution state of the interface based on the dynamic contact gap; Specifically, a judgment logic is established that compares the dynamic contact gap with a preset zero-point tolerance. When the dynamic contact gap is equal to or less than the zero tolerance, the evolution state is determined to be a contact state; when the dynamic contact gap is greater than the zero tolerance, the evolution state is determined to be a delamination state. The result of this state determination directly determines the mechanical support characteristics and medium conductivity properties at the interface.

[0051] Step 106: Perform gated feedback calculations based on the evolution state and generalized jump amount to obtain the interface equivalent thermal impedance and interface hydraulic conductivity in order to construct a bidirectional feedback interface term. In this embodiment, the interface is no longer considered a static, constant boundary, but is designed as a dynamic physical gating mechanism driven by contact states. When a contact state is determined, the heat and moisture at the interface are directly conducted through the solid skeleton and porous medium. When a void state is determined, the physical gating triggers nonlinear impedance adjustment, using dynamically calculated nonlinear impedance to cut off the heat conduction path and amplifying the hydraulic conduction channel based on the void gap width. The equivalent thermal impedance and hydraulic conduction capacity of the interface obtained through the above process are assembled at the equation level into a bidirectional feedback interface term with bidirectional coupling characteristics.

[0052] Step 107: Construct augmented control equations based on the multiphysics domain model, bidirectional feedback interface terms, and preset contact constraints. Specifically, the bulk stiffness matrix characterizing the macroscopic continuous medium, the diagonal feedback flux matrix characterizing the interface gating effect, and the preset contact constraints used to limit the mutual penetration of mesh nodes are assembled. This integrates the independent moisture field, temperature field, displacement field, and discontinuous interface states to generate globally aligned augmented governing equations.

[0053] Step 108: Solve the augmented control equations nonlinearly using simultaneous iterations to output the multiphysics distribution and incoordination deformation state of the high panel rockfill dam.

[0054] Furthermore, a numerical algorithm employing global tangent differentiation and iterative approximation is used to solve the augmented governing equations. After the residuals meet the convergence requirements, the temperature gradient field, ice content distribution field, and stress-strain distribution field within the entire dam body space are obtained. Simultaneously, the incoordination deformation state at the interface is output, realizing a closed-loop deduction from meteorological and hydrological driving forces to the internal physical response of the dam.

[0055] Compared to calculations using traditional node-binding or fixed penalty function contact simplification, this method more accurately reflects the temperature field isolation effect and concentrated moisture leakage paths in the voided area of ​​the panel. Specifically, traditional methods, assuming complete interface continuity, cannot capture the abrupt changes in interfacial heat flux density and the redistribution of the seepage field after voiding occurs. In contrast, this method achieves the dynamic evolution of physical flux with contact state through gated feedback interface terms, thus achieving higher computational fidelity in interfacial temperature gradient prediction and voiding range identification.

[0056] In one possible implementation, after outputting the multiphysics distribution and incompatible deformation state of the high-panel rockfill dam, the following steps are also included: Step 109: Based on the output of the non-coordinated deformation state, extract the spatial distribution range of the panel void area and the interfacial frost heave pressure of the padding material on the panel. Specifically, the incompatible deformation state contains rich information on both microscopic and macroscopic defects. By traversing the node state set of the mismatched mesh interface, the geometric coordinates of nodes currently in a voided state are selected and compiled to generate the spatial distribution range of the voided area of ​​the panel. At the same time, the normal compressive stress tensor of the interface is extracted, and the additional compressive stress component caused by ice crystal expansion is separated and extracted as the interface frost heave pressure.

[0057] Step 110: Combining the pre-stored structural strength threshold and interface frost heave pressure of the high-face rockfill dam, a cracking risk assessment is conducted on the concrete panel to obtain the cracking risk assessment results.

[0058] In this step, the safety of the dam structure is assessed using the extracted mechanical loads. The pre-stored structural strength threshold is specifically the allowable tensile and compressive strength limits of the concrete material. Assuming that the normalized value of the currently extracted interfacial frost heave pressure is 0.85 after normalization, and the pre-configured structural strength threshold is set to 0.70, since 0.85 > 0.70, a high risk of local stress overload is determined at the current dam node, thus triggering a panel cracking risk warning for that area.

[0059] Step 111: Based on the cracking risk assessment results and the spatial distribution range of the voided area, generate dam operation safety evaluation indicators and issue early warnings.

[0060] In the specific implementation, the high-risk coordinate point set of the cracking risk assessment results is spatially overlaid and mapped with the spatial distribution range of the voided area. If a high-risk cracking point is located exactly at the voided boundary, it indicates a tendency for shear failure caused by voiding. Based on this overlay result, a structured dam operation safety evaluation index is generated by calculating the risk area ratio or stress concentration coefficient. Finally, this evaluation index is converted into a visualized alarm layer or early warning log stream output, providing direct data support and decision-making basis for the safe operation and scheduling of water conservancy projects during winter.

[0061] According to one aspect of this application, it further includes: spatially superimposing and mapping the spatial distribution range of the cracking risk area and the void area, calculating the void area ratio and / or stress concentration coefficient as dam operation safety evaluation indicators, and issuing an early warning when the safety evaluation indicators exceed the preset safety warning threshold.

[0062] This embodiment addresses several issues in existing technologies, including: Existing simulation methods often struggle to accurately reflect the transient transfer mechanisms of interface physical properties when dealing with large, discontinuous deformations caused by extreme conditions. In the critical region of multiphysics interactions, computational models often fail to balance the discontinuous evolution of macroscopic meshes with the conservation of microscopic flux, leading to predicted boundary physical responses deviating from actual engineering monitoring data. This results in unreliable temporal evolution support for refined structural safety early warning under harsh environments.

[0063] Example 2: Further explanation of the specific physical response mechanism and modified model in the water field equation containing the phase change source and the displacement field equation based on effective stress.

[0064] In the computational framework of this embodiment, the underlying physical computing environment is uniformly defined using the International System of Units (SI) for dimensional conversion and parameter matching.

[0065] For the existing fundamental governing equations in the multiphysics domain model, the conventional partial differential equations in this field are adopted: the temperature field is solved using the energy conservation equation based on the equivalent apparent heat capacity method; the moisture field is solved using the Harlan-type mass conservation equation that includes the ice-water phase transition source term and the Darcy constitutive model; and the displacement field is solved using the Bishop modified effective stress formula that includes pore water pressure and pore ice pressure, as well as the Clausius-Clapeyron equation that describes the thermodynamic equilibrium of the ice-water phase transition.

[0066] Based on the aforementioned conventional framework, this embodiment focuses on achieving gating and blocking of physical mechanisms and parameter correction through the following computational steps. In one possible implementation, the steps include: Step 201: Based on the calculation results of the temperature field in the multiphysics domain model, obtain the volumetric ice content and unfrozen water content of the soil, and obtain the pre-configured porosity. Specifically, the volumetric ice content and unfrozen water content are direct functions of the temperature field. During the multiphysics iteration process, the current temperature value at each grid integration point is extracted, and the aforementioned moisture characteristic variables are solved in conjunction with material property parameters. In particular, for the calculation of unfrozen water content, a segmented judgment rule with the freezing temperature as the boundary is introduced. When the current temperature is determined to be greater than or equal to a pre-configured freezing temperature threshold, the unfrozen water content is directly equal to the initial volumetric water content; when the current temperature is determined to be strictly less than the freezing temperature threshold, based on the power function decay law, the ratio of the current temperature to the freezing temperature threshold is calculated, and this ratio is multiplied by the initial parameters after taking its negative exponent to obtain the unfrozen water content. This segmented processing avoids the numerical divergence problem that may occur when solving empirical formulas under positive temperature conditions.

[0067] Step 202: Based on the volumetric ice content, construct an exponential decay function to characterize the ice blockage effect, as the ice resistance factor; In this step, the ice resistance factor is used to characterize the nonlinear resistance effect of the gradual truncation of moisture transport channels as the volume of ice crystals in the pores increases. Specifically, the volumetric ice content is used as the volume fraction of the pore space occupied by ice, and an exponential decay model is constructed based on empirical material parameters. The specific formula is as follows: Ω i =10 (-Ω*θI) ; Among them, Ω i θ is the constructed ice impedance factor, Ω is the pre-configured empirical impedance parameter, and θ is the... I To obtain the volumetric ice content. Specifically, for coarse-grained riprap, the empirical impedance parameter Ω can be set in the range of 3 to 5; for fine-grained subbase material, this parameter can be set in the range of 5 to 8, in order to reflect the physical objective law that fine-grained soil is more easily blocked by ice crystals.

[0068] Step 203: Based on the unfrozen water content, the matrix potential-volume water content relationship function and the volumetric water content-permeability coefficient relationship function established by the VanGenuchten model are used to calculate the basic relative permeability coefficient, and the ice resistance factor is used to correct the basic relative permeability coefficient to obtain the equivalent relative permeability coefficient. In the specific implementation process, based on the Van Genuchten theoretical model, the basic relative permeability coefficient under ice-free conditions is calculated using the unfrozen water content. Subsequently, a multiplication operation is performed, multiplying the basic relative permeability coefficient by the ice resistance factor generated in step 202, thereby downgrading and correcting the original hydraulic conduction parameters.

[0069] Under one operating condition, the calculation formula is as follows: K eq =K base *Ω i ; Among them, Keq K is the equivalent relative permeability coefficient. base Based on the relative permeability coefficient, Ω i This represents the ice resistance factor. Through this multiplication operation, the physical cutoff effect of ice development on the seepage field is transformed into a coefficient attenuation matrix in numerical calculations.

[0070] Step 204: Substitute the equivalent relative permeability coefficient into the water field equation containing the phase change source to perform seepage calculation and obtain the pore water pressure distribution of the dam body.

[0071] By updating the equivalent relative permeability coefficient to the hydraulic conduction block matrix of the global stiffness matrix, the water field solver is driven to perform spatial seepage calculations for the next time step, thereby obtaining the pore water pressure distribution and water flux field distribution at the current moment.

[0072] Step 205: Based on the current iterative stress results of the displacement field equation, obtain the normal effective constraint compressive stress and the pre-configured reference constraint compressive stress at the interface of the cushion material. When analyzing the frost heave effect of the foundation material, the stress tensors transmitted to the grid nodes by the dam's self-weight and external water pressure are extracted and projected along the normal direction of the contact interface to obtain the effective normal constraint compressive stress. The reference constraint compressive stress represents the baseline stress level of the material's internal resistance to the expansion potential energy of microscopic ice crystals and belongs to the intrinsic strength property of the material.

[0073] Step 206: Based on the relative relationship between the effective normal constraint compressive stress and the reference constraint compressive stress, construct an exponential suppression function to characterize the gravitational constraint suppression effect; In practice, the amount of frost heave deformation is not solely determined by temperature. When there is significant compressive stress above the soil, the volume expansion of ice crystals is physically constrained. Therefore, a function characterizing this attenuation feature is constructed by calculating the ratio of the effective constrained compressive stress to the reference constrained compressive stress, taking a negative sign, and performing exponential mapping calculations.

[0074] In one example, the exponential inhibition function f n =exp(-σ n / σ ref ); Where, σ n For the effective normal confined compressive stress, σ ref The reference constraint compressive stress can be determined based on the initial swelling stress measured in indoor frost heave tests of the subbase material, or based on the recommended value of the initial constraint pressure for frost heave of the subbase material in relevant water conservancy industry standards. This function shows a monotonically decreasing trend with increasing normal compressive stress, thus quantifying the physical inhibition law of the overlying load on frost heave deformation.

[0075] Step 207: Based on the calculation results of the temperature field in the multiphysics domain model, obtain the volumetric ice content of the soil, and combine it with the pre-configured frost heave coefficient to convert the volumetric ice content into free frost heave strain. In this step, taking into account the physical fact that water increases in volume by 9% when it freezes, the obtained volume ice content is multiplied by the frost heave coefficient, which characterizes the volume expansion rate, to calculate the free frost heave strain without external force constraints.

[0076] Furthermore, due to the influence of the panel structure, the frost heave at the contact interface exhibits obvious anisotropy.

[0077] Therefore, an anisotropic distribution mechanism is introduced when distributing the three-dimensional strain tensor. For example, the frost heave strain distribution coefficient in the interface normal is set to 0.6, while the tangential distribution coefficient is set to 0.2, to reflect the physical tendency of water to preferentially accumulate and freeze along the direction of maximum temperature gradient (i.e., the normal).

[0078] Step 208: The free frost heave strain is attenuated and corrected using an exponential suppression function to obtain the confined frost heave strain, and the confined frost heave strain is introduced into the displacement field equation based on the effective stress.

[0079] Specifically, the free frost heave strain is multiplied by an exponential suppression function to calculate the partially suppressed restricted frost heave strain.

[0080] Furthermore, in the constrained computational logic, the volumetric expansion deformation suppressed by the function does not disappear. Instead, it is transformed into internal pore ice pressure through the volumetric modulus and introduced as an additional load vector into the displacement field equation based on effective stress for mechanical equilibrium calculation.

[0081] In an alternative implementation, the calculation of the frost heave effect on the subbase material further includes the following steps describing the cumulative effect of long-term freeze-thaw damage: Step 209: During the time-progression solution process, monitor the historical process of the soil temperature output by the multiphysics domain model crossing the preset freezing temperature threshold, and record the number of freeze-thaw cycles. In the specific implementation, at the end of each time evolution step of the solver, the current temperature is compared with the temperature of the previous time step. When it is detected that the temperature value changes from above the preset freezing temperature threshold to below the threshold, and then rises back to above the threshold in a complete trough cycle, the internal counter is triggered to perform an auto-increment operation, and the generated integer value is stored and recorded as the number of freeze-thaw cycles.

[0082] Step 210: Construct a logarithmic cumulative factor based on the number of freeze-thaw cycles; After repeated freeze-thaw cycles, the residual deformation caused by the loosening of the soil structure in nature exhibits a nonlinear convergence trend. Therefore, the natural logarithm of the result is calculated by adding the number of freeze-thaw cycles to a constant 1, generating a logarithmic cumulative factor.

[0083] Step 211: Calculate irreversible residual frost heave strain using a logarithmic cumulative factor and a pre-configured cumulative frost heave coefficient; In the specific implementation process, the residual frost heave strain ε r =ε1*α c *ln(1+Nf); Where ε1 is the initial frost heave strain recorded when the soil undergoes its first freeze-thaw cycle, and α c Nf represents the pre-configured cumulative frost heave coefficient obtained from the indoor cyclic test fitting, and Nf represents the number of recorded freeze-thaw cycles.

[0084] Step 212: The residual frost heave strain and the confined frost heave strain are superimposed to calculate the total frost heave strain, which is then used to replace the confined frost heave strain in the calculation of the displacement field equation based on the effective stress.

[0085] In this step, the residual frost heave strain vector with irreversible damage characteristics is vector-added with the constrained frost heave strain vector obtained from the real-time calculation of the current temperature field. After generating the total frost heave strain, it is used as the initial thermal strain condition for driving the displacement field, updating the global mechanical residual vector, and thus promoting the calculation of macroscopic non-coordinated deformation of the dam's solid skeleton.

[0086] Example 3: Further explanation of the specific calculation process of interface hydraulic conductivity and interface equivalent thermal impedance in the gated bidirectional feedback mechanism.

[0087] In the de-void state, the physical behavior at the interface of the dam structure will transform from a continuous medium model to a discontinuous gap channel model. One possible implementation includes the following steps: Step 301: Extract the water potential jump from the generalized jump; Since the generalized jump quantity output by the preceding steps exists in the form of a column vector matrix formed by packaging the degrees of freedom of multiple physics fields, the component representing the difference in potential energy of the water field is isolated separately by performing a preset row index matrix slicing operation in the program memory, and used as the water potential jump quantity for subsequent boundary flux calculation.

[0088] Step 302: Determine the conditions for the evolution state and trigger the corresponding hydraulic calculation path.

[0089] In the first operating condition, when the evolution state is a contact state, the interfacial hydraulic conductivity is calculated based on the pre-configured interfacial leakage coefficient and water potential jump. Specifically, when the subbase material and the panel are in a tightly compacted and closed geometric distribution, there are no macroscopic physical gaps at the grid interface. At this time, the cross-interfacial transfer of water depends on the micropores of the matrix soil on both sides. Using the interfacial leakage coefficient provided by the constant tensor matrix, a matrix multiplication is performed directly with the currently extracted water potential jump to obtain the interfacial water exchange flux under normal conditions.

[0090] In the second operating condition, when the evolutionary state is a voided state, based on the cubic law of fluid dynamics, the dynamic contact gap is converted into an equivalent hydraulic conductivity coefficient to describe the permeability within the voided fissure surface. Based on the equivalent hydraulic conductivity coefficient and the dynamic contact gap, an additional water flow recharge term is calculated. The interface leakage coefficient and the additional water flow recharge term are then superimposed, and combined with the water potential jump, to calculate the interface hydraulic conductivity under voided conditions. Specifically, voiding leads to the creation of a macroscopic physical space with significant thickness on the back of the panel. The flow of reservoir water along this space no longer follows the seepage law of porous media, but transforms into a fissure channel flow mechanism.

[0091] Based on the principles of fluid mechanics, the equivalent hydraulic conductivity coefficient K gap =(ρ w *g*g n 2 ) / (12*μ w ); Where, ρ w Let g be the physical density constant of liquid water, and g be the gravitational acceleration constant. n For dynamic contact gap, μ w Let K be the dynamic viscosity coefficient of liquid water. Obtain K. gap Then, combine it with g n Perform multiplication and integral reduction transformation along the interface tangential vector to generate additional water flow replenishment term.

[0092] The original interfacial leakage coefficient was used as the background conductivity, and it was superimposed with the additional water flow recharge term using a matrix addition method, then substituted into the solution equation. Through this superposition mechanism, the parallel flow superposition effect of water flow in micropores and macroscopic fissures was characterized, and the opening calculation of the contact-gated channel was realized.

[0093] Step 303: Extract the interface temperature jump from the generalized jump; Similar to the moisture field processing logic, the generalized jump amount is extracted for the temperature field degree of freedom to obtain the spatial scalar field distribution data that characterizes the temperature difference between the panel side node and the subbase material side node, and this data is used as the interface temperature jump amount.

[0094] Step 304: Determine the conditions for the evolution state and trigger the corresponding thermal impedance calculation path.

[0095] In the first operating condition, when the evolution state is the contact state, the harmonic average thermal conductivity of the solid interface is determined based on the heat transfer properties of the interface material. The harmonic average thermal conductivity of the solid interface is used as the equivalent thermal impedance of the interface at this time, and the interface heat flux feedback term is generated by using the equivalent thermal impedance of the interface and the interface temperature jump.

[0096] Specifically, since the thermal conductivity of the concrete panel is much greater than that of the subbase material, to avoid numerical stiffness singularities caused by excessively large material property differences during equation assembly, the thermal conductivity of both the panel and the subbase material is extracted, their reciprocals are calculated, and the sum of their reciprocals is taken as the harmonized average thermal conductivity of the solid interface. This treatment is equivalent to the contact thermal conductivity calculation in heat conduction theory.

[0097] In the second operating condition, when the evolution state is the voided state, the pre-configured thermal conductivity of the gap-filling medium and the pre-configured interface reference thickness are obtained; combined with the solid interface harmonized average thermal conductivity corresponding to the interface reference thickness and the gap-filling medium thermal conductivity corresponding to the dynamic contact gap, a series equivalent thermal resistance is constructed; the series equivalent thermal resistance is used as the interface equivalent thermal impedance in the voided state, and is used to calculate the interface heat flux feedback term in combination with the interface temperature jump, and the interface heat flux feedback term is used as the thermodynamic component of the bidirectional feedback interface term.

[0098] Specifically, the voided area forms a dielectric interlayer, which blocks the direct solid-phase heat conduction path between the panel and the pad, forming a series thermodynamic loop.

[0099] Series equivalent thermal resistance R eq =(d0 / λ c )+(g n / λ g ); Where d0 is the pre-configured interface reference thickness, i.e., the average feature size of the mesh elements at the non-matching mesh interface, which can be automatically extracted by the finite element preprocessing program, and λ c For the determined harmonic average thermal conductivity of the solid interface, g n For dynamic contact gap, λ g To obtain the thermal conductivity of the gap-filling medium, calculate R. eq Then, its reciprocal is taken to generate the equivalent thermal impedance of the interface, which characterizes the isolation effect of the geometric barrier generated by the de-void on the evolution of the temperature field.

[0100] In another specific embodiment, to obtain the thermal conductivity of the gap-filling medium under the above working conditions, it is necessary to perform dynamic determination of the medium phase state, which specifically includes the following steps: Step 305: Based on the current iteration result of the moisture field equation containing the phase change source, extract the current moisture phase state of the interface and obtain the interface temperature jump amount. Within each time step of the nonlinear iteration, the solver cache variable space is accessed to obtain the pore water content and ice content values ​​of the nodes in the vacuolated region to confirm the water phase state. At the same time, the separated interface temperature jump and absolute temperature scalar are retrieved.

[0101] Step 306: Determine the medium filling type in the current void gap based on the water phase state and the amount of interface temperature jump. The medium filling type includes air filling, liquid water filling and ice filling. In the specific implementation, a multi-logic decision tree is constructed: when the extracted temperature value is consistently lower than the pre-configured freezing point threshold, and the moisture phase indicates that there is accumulated moisture from previous leakage in the area, it is determined to be ice filling; when the temperature is higher than the freezing point threshold, and the local water potential index indicates that there is a sufficient reservoir water supply source, it is determined to be liquid water filling; when the local water head calculation result is close to the residual limit, indicating that the gap is in a state of drainage dryness, it is determined to be air filling.

[0102] Step 307: Based on the identified medium filling category, match the corresponding heat transfer coefficient value from the pre-built medium thermal conductivity property library as the thermal conductivity of the gap filling medium corresponding to the current iteration step.

[0103] Specifically, the thermal conductivity property library is a constant configuration table stored in memory in the form of key-value pairs. According to common physical laws, air has extremely low thermal conductivity, and ice has a much higher thermal conductivity than liquid water.

[0104] Based on the above identification results, index matching is performed. For example, when it is determined to be liquid water filling, the thermal conductivity value is obtained as 0.6; when it is determined to be ice filling, the thermal conductivity value is obtained as 2.2; and when it is determined to be air filling, the thermal conductivity value is obtained as 0.025. The specific numerical value obtained from the matching is assigned to the variable λ. g This ensures that the input parameters of the series equivalent thermal resistance formula fit the current coupled field evolution state.

[0105] This scheme constructs a gated feedback mechanism based on dynamic contact gaps, dynamically mapping the geometric characteristics of physical gaps with hydraulic conductivity and series thermal impedance. This accurately describes the nonlinear conduction effect of macroscopic void channels on reservoir water leakage, as well as the physical blocking effect of interfacial heat transfer caused by the change of media within the gaps. It effectively avoids the numerical flux assessment errors caused by the traditional approach of treating the interface as a fixed, continuous boundary.

[0106] Example 4: Further explanation of the multi-field unified mapping projection calculation and nonlinear simultaneous solution algorithm architecture. Due to the involvement of multi-field strongly coupled calculations at non-matched mesh interfaces, traditional sequential alternating solution methods are prone to numerical oscillations in phase transition regions or abrupt transition points. Therefore, this example adopts a homobasic mapping based on a biorthogonal projection matrix and a semi-smooth Newton algorithm to construct a globally consistent equation solver.

[0107] In one possible implementation, the following steps are included: Step 401: Extract the biorthogonal projection matrix from the mechanical discrete model of the non-matching mesh interface and use it as a unified projection basis operator; Specifically, since the mesh nodes of the panel and the subbase material cannot be perfectly overlapped at the interface, the Mortar contact mechanics method is introduced to extract the integral parameters of the dual shape functions corresponding to the principal and secondary basis functions, constructing a set of transformation matrices that satisfy the biorthogonality property. Utilizing the algebraic characteristics of the diagonalization or quasi-diagonalization of this biorthogonal matrix, high-precision mapping of unmatched boundary field variables can be achieved without the computationally complex global matrix inversion operation. This is then decoupled from the purely mechanical field calculation module and provided as a reusable unified operator for multiphysics calculation processes.

[0108] In some embodiments, the unified projection basis operator includes: determining the secondary surface element as the integration surface; for each secondary surface integration point, searching for its projection position on the principal surface; performing Gauss integration on the contact segment using the shape functions of the secondary surface and the principal surface; assembling the Mortar integration matrices D and M; and simplifying matrix inversion using the diagonalization property of the biorthogonal basis functions.

[0109] Step 402: Input the displacement parameters, temperature parameters and water potential parameters of the face nodes and main face nodes in the multiphysics domain model into the unified projection basis operator for projection calculation, and obtain the displacement jump, interface temperature jump and water potential jump respectively. In the specific calculation process, physical parameters are used as a general column vector of independent variables, and linear algebraic multiplication is performed on them using the pre-generated dual projection matrix. The specific calculation logic is as follows: φ jump =φ slave -D inv *M*φ master ; Where, φ jump φ is the parametric jump of a specified physical field obtained after projection calculation. slave To obtain the physical parameter matrix of the face nodes, D inv M is the inverse matrix of the extracted bioorthogonal basis function integral matrix, M is the master-slave cross basis function integral matrix, and φ is the inverse matrix of the extracted bioorthogonal basis function integral matrix. masterThis is the physical parameter matrix for the main face nodes. By using this general formula, the physical field identifiers are successively replaced with displacement vectors, temperature scalars, and water potential scalars, allowing for lossless calculation of the corresponding displacement jumps, interface temperature jumps, and water potential jumps.

[0110] Step 403: Project the displacement jump along the interface normal to obtain the dynamic contact gap; Obtain the unit normal vector of the mesh interface in the local coordinate system, and perform a vector dot product operation between the displacement jump and the unit normal vector. The resulting scalar result is the dynamic contact gap, which reflects the relative normal opening and closing physical distance of the interface.

[0111] Step 404: Package the displacement jump, interface temperature jump, and water potential jump as a generalized jump output.

[0112] During the memory data processing stage, the three sets of physical jump quantities obtained from the above calculations are assembled into a high-dimensional column vector according to the arrangement rules of the degrees of freedom, thereby forming a set of data blocks containing relative discontinuities in displacement, temperature, and moisture, which are then passed to the subsequent gating feedback judgment module for invocation.

[0113] Based on this, in order to achieve highly convergent solutions for multiphysics, the nonlinear simultaneous iterative solution of the augmented governing equations includes the following steps: Step 405: Reconstruct the preset contact constraint conditions using nonlinear complementary functions to uniformly express the unilateral discontinuity of the interface normal contact state. In traditional contact mechanics, non-penetration conditions are typically presented as a set of inequalities, which conflicts with the format of the governing equations of the finite element system. Therefore, a nonlinear complementary function is introduced to transform this inequality constraint into an equivalent non-smooth equality constraint. The specific reconstruction calculation formula is as follows: C n =p n -max(0,p n -c n *g n )=0; Among them, C n For the residual term of the constructed nonlinear complementary function, p n Let c be the normal contact compressive stress at the interface, max be the function to maximize the stress, and c be the maximum value. n For the pre-configured non-negative penalty parameters, g n This represents the dynamic contact gap obtained from the preceding steps. When the interface becomes vacant, the compressive stress returns to zero; when the interface closes, the gap returns to zero. This formula realizes a unified characterization of the opening and closing evolution logic of the interface state within a single equation.

[0114] Step 406: Calculate the physical parameter partial derivatives of the corresponding variables in the bidirectional feedback interface term and the multiphysics domain model, and extract the interface thermal partial derivative block, interface hydraulic partial derivative block and volume coupling partial derivative block. Using the rules of multivariable calculus, partial derivatives of each parameter in the system equations are calculated with respect to the independent variables. Sub-matrices representing changes in heat flux caused by temperature jumps are extracted from the partial derivative matrices to form interface thermal partial derivative blocks; sub-matrices representing changes in hydraulic conduction caused by gap variations are extracted to form interface hydraulic partial derivative blocks; and partial derivative matrices of the cross terms of the fundamental governing equations for each field within the dam soil are extracted as volume coupling partial derivative blocks.

[0115] Step 407: Merge and assemble the nonlinear complementary function, the interface thermal partial guide block, the interface hydraulic partial guide block, and the volume coupling partial guide block to generate the global tangent stiffness matrix. Based on the superposition rule of the global stiffness matrix in the finite element method, the various partial derivative matrix blocks generated above are assembled into corresponding positions in a unified large sparse matrix according to the global node index. This operation concentrates all nonlinear coupling behaviors into a single matrix system.

[0116] Step 408: Based on the global tangent stiffness matrix and the semi-smooth Newton method, the augmented governing equations are updated and uniformly monotonically converged, and the contact activation set of the current iteration step is output simultaneously.

[0117] The generated tangent stiffness matrix is ​​inverted and its increment is updated using a numerical algebra solver. Since the nonlinear complementary function contains non-differentiable extremum terms, a semi-smooth Newton method is used to construct the generalized Jacobian matrix for solution. During this process, the integration points in the current iteration step that are in a closed-loop bearing state are extracted and added to a dynamic set list, defined as the contact activation set.

[0118] Furthermore, to prevent spurious convergence caused by repeated jumps in interface states during the solution process, the iterative convergence exit conditions for the nonlinear simultaneous iterative solution of the augmented control equations include: Step 409: Obtain the multi-field residuals and relative values ​​of displacement increments for the current iteration step; Step 410: When the relative value of the displacement increment and the multi-field residual are both lower than the pre-configured tolerance threshold, and the contact activation set remains unchanged in two consecutive iterations, the overall iteration convergence is confirmed, and the simultaneous iteration solution is terminated.

[0119] Specifically, a dual-control mechanism is employed to verify the reliability of the calculation results. The first layer of control requires that the update magnitudes of the mechanical, temperature, and seepage residual vector norms, as well as the solution, be less than a given accuracy requirement. The second layer of control requires performing a logical comparison, comparing the contact activation set list of the current iteration step with the list content of the previous historical iteration step. Only when the node index numbers contained in the two sets of lists are completely identical is it confirmed that the spatial topology of the vacant region has stabilized, allowing the solver to exit the inner nonlinear loop of the current time step.

[0120] In some optional implementations, for situations where the dam's state rapidly evolves due to sudden changes in natural temperature or rapid water storage, an adaptive time step strategy is employed for time evolution control when solving the augmented control equations nonlinearly and iteratively. This specifically includes: Step 411: Obtain the actual number of iterations for the previous time step and the pre-configured target number of iterations; Step 412: Based on the deviation ratio between the actual number of iterations and the target number of iterations, dynamically scale the current time step to obtain the updated time step; In the specific implementation process, the time increment is dynamically adjusted proportionally according to the convergence efficiency of Newton's iteration method. The calculation logic is as follows: dt next =dt current *(N target / N actual ) γ ; Among them, dt next dt represents the next time step obtained after scaling calculation. current N represents the current time step being used. target N is the pre-configured target number of iterations. actual γ represents the actual number of iterations recorded by the system monitoring system, and γ is a pre-configured smoothing scaling exponent constant. When the actual number of iterations is too high, the deviation ratio is less than the constant 1, and the system will automatically reduce the time step to reduce the nonlinear intensity.

[0121] Step 413: During the solution process, when the contact activation set is detected to have a high frequency (more than a preset number of) state switching in continuous iterations, the time step shrinking mechanism is triggered to reduce the updated time step and use the reduced time step for the solution calculation of the next time step to ensure the convergence stability of the strongly coupled solution.

[0122] Specifically, if the contact activation set oscillates between adding and removing nodes in three or more consecutive nonlinear iterations within the same time step, it indicates that the current load increment exceeds the analytical capability of the interface gating feedback response. In this case, the conventional dynamic scaling logic is blocked, and the current time step is directly forced to be halved. The state is then rolled back and the equations are recalculated to prevent numerical divergence in the evolution of the physical field at such abrupt critical points.

[0123] This scheme utilizes asymmetric unified projection basis operators to perform high-fidelity cross-plane mapping. This method can extract generalized physical jumps without complex global interpolation corrections, thereby ensuring that the energy and mass of multiple field parameters such as heat, water, and force meet the conservation requirements when experiencing discontinuous abrupt changes in spatial scale.

[0124] Example 5: In the previous examples, the Lagrange multipliers were introduced as independent degrees of freedom into the global equations to apply contact constraints precisely. However, this would result in the global tangent stiffness matrix having zero principal component regions (i.e., the matrix is ​​not positive definite), thereby increasing the risk of divergence in the numerical solver's inversion calculation.

[0125] To eliminate this potential problem, this embodiment provides a local basis transformation alternative method based on the multiplier dimensionality reduction decoupling mechanism.

[0126] Specifically, for the preset contact constraints and the corresponding contact mechanics calculations, this method includes the following steps: Step 501: Establish the local projected coordinate system of the master and slave surface nodes using contact constraints, and define unconstrained generalized displacement variables. In a specific implementation, the normal and tangential geometric topological information of the mismatched mesh interface is extracted. Based on the displacement compatibility constraint equation where the normal cannot penetrate, the constrained slave node displacements are separated into two parts: one part is a non-independent component dominated by the master node displacement, and the other part is an independent slip / free-slip component that slides along the tangential direction or separates along the normal direction. This independent degree of freedom is combined with the unconstrained degree of freedom of the master node and defined as a reduced-order column vector, which is used as the unconstrained generalized displacement variable. Since the non-independent parameters are extracted, the dimension of this generalized displacement variable is necessarily smaller than the original node displacement dimension of the system.

[0127] Step 502: Construct a global transformation matrix based on the mapping relationship between the generalized displacement variables and the original nodal displacements; In this step, a mapping operator is constructed based on the kinematic dependency constraints established in the previous step. Specifically, the partial derivative coefficients in the master-slave node constraint equations are extracted, and through matrix transpose and recombination operations, a non-square projection matrix of dimension N rows by M columns is generated (where N is the total number of original node displacements and M is the total number of generalized displacement variables). This matrix is ​​the global transformation matrix. This matrix has the mathematical transfer function of mapping a low-dimensional unconstrained space to a high-dimensional physical coordinate space with contact constraints.

[0128] Step 503: Perform an orthogonal transformation operation on the original stiffness matrix of the augmented governing equations, and use the global transformation matrix to generate an equivalent unconstrained stiffness matrix to eliminate the contact force Lagrange multipliers. The calculation logic of the orthogonal transformation operation is to multiply the transpose of the global transformation matrix, the original stiffness matrix and the global transformation matrix in sequence. Specifically, the original augmented governing equations, due to the inclusion of Lagrange multipliers corresponding to contact forces, exhibit a typical saddle-point equation structure, making them difficult to process by conventional iterative solvers. Therefore, an orthogonal transformation operation is performed to algebraically remove the Lagrange multipliers from the solution objective. The specific matrix transformation calculation formula is as follows: K reduced =P T *K origin *P; Among them, K reduced For the generated equivalent unconstrained stiffness matrix, P T K is the transpose of the global transformation matrix. origin Let P be the original stiffness matrix of the augmented governing equations containing Lagrange multiplier degrees of freedom, and let P be the global transformation matrix. Through the above multiplication calculations, the saddle point matrix, originally containing multiplier degrees of freedom, is projected into a reduced subspace that does not contain constraint multipliers. Since the constraints are implicitly contained in the matrix operator P, the generated K... reduced It restores the standard symmetric positive definite or semi-positive definite properties.

[0129] Step 504: After solving the displacement principal variables based on the equivalent unconstrained stiffness matrix, the interface contact force is re-extracted as the updated interface gating feedback input through the reverse recovery step.

[0130] In this step, a numerical solver is invoked to invert or iteratively solve the scaled equation system generated in step 503, directly obtaining the unconstrained generalized displacement components. Subsequently, inverse algebraic operations are performed, multiplying the obtained generalized displacements by the aforementioned global transformation matrix to reconstruct the full-scale original nodal displacement distribution of the dam. After displacement reconstruction, small residual terms that do not satisfy the contact compatibility equations are extracted, and internal forces are inversely calculated using the system's interface stiffness coefficients. The contact compressive stress distribution data in the interface normal direction is extracted and used as the physical forces corresponding to the Lagrange multipliers, then reinjected into the subsequent multiphysics interface-gated feedback calculation loop. This alternative scheme preserves the mathematical rigor of the constraint equations while avoiding the risk of computational collapse due to the non-positive definiteness of the stiffness matrix.

[0131] According to one aspect of this application, when performing secondary development using a commercial engineering software platform, in order to be compatible with existing solver architectures and improve the convergence stability of large matrices, an augmented Lagrangian method can be used to construct a nonlinear solution mechanism for contact states. Specifically, this includes the following steps: Step 601: Define the concrete panel mesh surface as the master surface and the subbase material mesh surface as the slave surface to generate a master-slave contact pair; Specifically, after discretizing the solid dam interface in finite element space, the outer boundary node sets of both sides of the interface are extracted. Based on the differences in elastic modulus and stiffness of the materials, the boundary node set of the concrete panel with high stiffness and dense mesh is selected as the master surface, and the boundary node set of the relatively loose and less stiff subgrade material is selected as the slave surface. Using this geometric topological relationship, a constraint pairing relationship for searching relative normal distances is constructed in memory, forming master-slave contact pairs. The purpose of this definition is to solidify the reference surface for contact search on a surface with high rigidity, thereby reducing the geometric penetration identification error when determining the contact state.

[0132] Step 602: Introduce an augmented Lagrangian function with a penalty factor between the master-slave contact pairs to transform the physical impenetrability condition into an energy functional variational problem with a penalty term. In this step, the pure Lagrange multiplier method introduces additional unknowns and violates the positive definiteness of the system, while the pure penalty function rule suffers from numerical defects such as generating ill-conditioned matrices and uncontrollable penetration. This scheme integrates the two, adding both contact penalty functional terms and multiplier work terms to the overall potential energy control equation of the system. The specific mathematical update logic for contact compressive stress is as follows: T n =λ n +ε n *g pen ; Among them, T n λ is the contact compressive stress calculated in the current iteration step. nThe Lagrange multiplier components saved in the previous outer iteration step, ε n For the introduced penalty factor, g pen This represents the currently calculated geometric penetration amount. In the standard parameter configuration, to ensure sufficient initial constraint stiffness, the penalty factor ε is set... n The initial value is pre-configured as 10. 9 By constructing the augmented function described above, and utilizing the variational extremum principle, the absolute impenetrability constraint in physical space is transformed into a differentiable partial differential penalty term within the energy functional of the governing equations.

[0133] Step 603: Using the contact solver built into the commercial engineering software, execute the alternating update algorithm of Lagrange multipliers and penalty factors until the interface intrusion is lower than the set mesh tolerance threshold, thereby generating a stiffness matrix block equivalent to applying rigid contact constraints.

[0134] In the specific implementation, a double-nested iterative logic is used to solve the above functional equations. In the outer loop, the currently stored Lagrange multiplier variables are fixed, and the global displacement vector is solved using the Newton-Raphson iteration method. After the displacement solution is completed, the inner loop is entered to extract the actual interface intrusion amount of the contact interface nodes and determine whether it is lower than the mesh tolerance threshold. Here, the mesh tolerance threshold is usually pre-configured as 1% of the minimum feature size of the mesh at the contact interface. If the intrusion amount is determined to be greater than this tolerance, the Lagrange multipliers are updated using the formula in step 602, and the penalty factor is dynamically amplified based on the set amplification factor. Then, the outer loop is forcibly returned to redistribute the displacement field. When it is detected that the interface intrusion amount of all face nodes meets the tolerance threshold, the alternating update loop is stopped, and the tangent stiffness increment matrix caused by the contact nonlinearity is extracted and output as a stiffness block equivalent to the application of rigid contact constraints to the global equation system.

[0135] In some alternative implementations, when a native contact module is lacking and a partial differential equation compiler is used to perform the simulation, the preset contact constraints can also be implemented through weak-form substitution of partial differential equations. Specifically, the nonlinear complementary function residual term from Example 4 is extracted and algebraically rewritten using conditional logic instructions supported by the programming language or smoothed step function instructions. After rewriting, it is directly incorporated as an independent boundary flux integral term or intrinsic term into the weak-form residual integral expression of the global stiffness matrix. The residual term is implicitly processed by the underlying linear equation solver, which can also realize the opening and closing abrupt change logic operation of the interface, thereby avoiding the geometric surface search process of complex contact pairs.

[0136] By reconstructing the unilateral contact constraints of the interface using nonlinear complementary functions, a globally consistent nonlinear strongly coupled solution framework is established. Simultaneously, by combining a time-step adaptive algorithm based on physical residuals with a low-level constraint dimensionality reduction and decoupling mechanism, the potential for matrix nondefiniteness is eliminated from the algebraic equation structure, enhancing the model's computational robustness under extreme loading and alternating phase transition conditions.

[0137] like Figures 2 to 7 As shown, this method can clearly reveal the discontinuous temperature field distribution and concentrated moisture leakage paths behind the voided area of ​​the panel. Compared with traditional node binding methods, the temperature gradient and displacement distribution at the interface output by this method better reflect the impact of changes in contact state on the multiphysics conduction path in actual engineering.

[0138] Those skilled in the art should understand that, in practical applications, some of the preset parameters, threshold conditions, and inherent attribute values ​​involved in the embodiments of the present invention can be obtained or adjusted by engineers through conventional technical means, experimental statistics, factory calibration, or by directly consulting data manuals and other existing methods. For example, the inherent physical properties of materials such as liquid water density, gravitational acceleration constant, porosity, and the thermal conductivity coefficient of the panel and subbase materials can be directly obtained by consulting standard physical data handbooks, material factory calibration documents, or through conventional geotechnical testing. Model evolution correction parameters such as ice resistance empirical parameters, anisotropic frost heave distribution coefficient, cumulative frost heave coefficient, and reference constraint compressive stress can be obtained through statistical fitting based on indoor freeze-thaw cycle mechanical test data, consulting historical experience databases of geotechnical engineering in high-altitude and cold regions, or calculation according to relevant water conservancy industry standards. Algorithm control parameters and business discrimination indicators such as freezing temperature threshold, target iteration number for simultaneous iterative solutions, smoothing scaling exponential constant, initial augmentation penalty factor, mesh tolerance threshold, and structural strength allowable threshold for risk warning can be adaptively evaluated and calibrated by engineering designers based on the specific dam project scale, finite element mesh discretization accuracy level, dam safety operation design specifications, and convergence and debugging experience from numerical simulation calculations. The specific values ​​of these parameters will be adaptively adjusted according to changes in the physical characteristics of the dam, the evolution of meteorological environmental conditions, and the selection of software and hardware computing platforms. The process of obtaining and adjusting these parameters is common knowledge to those skilled in the art and should not be regarded as a limitation on the integrity and scope of protection of the core technical solution of this invention.

Claims

1. A method for calculating the mismatch deformation between the subbase material and the face panel of a high-altitude rockfill dam in cold regions, characterized in that, include: Obtain the structural parameters and operational boundary conditions of a high-face rockfill dam; A multiphysics domain model with a non-matching mesh interface is established based on structural parameters. The multiphysics domain model includes a moisture field equation with a phase change source, a temperature field equation with latent heat, and a displacement field equation based on effective stress. Construct a unified projection basis operator based on contact mechanics at the interface of the mismatched mesh; By using the unified projection basis operator to perform cross-interface mapping of the physical parameters of displacement field, temperature field and moisture field, the generalized jump quantity and dynamic contact gap of the interface are obtained. The evolution state of the interface is determined based on the dynamic contact gap; Gated feedback calculations are performed based on the evolutionary state and generalized jump amount to obtain the interface equivalent thermal impedance and interface hydraulic conductivity, so as to construct a bidirectional feedback interface term. Augmented control equations are constructed based on a multiphysics domain model, bidirectional feedback interface terms, and preset contact constraints. The augmented control equations are solved nonlinearly and iteratively to output the multiphysics distribution and incoordination deformation state of the high panel rockfill dam.

2. The method according to claim 1, characterized in that, The water field equations containing phase change sources include: Obtain the volumetric ice content and unfrozen water content of the soil; Based on the volumetric ice content, an exponential decay function is constructed to characterize the ice blockage effect, serving as the ice resistance factor. Based on the unfrozen water content, the matrix potential-volume water content relationship function and the volumetric water content-permeability coefficient relationship function established by the VanGenuchten model were used to calculate the foundation relative permeability coefficient, and the ice resistance factor was used to correct the foundation relative permeability coefficient to obtain the equivalent relative permeability coefficient. The head boundary conditions are determined based on the reservoir water level time series in the operational boundary conditions. The equivalent relative permeability coefficient and the head boundary conditions are then substituted into the water field equation containing the phase change source to perform saturated-unsaturated seepage calculations.

3. The method according to claim 1, characterized in that, The displacement field equation based on effective stress includes the calculation of the frost heave effect on the subbase material. The specific process includes: Obtain the normal effective constraint compressive stress and the pre-configured reference constraint compressive stress at the interface of the subbase material; Based on the relative relationship between the effective normal constraint compressive stress and the reference constraint compressive stress, an exponential suppression function is constructed to characterize the gravitational constraint suppression effect. The volumetric ice content of the soil is obtained, and the volumetric ice content is converted into free frost heave strain in combination with the pre-configured frost heave coefficient. The free frost heave strain is attenuated and corrected using an exponential suppression function to obtain the confined frost heave strain, which is then incorporated into the displacement field equation based on the effective stress.

4. The method according to claim 1, characterized in that, Construct a unified projection basis operator based on contact mechanics at the interface of the mismatched mesh, including: The bioorthogonal projection matrix is ​​extracted from the mechanical discrete model of the non-matching mesh interface and used as a unified projection basis operator. Accordingly, the physical parameters of the displacement, temperature, and moisture fields are mapped across the interface using unified projection basis operators to obtain the generalized jump quantity and the dynamic contact gap of the interface, including: The displacement parameters, temperature parameters, and water potential parameters of the face nodes and the main face nodes are respectively input into the unified projection basis operator for projection calculation, and the displacement jump, interface temperature jump, and water potential jump are obtained accordingly. The dynamic contact gap is obtained by projecting the displacement jump along the interface normal. The displacement jump, interface temperature jump, and water potential jump are packaged together and output as a generalized jump quantity.

5. The method according to claim 1, characterized in that, The pre-defined contact constraints in the augmented governing equations are implemented using the multiplier transformation method, specifically including: Define the non-penetration condition in the interface normal and the frictional sliding condition in the tangential direction; Based on local basis transformation, a generalized displacement is defined for interface nodes in the active contact state; Construct a global transformation matrix based on generalized displacement; By using a global transformation matrix, the Lagrange multiplier degrees of freedom introduced by solving the contact constraints are transformed into equivalent constraint terms, eliminating the Lagrange multipliers, completing the point-by-point decoupling of the augmented control equations, and satisfying the non-penetration condition and the friction sliding condition.

6. The method according to claim 1, characterized in that, The pre-defined contact constraints in the augmented governing equations are implemented using the augmented Lagrangian method, specifically including: Define the master face and slave face on both sides of the non-matching mesh interface to construct interface contact pairs; Set the initial value of the pre-configured penalty factor and the allowable penetration amount; In the nonlinear iterative solution of the augmented control equations, based on the interface contact pairs, the interface approximates the non-penetration condition and friction constraint condition through the alternating update calculation of the penalty factor and Lagrange multiplier, thereby achieving control of the contact constraints.

7. The method according to claim 1, characterized in that, The augmented governing equations are solved nonlinearly through simultaneous iterative solutions, including: The preset contact constraints are reconstructed using nonlinear complementary functions to uniformly express the unilateral discontinuity of the interface normal contact state. The physical parameter partial derivatives of the corresponding variables in the bidirectional feedback interface term and the multiphysics domain model are calculated respectively, and the interface thermal partial derivative block, interface hydraulic partial derivative block and volume coupling partial derivative block are extracted. The nonlinear complementary function, the interface thermal partial guide block, the interface hydraulic partial guide block and the volume coupling partial guide block are merged and assembled to generate the global tangent stiffness matrix. Based on the global tangent stiffness matrix and the semi-smooth Newton method, the augmented governing equations are updated and uniformly monotonically converged, and the contact activation set of the current iteration step is output simultaneously.

8. The method according to claim 7, characterized in that, The iterative convergence exit conditions for nonlinear simultaneous iterative solutions of augmented governing equations include: Obtain the multi-field residuals and relative values ​​of displacement increments for the current iteration step; When the relative value of the displacement increment and the multi-field residual are both lower than the pre-configured tolerance threshold, and the contact activation set remains unchanged in two consecutive iterations, the overall iteration convergence is confirmed, and the simultaneous iteration solution ends.

9. The method according to claim 8, characterized in that, When solving the augmented control equations nonlinearly using simultaneous iterative methods, an adaptive time step strategy is employed for time evolution control, specifically including: Obtain the actual number of iterations at the previous time step, and the pre-configured target number of iterations; The current time step is dynamically scaled based on the deviation ratio between the actual number of iterations and the target number of iterations. During the solution process, when a state switch is detected in the contact activation set after more than a preset number of consecutive iterations, a time step shrinking mechanism is triggered to reduce the current time step in order to ensure the convergence stability of the strongly coupled solution.

10. The method according to claim 1, characterized in that, After outputting the multiphysics distribution and incoordination deformation state of the high-face rockfill dam, the following is also included: Based on the output of the non-coordinated deformation state, the spatial distribution range of the panel void area and the interfacial frost heave pressure of the padding material on the panel are extracted. Based on the pre-stored structural strength threshold and interfacial frost heave pressure of the high-face rockfill dam, a cracking risk assessment of the concrete face is conducted. The spatial distribution ranges of cracking risk areas and voided areas are spatially superimposed and mapped to calculate the proportion of voided area and / or stress concentration coefficient as dam operation safety evaluation indicators, and an early warning is issued when the safety evaluation indicators exceed the preset safety warning threshold.

Images