A method for simulating thermodynamics of dynamic terrain based on fluid-structure coupling framework

Abstract

The application discloses a dynamic topography thermodynamic simulation method based on a fluid-solid coupling framework, and comprises the following steps: at the beginning of each time step, temperature extrapolation is performed by using temperature fields of two previous time points to obtain a predicted temperature field; the predicted temperature field is used as a thermodynamic initial value, fluid equations, viscosity updating, solid equations and thermodynamic equations are sequentially solved in the same Newton framework, and velocity continuity is applied on a fluid-solid interface to form a nonlinear block coupling iteration; when residual convergence is achieved, a grid and a field quantity are updated according to the converged solution, and a step length is adaptively determined, and the cycle is repeated until the end. By establishing a heat-force-flow single-block coupling variation model under a unified arbitrary Lagrange-Euler framework and jointly solving, the physical fidelity of mantle convection driven rock circle deformation simulation is improved, and the problem that a fluid-solid-thermal closed-loop coupling feedback cannot be accurately described is solved.

Description

Technical Field

[0001] This invention relates to the field of data simulation technology, and in particular to a dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework. Background Technology

[0002] In fields such as geodynamics, geological engineering, and earthquake prediction, accurately simulating the long-term deformation, thermal convection, and faulting activity of the Earth's crust and mantle over millions of years is of fundamental scientific significance for understanding major geological processes such as plate tectonics, orogeny, basin evolution, and periodic seismic activity. The Earth's internal materials exhibit complex fluid-structure interaction and thermodynamic coupling characteristics: mantle convection is essentially a thermal convection process of high-viscosity fluids driven by temperature gradients, while the crust and lithosphere exhibit elastic or elastoplastic solid behavior under long-term tectonic stress. The interaction between the two—that is, mantle convection driving lithosphere deformation, and lithosphere deformation in turn influencing mantle convection patterns and heat transport processes—is the fundamental physical mechanism driving core geological phenomena such as plate tectonics, orogeny, and seismic activity. Therefore, establishing a numerical simulation framework capable of uniformly describing the aforementioned fluid-structure interaction, thermal convection, and complex material rheological behavior has long been a core goal pursued in the field of computational geodynamics.

[0003] In existing numerical simulation methods, early geotectonics simulation codes, such as the SOPALE program developed by Fullsack in 1995, employed the arbitrary Lagrange-Euler method to solve the Stokes equations under the "creep flow" limit, achieving the simulation of large deformation behavior of viscoplastic materials. However, such methods typically treat the crust and mantle as fluid media with different high viscosities, making it difficult to accurately characterize the strong coupling effects at the fluid-solid interface, especially the Mohorovičić discontinuity between the crust and mantle, and also failing to fully describe the elastic or elastoplastic finite deformation behavior of the solid lithosphere. In recent years, significant progress has been made in numerical methods for fluid-solid coupling based on the arbitrary Lagrange-Euler framework, exemplified by the single-phase fluid-solid coupling solver implemented by Wick in the deal.II finite element library in 2013. This method uses biharmonic equations to construct fluid mesh motion to handle large deformations and solves strongly coupled nonlinear systems using Newton's method, demonstrating good applicability in fluid-solid coupling problems in biomedical and engineering fields. However, the design goals of this solution framework are mainly for biomedical and engineering applications. It lacks specific optimization for the long-term, large-scale three-dimensional simulation requirements of geodynamic problems, and it does not integrate thermodynamic equations, so it cannot describe the mutual feedback between the temperature field and the fluid-structure interaction process.

[0004] In thermodynamic simulations, the step-32 tutorial program in the deal.II finite element library is a representative example, demonstrating how to solve the thermal convection equations describing mantle convection on a large-scale parallel basis using the Boussinesq approximation. This program can effectively simulate the evolution of the temperature and velocity fields of mantle thermal convection and has good parallel scalability. However, it is not coupled with a solid deformation module, simplifying the lithosphere as the upper boundary condition of the convection system rather than a deformable solid. Therefore, it cannot fully describe the closed-loop coupled process of mantle convection driving lithosphere deformation and lithosphere deformation influencing mantle convection.

[0005] In summary, the main shortcomings of existing technologies are as follows: Firstly, existing fluid-structure interaction (FSI) solution frameworks lack deep integration with thermodynamic modules, making it impossible to establish a unified numerical model with complete thermo-mechanical-fluid coupling. Secondly, existing mantle thermal convection solvers are not coupled with solid mechanics modules, making it difficult to simulate key physical processes such as thermal expansion, thermal stress, and the feedback of deformation on heat transport. These shortcomings result in the current lack of a dynamic topographic thermodynamic simulation method that can solve fluid dynamics, solid mechanics, and thermodynamic equations in a single-block coupled manner within a unified arbitrary Lagrange-Euler variational framework, making it difficult to simulate the complete physical process of mantle convection driving lithospheric deformation with high accuracy and efficiency. Summary of the Invention

[0006] The purpose of this invention is to provide a dynamic topographic thermodynamic simulation method based on a fluid-structure interaction framework. By establishing a single-block coupled variational model of the thermo-mechanical-fluid multiphysics field under a unified arbitrary Lagrange-Euler framework and solving it jointly, the physical fidelity of the simulation of mantle convection-driven lithospheric deformation is improved. This solves the technical problem that existing technologies cannot accurately characterize the closed-loop coupled feedback between fluid motion, solid deformation and heat transport in the same numerical framework.

[0007] To address the aforementioned technical problems, a first aspect of this invention provides a method for simulating dynamic terrain thermodynamics based on a fluid-structure interaction framework, comprising the following steps: Step S1: Load the initial temperature field of the target geological area. At the beginning of each time step, perform temperature extrapolation using the temperature fields of the previous two known times to obtain the predicted temperature field of the current time step. The temperature extrapolation is an operation that uses the temperature field at the known time as input and constructs an approximate value of the temperature field at the current time through time extrapolation. Step S2: Using the predicted temperature field as the initial thermodynamic value, under the same Newton iteration framework, sequentially solve the momentum and mass conservation equations in the fluid domain, update the viscosity parameters, solve the finite deformation equations in the solid domain, and solve the global thermodynamic equations. Apply the velocity continuity condition at the fluid-solid interface to form a nonlinear block coupling iteration that solves all physical field variables jointly. Step S3: Determine whether the residual of the current iteration has converged. If it has not converged, adjust the Newton iteration step size and return to step S2 to continue iterating. If it has converged, update the mesh and field quantity according to the converged solution, adaptively determine the step size of the next time step, perform velocity extrapolation and enter the next time step loop until the simulation ends, and output the dynamic terrain evolution results.

[0008] Furthermore, the joint solution of all physical field variables described in step S2 includes: S21: Identify the block sparse structure of the algebraic equation system in the nonlinear block coupling iteration, and divide the unknowns into the first variable group and the second variable group according to the physical field type; S22: Construct the Schur complement matrix for the divided variable groups to form a reduced system with the first variable group as the unknown; S23: Solve the reduced system to obtain the correction amount for the first set of variables; S24: Substitute the correction amount of the first variable group back into the second variable group to calculate the correction amount of the second variable group; S25: Update all physical field variables together with the correction values ​​of the first variable group and the second variable group to complete this joint solution.

[0009] Further, the step S21 of identifying the block sparse structure of the algebraic equation system in the nonlinear block coupling iteration includes: S211: The algebraic equations are divided into blocks according to the four physical field variables of fluid velocity, pressure, solid displacement and temperature to obtain a block sparse matrix; S212: Identify the symmetric positive definiteness of the fluid velocity submatrix and the block diagonality of the solid displacement submatrix in the block sparse matrix, wherein the symmetric positive definiteness is used to guide the selection of a geometric multigrid preprocessor and the block diagonality is used to guide the selection of an algebraic multigrid preprocessor.

[0010] Furthermore, the method of classifying by physical field type in step S21 is as follows: fluid-solid interface displacement and temperature are used as the first variable group, and fluid velocity, pressure and solid displacement other than interface displacement are used as the second variable group. This allows the key physical quantities controlling fluid-solid coupling and thermo-mechanical coupling to be assigned to the first variable group, and the subordinate physical quantities that can be driven by the first variable group through coupling relationship to be assigned to the second variable group.

[0011] Furthermore, in step S221, the submatrix of the second variable group is constructed using a hybrid multigrid preprocessing method to approximate its inverse: S221a: Apply geometric multigrid preprocessing to the sub-matrix block corresponding to fluid velocity and pressure to obtain the approximate inverse of the sub-matrix block; S221b: Apply algebraic multigrid preprocessing to the sub-matrix block corresponding to the solid displacement to obtain the approximate inverse of the sub-matrix block; S221c: Assemble the approximate inverses obtained in the two steps according to the block diagonal structure to form the overall approximate inverse of the second variable group submatrix.

[0012] Furthermore, the adaptive determination of the step size for the next time step in step S3 includes: S31: Calculate candidate step sizes based on numerical stability conditions; S32: Determine the order of magnitude range of the current simulation time, apply logarithmic time scale adjustment to the candidate step size according to the order of magnitude range, and obtain the final next time step size to be adopted; S33: Update the current time step with the final adopted next time step size, and use it as the time advance step size for the next time step cycle.

[0013] Further, step S32, which involves determining the order-of-magnitude interval of the current simulation time and applying a logarithmic time scale adjustment to the candidate step size based on the order-of-magnitude interval to obtain the final adopted next time step size, includes: S321: Divide the total simulation duration into multiple order-of-magnitude intervals; S322: Determine the order of magnitude interval of the current simulation time. When the simulation process moves from the previous order of magnitude interval to the next order of magnitude interval, adjust the candidate step size according to the preset interval step size mapping relationship so that the logarithm of the adjusted step size is linearly related to the simulation process, and obtain the final next time step size.

[0014] Furthermore, in step S2, the momentum and mass conservation equations for the fluid domain are adopted as the incompressible Navier-Stokes equations approximated by Boussinesq, and the finite deformation equations for the solid domain are adopted as hyperelastic or viscoelastic-plastic constitutive equations described by the Piola-Kirchhoff stress tensor. In step S3, the mesh motion in the updated mesh and field quantities is controlled by a biharmonic equation, which is a mesh motion equation with the composite of the Laplace operator as the control operator.

[0015] Furthermore, the momentum and mass conservation equations in the fluid domain, the finite deformation equations in the solid domain, and the global thermodynamic equations are constructed as a unified variational weak form under the reference configuration. The unified variational weak form uses fluid velocity, pressure, solid displacement, and temperature as unknowns and includes fluid momentum terms, fluid convection terms, fluid stress terms, pressure divergence terms, solid inertia terms, solid stress terms, thermal transient and convection terms, thermal diffusion terms, as well as buoyancy terms, volume force terms, and internal heat source terms. The velocity continuity and stress equilibrium conditions at the fluid-solid interface are satisfied in the unified variational weak form.

[0016] Furthermore, during the time progression in steps S1 to S3, a multi-mechanism numerical stabilization process is applied, which includes: discretizing the time derivative term using a backward second-order difference scheme in step S2; performing pressure regularization when solving the fluid domain equation in step S2; and applying artificial viscosity based on entropy viscosity when solving the thermodynamic equation in step S2.

[0017] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework.

[0018] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework.

[0019] The above-described technical solutions of the embodiments of the present invention have the following beneficial technical effects: 1. By identifying the block sparse structure of multi-physics coupled algebraic equation systems, the unknowns are divided into a first variable group containing fluid-solid interface displacement and temperature and a second variable group containing fluid velocity, pressure and other solid displacements according to the type of physical field. Geometric multigrid and algebraic multigrid preprocessing are applied to the sub-matrix blocks corresponding to different physical fields, which significantly accelerates the convergence speed of large-scale nonlinear systems and achieves higher computational efficiency than general solvers. This solves the problem of low computational efficiency caused by the failure to fully utilize the system matrix structure characteristics when dealing with fluid-solid-thermal three-way coupling problems in existing technologies. 2. By constructing the fluid dynamics equations, solid mechanics equations, and thermodynamics equations into a single-block coupled variational model with fluid velocity, pressure, solid displacement, and temperature as unknowns within a unified arbitrary Lagrange-Euler framework, the velocity continuity and stress equilibrium conditions at the fluid-solid interface are automatically satisfied in this model. This enables a closed-loop coupled simulation of thermal driving, force transmission, and deformation feedback during mantle convection-driven lithospheric deformation, fundamentally overcoming the shortcomings of existing technologies that cannot fully characterize the thermo-mechanical-fluid feedback mechanism due to separating fluid, solid, and thermodynamic processes. 3. By combining candidate step size calculation based on numerical stability conditions with logarithmic time scale adjustment during the time progression, the step size is adaptively adjusted according to the order of magnitude range of the simulation process. At the same time, backward second-order difference scheme time discretization, pressure regularization processing and artificial viscosity processing based on entropy viscosity are applied to ensure the numerical stability and accuracy of the integration process with a time span of millions of years. This solves the technical problem that existing technologies cannot balance long-term integration efficiency and numerical stability of multi-physics coupling. Attached Figure Description

[0020] Figure 1 This is a flowchart of the dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework provided in this embodiment of the invention; Figure 2 This is a logical schematic diagram of the dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework provided in an embodiment of the present invention; Figure 3 This is an adaptive time step control flowchart of the dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework provided in this embodiment of the invention. Figure 4 This is a logical architecture diagram of the fluid-structure-thermal coupling solution for the dynamic terrain thermodynamic simulation method based on the fluid-structure coupling framework provided in this embodiment of the invention. Detailed Implementation

[0021] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments and the accompanying drawings. It should be understood that these descriptions are merely exemplary and not intended to limit the scope of the invention. Furthermore, descriptions of well-known structures and techniques are omitted in the following description to avoid unnecessarily obscuring the concept of the invention.

[0022] Please refer to Figure 1 The first aspect of this invention provides a dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework, comprising the following steps: Step S1: Load the initial temperature field of the target geological region. At the beginning of each time step, perform temperature extrapolation using the temperature fields of the previous two known times to obtain the predicted temperature field for the current time step. Temperature extrapolation is an operation that uses the temperature fields at known times as input and constructs an approximate value of the temperature field at the current time through time extrapolation.

[0023] Starting with the initial temperature field of the target geological region, the method determines the initial thermal state for the entire long-term dynamic evolution calculation. At the beginning of each time step, the method first performs a temperature extrapolation operation. That is, using the temperature field information from the previous two known time steps, an approximate value of the temperature field at the current time step is constructed through time extrapolation. This approximate value is used as the predicted temperature field for the current time step. This prediction step provides the necessary and reasonable thermodynamic initial conditions for the subsequent coupled solution process, ensuring that the temperature field already has a distribution pattern that matches the evolution trend of the previous time step before the formal solution begins. This reduces the initial value deviation in the nonlinear iteration process and helps to improve the overall solution efficiency.

[0024] Step S2: Using the predicted temperature field as the initial thermodynamic value, under the same Newton iteration framework, sequentially solve the momentum and mass conservation equations in the fluid domain, update the viscosity parameters, solve the finite deformation equations in the solid domain, and solve the global thermodynamic equations. Apply the velocity continuity condition at the fluid-solid interface to form a nonlinear block coupled iteration that solves all physical field variables jointly.

[0025] After obtaining the predicted temperature field, the core multiphysics coupling solution process begins. Using the predicted temperature field as the initial thermodynamic value, four key solution actions are executed sequentially within the same Newtonian iterative framework: solving the momentum and mass conservation equations for the fluid domain to obtain the velocity and pressure fields of the fluid region; updating the viscosity parameters of the fluid and solid materials based on the currently calculated velocity field, strain rate, and temperature distribution to reflect the rheological properties of deep Earth rocks under high temperature and high pressure conditions; solving the finite deformation equations for the solid domain to obtain the displacement field of the lithosphere under the combined action of gravity load and fluid drag force; and solving the global thermodynamic equations, covering heat transport in the mantle convection region, heat conduction within the crust, and internal heat source effects such as shear heat generation and radioactive decay. During the above solution process, a velocity continuity condition is applied at the fluid-solid interface, ensuring a strict match between the fluid velocity and the solid deformation rate at the interface, guaranteeing the coordinated transfer of kinematic quantities. The four solution actions and interface conditions described above together constitute a nonlinear block-coupled iteration that jointly solves for all physical field variables. That is, within each Newton iteration step, the four unknown fields—fluid velocity, pressure, solid displacement, and temperature—are updated synchronously within a unified variational framework, rather than being solved independently and then loosely coupled iterated. This block-based solution strategy fundamentally avoids the additional mass effects and slow convergence numerical stability problems that may occur when using partitioned solution methods to handle fluid-structure systems with similar densities.

[0026] Step S3: Determine whether the residual of the current iteration has converged. If it has not converged, adjust the Newton iteration step size and return to step S2 to continue iterating. If it has converged, update the mesh and field quantity according to the converged solution, adaptively determine the step size of the next time step, perform velocity extrapolation and enter the next time step loop until the simulation ends, and output the dynamic terrain evolution results.

[0027] After completing the multiphysics solution within a time step, the system enters the residual judgment and iteration control phase. The system determines whether the residual achieved by the current nonlinear iteration is lower than the preset convergence tolerance. If the residual has not yet converged, the Newton iteration step size or correction line search parameters are adjusted, and the system returns to the coupled solution step to continue iteration until the residual meets the accuracy requirements. If the residual has converged, it indicates that the nonlinear coupled solution for the current time step has been completed. The system then updates the fluid-structure interface position and fluid domain mesh configuration based on the converged displacement solution, and calculates the mesh motion velocity satisfying the biharmonic equations based on the new mesh position to ensure that the internal mesh nodes adjust smoothly and without distortion with boundary deformation. Simultaneously, based on the current flow field, temperature field, and mesh velocity, and combined with numerical stability conditions, the step size for the next time step is adaptively determined, and the velocity field is extrapolated to provide initial value predictions for the next time step. Finally, all field quantities and mesh information are updated, and the calculation is advanced to the next time step cycle. This adaptive step size control mechanism enables the method to maintain the stability and efficiency of numerical integration throughout the simulation process, spanning multiple time orders from earthquake cycles to orogenic movements.

[0028] Through the above three steps, this method unifies the equations of fluid mechanics, solid mechanics, and thermodynamics within an arbitrary Lagrange-Euler variational framework. By employing a single-block nonlinear solution strategy, it achieves closed-loop coupled simulation of heat-driven, force-transfer, and deformation feedback. This unified numerical model can fully characterize the core geodynamic process of mantle convection-driven lithospheric deformation and the feedback of deformation on heat transport, providing a high-fidelity numerical simulation tool for geological science problems such as plate tectonics simulation, orogenic belt evolution, and earthquake periodicity research.

[0029] like Figure 2 The diagram illustrates the multiphysics logic architecture for fluid-solid-thermal coupling solution of this invention. This method divides the computational domain into a fluid domain and a solid domain, applying the Stokes equations with the Boussinesq approximation and the solid equilibrium equations described by the Piola-Kirchhoff stress, respectively. At the interface, by forcibly satisfying the conditions of consistent interface traction force and consistent interface displacement, the two physical fields are assembled into a unified single-system matrix. This matrix is ​​then efficiently solved through algebraic multigrid preprocessing, ultimately outputting the dynamic terrain calculation results. Figure 3The “initial boundary conditions” in this invention refer to the initial field, extrapolated initial values, and various spatial boundary / interface constraints required for the numerical simulation of this invention. Among them, the most critical are the Moho surface fluid-structure continuity, normal stress balance, and biharmonic mesh motion boundary conditions.

[0030] In step S2, a block Schur complement decomposition-based acceleration strategy is adopted to jointly solve all physical field variables. The core idea of ​​this strategy is to utilize the naturally existing block sparse structure in the algebraic equations of multi-physics coupled systems. By grouping variables and reducing matrix order, the solution of large-scale coupled linear systems is transformed into a combination of a significantly reduced system solution and a low-cost matrix back substitution operation, thereby reducing the computational overhead of the linear solution stage in a single Newton iteration.

[0031] Specifically, step S2 involves jointly solving for all physical field variables, including: S21: Identify the block sparse structure of algebraic equations in nonlinear block-coupled iterations, and divide the unknowns into the first variable group and the second variable group according to the physical field type.

[0032] This study identifies the sparse block structure of algebraic equations in nonlinear block-coupled iterations and divides all unknowns into a first variable group and a second variable group according to the type of physical field. In geodynamic multiphysics simulations, each Newton iteration step requires solving a linearized algebraic equation system containing four types of unknowns: fluid velocity, pressure, solid displacement, and temperature. The coefficient matrix of this equation system is not densely packed but exhibits a sparse characteristic of natural block partitioning according to physical fields: apart from the diagonal sub-blocks of each physical field itself, only the cross sub-blocks related to fluid-solid interfaces and thermodynamic coupling have non-zero elements, while some physical fields, such as pressure and solid displacement, and pressure and temperature, have no direct coupling relationship, and the corresponding positions are zero sub-blocks. Step S21 utilizes this structural feature to reorganize the four types of unknowns according to their physical roles in the coupled system: fluid-solid interface displacement and temperature are classified into the first variable group, which are the key driving forces of multiphysics coupling, carrying the kinematic information and thermal buoyancy driving information of the fluid-solid interface; fluid velocity, pressure, and other solid displacements other than interface displacement are classified into the second variable group, which can be derived and solved by the first variable group through physical field equations. This division method provides a clear basis for the subsequent construction of the reduced system.

[0033] S22: Construct the Schur complement matrix for the partitioned variable groups to form a reduced system with the first variable group as the unknown.

[0034] The construction process is as follows: First, an approximate inverse is constructed for the submatrix block corresponding to the second variable group. Then, using this approximate inverse, all unknowns of the second variable group are eliminated from the original system through matrix row and column operations, resulting in a reduced linear system with only the first variable group as unknowns. The size of this reduced system depends only on the sum of the degrees of freedom at the fluid-structure interface nodes and the degrees of freedom at the global temperature nodes, which is much smaller than the original algebraic equations that simultaneously include all velocity, pressure, and displacement degrees of freedom. In typical three-dimensional simulation scenarios of geodynamics, if a million-level grid element is used for discretization, the number of unknowns in the original equations can reach tens of millions of degrees of freedom. However, the size of the system after Schur complement reduction can usually be reduced to one-tenth or even less of the original size, and the condition number of the reduced system is usually better than that of the original full system, which provides more favorable numerical conditions for subsequent linear solutions.

[0035] S23: Solve the reduced system to obtain the correction amount for the first variable group.

[0036] The reduced system is a linear system of equations with only fluid-solid interface displacement and temperature as unknowns, and its coefficient matrix is ​​the Schur complement matrix. Compared with the original full system of algebraic equations, the reduced system has two significant advantages. In terms of scale, the unknowns of the original system cover the velocities, pressures, displacements, and global temperatures of all nodes in both the fluid and solid domains, typically reaching millions or even higher degrees of freedom. The reduced system, however, retains only the displacement degrees of freedom related to the fluid-solid interface nodes and the global temperature degrees of freedom, significantly reducing the number of unknowns. In practical geological models, such as lithosphere-mantle coupling models spanning hundreds of kilometers and containing detailed fault geometry, the scale of the reduced system can be reduced to a fraction of the original system. In terms of numerical properties, the condition number of the Schur complement matrix is ​​usually superior to that of the coefficient matrix of the original full system. This is because eliminating the second variable group is essentially equivalent to performing an accurate block Gaussian elimination on the original matrix, resulting in a Schur complement matrix with better spectral properties, which is beneficial for the rapid convergence of the iterative solver. Based on the above two advantages, the solution of this reduced system can be completed by choosing an efficient sparse direct solver or by using the preprocessed Krylov subspace iteration method. Regardless of the method used, the computational resources and computation time occupied by the solution process are significantly reduced compared to directly solving the original whole system, and the loss of solution accuracy is within the order of magnitude of numerical rounding error, which will not have an adverse effect on the overall convergence of Newton's iteration.

[0037] S24: Substitute the correction amount of the first variable group back to calculate the correction amount of the second variable group.

[0038] After obtaining the corrections for the fluid-solid interface displacement and temperature, these known quantities are substituted into the second variable group subsystem of the original equation system. The corrections for fluid velocity, pressure, and other solid displacements are then directly calculated using a single matrix-vector multiplication and a single vector subtraction. This back-substitution process is essentially an explicit calculation, not involving the solution of any linear equation system, and its computational cost is far lower than that of iterative solutions. The reason the calculation of the second variable group corrections can be completed at such a low cost is that the variable grouping strategy in step S21 and the Schur elimination process in step S22 have fully encoded the coupling relationship between the two variable groups into matrix operations. The back-substitution step is simply a direct execution of this algebraic relationship. Although the variables in the second variable group are physically driven by the first variable group, their corrections are equally necessary to ensure the convergence accuracy of the entire nonlinear system. In geodynamic simulations, the velocity field of mantle convection determines heat transport efficiency, and the stress distribution within the lithosphere determines fault activity. Accurate solutions for these physical quantities depend on a complete update of all unknowns in each Newton iteration step, rather than just updating the critical subset of interface displacement and temperature.

[0039] S25: Update all physical field variables together with the correction values ​​of the first variable group and the second variable group to complete this joint solution.

[0040] The correction values ​​for the first set of variables obtained in step S23 and the correction values ​​for the second set of variables obtained in step S24 are simultaneously applied to the current iteration value. The fluid velocity field, pressure field, solid displacement field, and temperature field are updated at the same time, forming new iteration values ​​for residual judgment or for proceeding to the next Newton step. This joint update method ensures the synchronous advancement of the four types of physical field variables in each Newton iteration step, avoiding the interface constraint mismatch and additional mass effect introduced by the asynchronous updates of fluid and solid in the partitioned solution method, and also avoiding the time lag error between temperature and deformation when the thermodynamic field and the mechanical field are solved separately. After a single joint solution is completed, the updated physical field is used as the initial value for the next Newton iteration step to continue participating in the nonlinear iteration until the residual of the current time step converges to the preset tolerance. The above-described solution strategy, through in-depth utilization of the algebraic structure of multiphysics coupled systems, decomposes the problem of solving large-scale coupled linear systems into a combination of a significantly reduced-scale system solution and a low-cost matrix back substitution operation by using a Schur complement matrix reduction operation. While maintaining the accuracy and convergence robustness of the single-block coupled solution, it effectively reduces the computational cost of a single Newton iteration, providing key support for the engineering feasibility of large-scale three-dimensional dynamic terrain thermodynamic simulation under limited computational resources.

[0041] Furthermore, the identification of the block sparse structure of the algebraic equation system in the nonlinear block coupling iteration in step S21 is completed through two sequentially executed operations. This identification process is a prerequisite for the entire block Schur complement acceleration strategy, and its purpose is to extract structural information from the algebraic equation system formed after the multiphysics coupling system is discretized, which can be used for subsequent variable grouping and matrix order reduction.

[0042] Specifically, step S21, identifying the block sparse structure of the algebraic equation system in the nonlinear block-coupled iteration, includes: S211: The algebraic equations are divided into blocks according to the four physical field variables of fluid velocity, pressure, solid displacement and temperature to obtain a block sparse matrix.

[0043] In each Newton iteration step of the multiphysics simulation of geodynamics, the algebraic equations formed after finite element discretization and linearization contain all unknowns at all nodes. If arranged directly according to the node numbering order, the distribution of each physical field variable in the matrix will present an intertwined and mixed form, making it difficult to intuitively reveal the coupling relationship between different physical fields and the structural characteristics of the matrix. The block arrangement operation performed in step S211 groups the unknowns belonging to the same physical field type in the equation system together, forming four large variable blocks in the order of fluid velocity, pressure, solid displacement, and temperature. Correspondingly, the coefficient matrix is ​​also divided into sixteen sub-matrix blocks of four rows and four columns. This arrangement makes the algebraic coupling relationships between different physical fields clearly discernible: the fluid velocity sub-matrix block originates from the momentum conservation part of the Navier-Stokes equations, the pressure sub-matrix block from the mass conservation constraint, the solid displacement sub-matrix block from the discretization of the finite-deformation solid mechanics equations, the temperature sub-matrix block from the discretization of the heat conduction-convection equations, and the cross sub-matrix blocks at off-diagonal positions correspond to the coupling effects between different physical fields—for example, the cross block between fluid velocity and temperature corresponds to the thermal buoyancy term, and the cross block between fluid velocity and solid displacement corresponds to the velocity continuity constraint at the fluid-solid interface. The distribution pattern of the non-zero sub-blocks in the block sparse matrix obtained after the block arrangement directly reflects the existence and strength of the coupling relationships between the physical fields in the original physical problem.

[0044] Discretization yields a nonlinear algebraic system:

[0045] Where F is the discretized global residual function, which characterizes the unbalance of each physical field equation (fluid momentum / mass conservation, solid deformation, thermodynamics) after finite element discretization; This is the vector of unknowns at the (n+1)th time step, containing all nodal degrees of freedom; It is the set of unknown fluid velocities; For the set of unknown pressure quantities; Let T be the set of unknown quantities for solid displacement; T is the set of unknown quantities for temperature; the superscript T is transpose.

[0046] in Including all discrete residual terms, the system has a block sparse structure: .

[0047] The above formulas represent a system of block sparse algebraic equations for a fluid-solid-thermal coupled system, obtained through linearization during Newton's iterative solution process. The block matrix on the left is the system's Jacobian matrix, and the block subscripts are... v , p , u , T These represent the four physical fields: velocity, pressure, displacement, and temperature. A ij Indicates the first i The physical field equations for the first j A block of partial derivative submatrices of physical field variables (e.g., A vv It is the diagonal block corresponding to the fluid velocity, which comes from the time derivative term and viscosity term in the fluid momentum equation; A vp The velocity-pressure coupled block is derived from the pressure gradient term and divergence constraint. A vT The velocity-temperature coupling block is derived from the thermal buoyancy term in the Boussinesq approximation; A uu The diagonal block of solid displacement is derived from the stress term in the finite deformation equation of solids; A TT The temperature diagonal block is derived from the heat conduction and convection equations; A uv It is a displacement-velocity coupling block, which originates from the drag force or stress balance at the fluid-structure interface; A uT It is a displacement-temperature coupling block, originating from thermal stress or thermal expansion effects; A Tv It is a temperature-velocity coupled block, originating from the thermal convection term; A Tu This is a temperature-displacement coupled block, originating from the convection correction caused by mesh motion under the ALE framework; meanwhile, the zero blocks (0) in the block matrix indicate that the corresponding physical fields do not have a direct algebraic coupling relationship in the current equation form, such as the pressure equation which does not explicitly include displacement and temperature variables); the Δ in the left column vector V ,Δ P ,Δ U ,Δ TThese represent the increments of the fluid velocity, pressure, solid displacement, and temperature to be solved (i.e., the corrections in Newton's iteration), with the symbol Δ indicating the increment; the column vectors on the right... Rv , Rp , Ru , RT These represent the residual vectors of the fluid momentum conservation equation, fluid mass conservation equation, solid finite deformation equation, and global thermodynamic equation in the current iteration state, respectively. The negative sign on the rightmost side of the equation indicates that Newton's method is solving the incremental equation. J Δ X =- R That is, by solving these increments, the current variable values ​​are corrected so that the total residual of each physical field equation approaches zero.

[0048] S212: Identify the symmetric positive definiteness of the fluid velocity submatrix and the block diagonality of the solid displacement submatrix in a block sparse matrix. The symmetric positive definiteness is used to guide the selection of a geometric multigrid preprocessor, and the block diagonality is used to guide the selection of an algebraic multigrid preprocessor.

[0049] Identifying the numerical properties of each submatrix block within the sparse block matrix is ​​crucial for determining the selection of subsequent preprocessors. The fluid velocity submatrix originates from the discretization of the viscous terms in the Stokes or Navier-Stokes equations. After applying the standard Galerkin finite element method and appropriate pressure stabilization, this submatrix exhibits symmetric positive definiteness. This property means that a geometric multigrid preprocessor can be used for the fluid velocity-related submatrix during subsequent solution processing. The geometric multigrid method utilizes smooth iteration between nested mesh levels and coarse mesh correction, offering optimal or near-optimal computational complexity for symmetric positive definite elliptic problems, achieving iterative convergence at a cost that is linearly proportional to the number of unknowns. The solid displacement submatrix originates from the discretization of finite deformation elastic or viscoelastic-plastic equations. In large-scale three-dimensional discretized models, due to the coupling between displacement components in different spatial directions and the stiffness differences between different material regions, this submatrix exhibits a block diagonal structure. That is, the displacement degrees of freedom of each node form dense small blocks locally, while the stiffness coupling between different nodes forms sparse off-diagonal blocks. The identification of this property implies that the structure of the solid displacement submatrix is ​​suitable for an algebraic multigrid preprocessor. The algebraic multigrid method does not rely on geometric grid hierarchy information but directly constructs coarsening levels and interpolation operators based on the matrix's algebraic connectivity and coefficient values. It exhibits good preconditioning effects for matrices with block diagonal structures derived from elliptic partial differential equation discretization. The identification results of these two submatrix properties directly provide a basis for the design of subsequent hybrid multigrid preprocessing strategies. Within a unified linear solution framework, the most suitable preprocessing method for the numerical characteristics of the submatrix blocks corresponding to different physical fields is applied, thereby maximizing the overall solution efficiency.

[0050] Furthermore, in step S21, the method of dividing according to the physical field type is as follows: fluid-solid interface displacement and temperature are used as the first variable group, and fluid velocity, pressure and solid displacement other than interface displacement are used as the second variable group. This allows the key physical quantities controlling fluid-solid coupling and thermo-mechanical coupling to be assigned to the first variable group, and the subordinate physical quantities that can be driven by the first variable group through coupling relationship to be assigned to the second variable group.

[0051] The method of dividing variables into groups based on physical field type is not a simple mechanical splitting of four categories of unknowns, but rather a functional grouping based on the dynamic role of each physical quantity in the multiphysics coupling system. Under this division, the first variable group includes fluid-solid interface displacement and temperature, while the second variable group includes fluid velocity, pressure, and other solid displacements besides interface displacement. In this grouping strategy, fluid-solid interface displacement directly determines the geometric configuration and motion state of the interface between mantle convection and lithosphere deformation, serving as the geometric carrier for the mechanical coupling between the fluid and solid domains. Temperature, on the other hand, is the core thermodynamic variable driving the evolution of the entire system. The temperature gradient within the mantle, formed by radioactive decay and residual heat, generates buoyancy through thermal expansion, driving mantle convection. The convection pattern, in turn, determines the efficiency of heat transport to the bottom of the lithosphere, thus affecting the thermal structure and rheological intensity of the lithosphere. These two types of physical quantities together constitute the key driving information controlling fluid-solid coupling and thermodynamic coupling. The updated quantities are solved preferentially in each Newton iteration, capturing the core nonlinear changes in multiphysics coupling.

[0052] Furthermore, the core task of step S22 is to transform the original large-scale algebraic equation system into a significantly reduced system after the variable group partitioning is completed, using the block elimination technique of linear algebra. This is achieved through two sequentially executed operations: first, constructing the approximate inverse of the submatrix of the second variable group, and then using this approximate inverse to eliminate all the unknowns of the second variable group.

[0053] Specifically, in step S22, a Schur complement matrix is ​​constructed for the partitioned variable groups to form a reduced system with the first variable group as the unknown, including: S221: Construct an approximate inverse for the submatrix of the second variable group.

[0054] In the implementation of the block Schur complement method, constructing the approximate inverse of the second variable set submatrix is ​​the key operation of the entire order reduction strategy. The second variable set includes fluid velocity, pressure, and solid displacements other than interface displacement. These unknowns constitute the largest part of the original algebraic equations, and the corresponding submatrix is ​​a large-scale sparse matrix block. Solving for its exact inverse matrix would be computationally expensive, equivalent to directly solving the original entire system, completely negating the efficiency gains of the order reduction strategy. Therefore, an approximate inverse is constructed by using appropriate preprocessing techniques to generate an approximate inverse operator for this submatrix. This approximate inverse operator maintains sufficient accuracy to ensure the accuracy of the Schur complement matrix calculation while being obtained at a much lower computational cost than exact inversion. In practice, this approximate inverse can be obtained by applying several multi-grid iterations or Krylov subspace iterations to the second variable set submatrix. The core principle is to exchange significant computational savings for controllable accuracy loss.

[0055] S222: By using approximate inverses, the second set of variables is eliminated through matrix operations, resulting in a reduced linear system with only the first set of variables as unknowns. The scale of the reduced system is smaller than the original algebraic equations, and the solution of the reduced system is directly the correction amount of the first set of variables.

[0056] The specific elimination process is as follows: by using the product of the approximate inverse operator and the cross-coupled submatrix, the unknowns of the second variable group are completely eliminated from the original equation system through matrix subtraction and row and column operations, forming the Schur complement matrix containing only the unknowns of the first variable group and the corresponding right-hand side terms.

[0057] Linear solution acceleration based on block Schur complement decomposition: First, solve the reduced system,

[0058] in, S The Schur complement matrix represents the reduced system coefficient matrix after eliminating the second set of variables, retaining only the first set (interface displacement, temperature); Δ Y This is the correction amount for the first set of variables (i.e., the increment of interface displacement and temperature in each Newton iteration). R Y For the residuals (right-hand side) corresponding to the first variable group, Ry is also the Schur complement of the original residuals; Y The first set of variables includes fluid-solid interface displacement and temperature; X The second variable group includes fluid velocity, pressure, and solid displacement excluding interface displacement; A YY This is the submatrix block corresponding to the first variable group itself; A YX This is the coupling submatrix between the first variable group and the second variable group; A XY This is the coupling submatrix between the second variable group and the first variable group; A XX This is the submatrix block corresponding to the second variable group itself; the superscript -1 indicates matrix inversion.

[0059] Utilizing the characteristics of geophysical problems: matrices Symmetric positive definite Diagonal blocks can be solved in parallel; The preprocessor is designed for: geometric multigrids for Algebraic multigrids are used for .

[0060] This elimination process is mathematically a block Gaussian elimination, compressing the original large-scale linear system with four coupled physical field variables into a small-scale linear system containing only fluid-structure interface displacement and temperature degrees of freedom. The size of the reduced system can typically be reduced to a fraction of the original system. More importantly, the solution of this reduced system, i.e., the vector output by the linear solver, is directly the correction value for the first variable group, without any further post-processing or variable transformation steps. This property allows the solution result of the reduced system to be directly used in step S23, providing a clear data interface for subsequent back-substitution calculations.

[0061] Furthermore, in step S221, the process of constructing the approximate inverse of the second variable group submatrix employs a hybrid multigrid preprocessing strategy. This strategy addresses the fundamental differences in numerical properties and matrix structure between different physics field submatrices in a multiphysics coupled system, selecting the most suitable preprocessing technique within the same solution framework, thereby achieving a better approximate inverse construction effect overall than a single preprocessing method.

[0062] Specifically, in step S221, the submatrix of the second variable group is constructed using a hybrid multigrid preprocessing method to approximate its inverse: S221a: Apply geometric multigrid preprocessing to the sub-matrix block corresponding to fluid velocity and pressure to obtain the approximate inverse of the sub-matrix block.

[0063] The fluid velocity and pressure sub-matrix blocks are derived from the finite element discretization of the incompressible Navier-Stokes equations or Stokes equations under the Boussinesq approximation. When spatial discretization is performed using a hybrid finite element scheme, this sub-matrix block exhibits a saddle point structure, where the velocity-velocity diagonal blocks correspond to the elliptic operator discretization of the viscous terms, possessing symmetric positive definiteness; the pressure-related blocks correspond to the divergence operator discretization of the incompressible constraints. Geometric multigrid methods have a natural advantage for such matrices derived from elliptic partial differential equation discretization. This method explicitly constructs a nested mesh hierarchy from fine to coarse meshes, performs smooth iterations on the fine mesh to eliminate high-frequency error components, and performs residual correction on the coarse mesh to eliminate low-frequency error components. Through recursive loops between multi-level meshes, it achieves computational complexity approximately linear with the number of unknowns. In the implementation of actual geodynamic models, the mesh hierarchy of geometric multigrids can directly inherit the mesh refinement hierarchy information generated during finite element discretization without additional construction. The coarsening operator adopts natural injection or piecewise constant constraint, and the interpolation operator adopts the extension method matching the finite element shape function. This makes the application of geometric multigrid preprocessing on this sub-matrix block both efficient and easy to implement.

[0064] S221b: Apply algebraic multigrid preprocessing to the sub-matrix block corresponding to the solid displacement to obtain the approximate inverse of the sub-matrix block.

[0065] Solid displacement submatrix blocks originate from the discretization of the finite-deformation elastic or viscoelastic-plastic solid mechanics equations. Unlike fluid velocity submatrices, solid displacement submatrices, in different material regions within the lithosphere (e.g., the brittle upper crust, ductile lower crust, and mantle lithosphere), exhibit significant differences in Lamé constants and yield strengths, resulting in stiffness matrix element values ​​that can span orders of magnitude. Furthermore, the local refinement and complex geometry of fault regions mean that solid domain meshes often lack the regular nested hierarchical structure of fluid domain meshes. Algebraic multigrid methods do not rely on explicit geometric mesh hierarchy information but instead directly determine coarsening strategies based on the algebraic connectivity and coefficient values ​​of the matrix, constructing virtual coarse mesh hierarchies, interpolation operators, and coarsening operators by analyzing the size and connectivity of matrix elements. This method naturally adapts to drastic jumps in matrix coefficients and effectively handles variations in stiffness matrix numerical properties caused by unstructured meshes and inhomogeneous material parameters, making it particularly suitable for the approximate inverse construction of solid displacement submatrix blocks.

[0066] S221c: Assemble the approximate inverses obtained in the two steps according to the block diagonal structure to form the overall approximate inverse of the second variable group submatrix.

[0067] Since the fluid velocity and pressure sub-matrix block and the solid displacement sub-matrix block correspond to different physical fields and different spatial regions, the coupling between them only occurs at the fluid-solid interface nodes. There is no direct algebraic coupling within the internal sub-matrixes of the second variable group. Therefore, the approximate inverses of the two sub-matrix blocks can be combined in a block-by-block diagonal manner to form the approximate inverse of the overall sub-matrix of the second variable group. This assembly method is mathematically equivalent to assuming decoupling between the two sub-blocks at the approximate inverse level, maintaining the independence and specificity of the preprocessing strategies for each sub-block while simplifying the construction logic of the overall approximate inverse. Through the aforementioned hybrid multi-grid preprocessing strategy, the approximate inverse of the second variable group sub-matrix can fully utilize the elliptical structure advantage and geometric grid hierarchy information of the fluid sub-blocks, while also adapting to the non-homogeneity and unstructured characteristics that may occur in the solid sub-blocks in complex geological structures. This provides a fundamental operator with both accuracy and computational efficiency for the subsequent Schur elimination process.

[0068] Furthermore, the Boussinesq approximation of the incompressible Navier-Stokes equations includes the ALE time derivative term, convection term, pressure gradient term, viscous diffusion term, and thermal buoyancy term, and satisfies the velocity divergence condition. The Boussinesq approximation refers to an approximation method in which density change only acts on the thermal buoyancy term, while the density in the other terms is regarded as a constant. The thermal buoyancy term is determined by the coefficient of thermal expansion, the difference between the temperature and the reference temperature, and the acceleration due to gravity.

[0069] The ALE time derivative term describes the rate of change of fluid velocity under a reference configuration when the grid moves in any way. The velocity in the convection term needs to be subtracted from the grid velocity to correct for the material transport effect caused by grid motion. The pressure gradient term and the viscous diffusion term together constitute the stress divergence inside the fluid. Pressure, as a Lagrange multiplier, forces the fluid to satisfy the incompressibility constraint, while the viscous term characterizes the momentum diffusion of mantle material during creep deformation. The thermal buoyancy term is determined by the coefficient of thermal expansion, the difference between the local temperature and the reference temperature, and the gravitational acceleration vector. It is the only volumetric force source driving mantle convection. When high-temperature, low-density mantle material is in a gravitational field, it generates an upward buoyancy force; conversely, low-temperature, high-density material generates a downward settling force. The core of the Boussinesq approximation is that density change only acts on this buoyancy term, while the density in the other terms of the momentum equation and the continuity equation is taken as the reference density. This treatment preserves the core physical mechanism of thermal drive in mantle convection simulation while significantly simplifying the mathematical structure of the equations.

[0070] In the finite deformation equations of the solid domain, the first Piola-Kirchhoff stress tensor is related to the second Piola-Kirchhoff stress tensor through the deformation gradient. The second Piola-Kirchhoff stress tensor adopts the Saint-Venant Kirchhoff hyperelastic model expressed in terms of Lamé constant and Green-Lagrange strain tensor. The Green-Lagrange strain tensor is composed of the gradient of the solid displacement and the product of its transpose and the gradient.

[0071] In finite deformation theory, the choice of stress metric directly affects the formulation of constitutive equations and the ease of numerical implementation. The first Piola-Kirchhoff stress tensor is a two-point tensor that maps a unit-area surface element under the reference configuration to the actual force vector under the current configuration. Its mechanical meaning is clear and intuitive, facilitating the establishment of momentum balance equations in the reference configuration and their direct correlation with external force conditions at the boundaries. However, the first Piola-Kirchhoff stress tensor itself is not a symmetric tensor and is unsuitable as a direct stress metric for the constitutive model. The constitutive equation requires that the stress metric and strain metric satisfy the principle of objectivity within the material framework, i.e., the constitutive response is unaffected by rigid body rotation. The second Piola-Kirchhoff stress tensor is introduced precisely to meet this requirement. It is fully defined in the reference configuration and, through the inverse transformation of the deformation gradient, pulls the actual stress in the current configuration back to the reference configuration and subtracts the influence of rigid body rotation. It possesses symmetry and can form a conjugate pair with the Green-Lagrange strain tensor, which is also defined in the reference configuration. The first Piola-Kirchhoff stress tensor and the second Piola-Kirchhoff stress tensor are related by the deformation gradient. The former is equal to the tensor product of the deformation gradient and the latter. This relationship provides a bridge for simultaneously handling the equilibrium equations and constitutive equations under the reference configuration.

[0072] The second Piola-Kirchhoff stress tensor is expressed using the Saint-Venant Kirchhoff hyperelastic model, which uses Lamé constants and Green-Lagrange strain tensors. The Saint-Venant Kirchhoff model is the most fundamental hyperelastic model in finite deformation elasticity theory. It assumes a linear relationship between the second Piola-Kirchhoff stress and the Green-Lagrange strain, with the proportionality constant determined by the two Lamé constants. This model is formally consistent with the generalized Hooke's law in small deformation linear elasticity theory. However, since the Green-Lagrange strain tensor itself contains a quadratic term of the displacement gradient, this model can actually describe the elastic response of materials under large displacement and rotation conditions. It is applicable to geological deformation scenarios where the lithosphere undergoes significant displacement under tectonic stress, but the strain remains at a moderate level. Hyperelasticity implies the existence of an elastic potential energy function in this model. The second Piola-Kirchhoff stress can be obtained by taking the partial derivative of this potential energy with respect to the Green-Lagrange strain. This ensures the thermodynamic consistency of the material constitutive relation, i.e., no dissipation occurs in a closed elastic deformation cycle. The Green-Lagrange strain tensor is composed of the gradient of the solid displacement and the product of its transpose and the gradient. Its definition includes a first-order term and a second-order term of the displacement gradient. The first-order term is the symmetric displacement gradient, which describes infinitesimal strain. The second-order term is the product of the displacement gradient transpose and the displacement gradient, which provides geometric nonlinear correction under large rotation or large displacement conditions, ensuring that the strain degree remains zero under rigid body rotation, thus satisfying the objectivity requirement.

[0073] Furthermore, in step S2, the momentum and mass conservation equations for the fluid domain adopt the incompressible Navier-Stokes equations approximated by Boussinesq. In the context of geodynamics, mantle convection is a low-speed thermal convection process driven by a temperature gradient. Its characteristic velocity is typically only on the order of a few centimeters per year, with Mach numbers far less than 1. Density changes are mainly caused by thermal expansion, and the relative change is usually no more than 10% to 20%. The Boussinesq approximation is a simplified model established for buoyancy-driven flow problems where density changes are far less than the average density. Its core approach is to use the reference density in all terms of the momentum equation except the buoyancy term, retaining only the linear change in density with temperature in the buoyancy term. This approach simplifies the mantle convection problem numerically into an incompressible flow problem with constant density, greatly reducing the nonlinear coupling of the equations and the difficulty of solving them, while fully preserving the core physical mechanism of thermal buoyancy-driven flow—that is, temperature anomalies generate volume forces through the coefficient of thermal expansion and gravitational acceleration, driving the mantle material to rise or sink. The incompressible constraint requires zero divergence in the velocity field, which physically corresponds to the observation that the density of mantle material changes very little over geological timescales. It also means that pressure acts as a Lagrange multiplier in the governing equations, enforced to satisfy the mass conservation condition rather than being independently given by the equations of state. The expression of this system of equations within an arbitrary Lagrange-Eulerian framework requires the introduction of grid velocities to correct the convection terms, ensuring that fluid motion can still be accurately described under the reference configuration even as the grid moves with the solid boundary.

[0074] Specifically, the governing equations for the fluid-structure interaction system are as follows: Fluid domain (mantle / lower crust): Incompressible Navier-Stokes equations using the Boussinesq approximation.

[0075] in ρ f This refers to the fluid density (density of mantle / lower crust material). v f For fluid velocity vector; t For time; The time derivative of ALE (Arbitrary Lagrange-Euler) represents the rate of change over time under the reference configuration; w This is the ALE mesh velocity vector, which represents the velocity of the mesh as it deforms along the solid boundary. For spatial gradient operators; p For the pressure field in the fluid domain; η For fluid dynamic viscosity; For the Laplace operator; α The coefficient of thermal expansion; T This represents the current temperature field.T 0 represents the reference temperature (initial temperature); g This is the gravitational acceleration vector.

[0076] The finite deformation equations for the solid domain employ hyperelastic or viscoelastic-plastic constitutive equations described by the Piola-Kirchhoff stress tensor. The deformation of the lithosphere under tectonic stress is often not considered a small deformation problem, especially in orogenic belts, subduction zones, and fault-prone areas, where the displacement of solid materials can reach tens of kilometers or even greater, far exceeding the applicability of the small deformation assumption of linear elasticity theory. Therefore, the solid domain adopts finite deformation equations, using the deformation gradient tensor as the basic kinematic measure. Rigid body rotation and pure deformation are separated through pole decomposition, and the magnitude of the finite deformation is measured using the Green-Lagrange strain tensor.

[0077] Solid domain (lithospheric / upper crust): using finite deformation elastic or viscoelastic-plastic equations

[0078] in, For divergence operators; P For the first Piola-Kirchhoff stress tensor (a two-point tensor that maps reference configuration elements to the actual forces of the current configuration); ρ s Density of the solid domain (lithospheric / upper crust); b It is a vector of body forces (such as gravity, mantle buoyancy).

[0079] For stress description, the Piola-Kirchhoff stress tensor system is adopted. The first Piola-Kirchhoff stress tensor maps surface elements under the reference configuration to the actual stress vector under the current configuration. It is a two-point tensor, suitable for establishing mechanical equilibrium between the reference and current configurations, facilitating stress equilibrium matching with the fluid equations expressed under the reference configuration in the ALE framework at the interface. The hyperelastic constitutive model assumes that stress is directly given by the partial derivative of the elastic potential function with respect to Green-Lagrange strain, and that deformation is thermodynamically reversible. It is suitable for describing the recoverable elastic deformation of the lithosphere before reaching the yield limit. The viscoelastic-plastic constitutive model further introduces time-dependent relaxation mechanisms and stress-dependent yielding and plastic flow mechanisms, enabling it to describe creep deformation of rocks under long-term geological stress and irreversible plastic strain when stress exceeds the yield surface. The choice between the two constitutive models depends on the temperature, pressure, and rock type of the simulated geological unit. Different regions can be specified within the same computational domain to achieve a continuous transition in material behavior from the brittle upper crust to the ductile lower crust.

[0080] Furthermore, the mesh motion in updating the mesh and field quantities in step S3 is controlled by a biharmonic equation, which is a mesh motion equation with the composition of the Laplacian operator as the control operator. Within any Lagrange-Euler framework, the motion of mesh nodes within the fluid domain needs to be adjusted according to the deformation of the fluid-solid interface. However, the mesh motion mode itself is not directly determined by the physical equations; auxiliary equations are needed to construct a mesh displacement field that satisfies the boundary conditions and has a reasonable distribution. The biharmonic equation uses the composition of the Laplacian operator as the control operator.

[0081] Specifically, the ALE mapping is constructed as follows: Fluid mesh motion is constructed using biharmonic equations:

[0082] Boundary conditions: , .

[0083] Where, Δ 2 It is a biharmonic operator, i.e., the square of the Laplace operator ( ); For fluid domain mesh displacement; Fluid domain under reference configuration; For solid displacement (as a forced boundary condition for mesh displacement at the interface); For fluid-solid interface (Mohorovičić discontinuity); This represents the boundary of the fluid domain.

[0084] This choice has important numerical properties: compared to simple Laplace smoothing, the solution of the biharmonic equation has higher smoothness, and its second derivative is continuous, effectively avoiding mesh line folding or excessive distortion in large deformation regions. Simultaneously, the biharmonic operator propagates boundary displacement information into the fluid domain in a more uniform manner, unlike the Laplace equation which concentrates boundary disturbances mainly in local regions near the boundary, thus producing a more uniform and coordinated mesh motion distribution throughout the entire fluid domain. In practical implementation, the biharmonic equation is usually transformed into a hybrid form, introducing the Laplace of mesh displacement as an auxiliary variable, reducing the fourth-order equation to two coupled second-order equations for finite element solution. This allows the use of the existing finite element discretization framework while reducing the continuity requirements of higher-order derivatives on shape functions. By controlling the mesh movement using biharmonic equations, the entire simulation can maintain the quality of the fluid domain mesh and the stability of the numerical solution even under significant tectonic deformations in the lithosphere, such as the significant depression of the Moho discontinuity caused by crustal shortening and thickening in orogenic belts or the geometrical abrupt change of the interface caused by fault displacement. This provides technical support at the mesh geometry level for long-term simulation of large-deformation dynamic terrain evolution.

[0085] The above heat conduction equation: ； in, ρ For material density (values ​​are taken according to the material properties of fluids or solids); Specific heat capacity (constant pressure); T For temperature field; t For time; v For velocity field (fluid velocity or solid material velocity, depending on the domain); For gradient operators; k This is the thermal diffusivity (thermal conductivity). Q It is an internal heat source (including source terms such as shear heat generation and radioactive decay heat).

[0086] Furthermore, the momentum and mass conservation equations in the fluid domain, the finite deformation equations in the solid domain, and the global thermodynamic equations are constructed in a unified variational weak form under the reference configuration. The unified variational weak form takes fluid velocity, pressure, solid displacement, and temperature as unknowns and includes fluid momentum terms, fluid convection terms, fluid stress terms, pressure divergence terms, solid inertia terms, solid stress terms, thermal transient and convection terms, thermal diffusion terms, as well as buoyancy terms, volume force terms, and internal heat source terms. The velocity continuity and stress equilibrium conditions at the fluid-solid interface are satisfied in the unified variational weak form.

[0087] In the finite element method, the strong form of partial differential equations requires that the unknown functions satisfy differential relations everywhere in the solution domain. The variational weak form relaxes this requirement to satisfy it in the sense of weighted integration. That is, both sides of the equation are multiplied by any test function satisfying homogeneous boundary conditions and integrated over the entire domain. By integrating by parts, the continuity requirement of higher-order derivatives is transferred to the test function, thereby reducing the smoothness requirement for the approximate function and making numerical discretization mathematically feasible and computationally flexible. This invention unifies the governing equations of the three physical fields into a single weak form under the ALE reference configuration.

[0088] The single-block variational form is as follows: Define function space Construct a unified variational problem: finding trial functions Make

[0089] For all test functions Established, among which For deformation gradient.

[0090] ( ) represents a variable under the reference configuration; for Space (the space of functions whose first derivatives are integrable by squares), used for velocity, displacement, and temperature; for Space (square-integrable function space) is used for pressure; Fluid domain under reference configuration; Solid domain under reference configuration; For the entire domain (the sum of the fluid domain and the solid domain); (•,•) For reference domain Inner product; Fluid velocity under reference configuration; Mesh velocity under reference configuration; Gradient operator for reference configuration; For the deformation gradient tensor, ( For grid displacement, I (unit tensor) The Cauchy stress tensor of the fluid domain under the reference configuration (after deformation gradient mapping); Pressure under reference configuration; For solid displacement; Solid stress under reference configuration (second Piola-Kirchhoff stress tensor); ρ f For fluid density; ρ s Density of solid; c p Specific heat capacity (constant pressure); k This is the thermal diffusivity (thermal conductivity). α The coefficient of thermal expansion; T 0 is the reference temperature; It is the vector of gravitational acceleration; b It is a force vector; Q Internal heat source items (including shear heat generation, radioactive decay heat, etc.); These are the test functions (shape functions) corresponding to velocity, pressure, displacement, and temperature, respectively.

[0091] Its terms comprehensively cover all mechanisms of multiphysics coupling: the fluid momentum term describes the time evolution of fluid velocity under the ALE framework; the fluid convection term includes grid motion correction to handle arbitrary moving boundaries; the fluid stress term and pressure divergence term jointly describe the momentum transfer and incompressible constraints within the fluid; the solid inertia term describes the dynamic response of the lithosphere under long-term tectonic loading; the solid stress term describes elastic or viscoelastic-plastic stress under finite deformation conditions; the thermal transient and convection term describes the evolution of the temperature field over time and thermal convection transport; the thermal diffusion term describes heat conduction; the buoyancy term describes the thermal drive generated by mantle temperature anomalies; the volume force term describes body force loads such as gravity; and the internal heat source term describes internal energy inputs such as radioactive decay heat and shear heat generation.

[0092] In this unified variational weak form, the velocity continuity and stress equilibrium conditions at the fluid-solid interface are no longer applied as explicit boundary constraints. Instead, these conditions are automatically satisfied in a variational sense by merging and eliminating area terms on the interface during the construction of the weak form. Specifically, during the derivation of the weak form, when integrating the momentum equations of the fluid and solid domains separately, an area term is generated at the fluid-solid interface, representing the virtual work done by the stresses acting on the interface on the fluid and solid sides in terms of virtual velocities or virtual displacements, respectively. When the weak forms of the fluid and solid equations are added together and the velocity continuity condition is applied, the stress vectors in these two area terms are equal in magnitude and opposite in direction, thus canceling each other out and naturally disappearing in the weak form without the need for explicit introduction of additional coupling mechanisms such as Lagrange multipliers or penalty functions. This approach embeds the interface coupling conditions within the variational mathematical structure, ensuring they are automatically satisfied during each Newton iteration step. This avoids the additional computational overhead of repeatedly exchanging interface data and iteratively coordinating as required by the partitioned coupling method, while also guaranteeing the strictness and consistency of the interface constraints. This provides a fundamental mathematical guarantee for the numerical accuracy and solution efficiency of the entire multiphysics simulation.

[0093] Spatial discretization: Hybrid finite element method: Q2 / Pdc element pair Q2 unit Material parameters can be assigned in zones, supporting special rheological properties in transition zones.

[0094] Specifically, in step S2, the fluid-solid interface is the Moho discontinuity, which is the discontinuity of seismic wave velocity between the crust and mantle, with a depth ranging from 35 km to 45 km. The fluid region represents the mantle or lower crust below the Moho discontinuity, while the solid region represents the lithosphere or upper crust above it. This interface setting directly corresponds the mathematical framework of fluid-solid coupling to the actual physical layering structure of the Earth's interior, giving the fluid and solid domains in the numerical model clear geological meanings. The Moho discontinuity is the most significant first-order velocity discontinuity in the Earth's interior. At this interface, the seismic P-wave velocity jumps from approximately 6.5 to 7.2 km / s at the bottom of the crust to approximately 7.8 to 8.2 km / s at the top of the mantle. This velocity jump reflects the fundamental changes in mineral composition and rock type between the upper and lower media. Above the interface, felsic and mafic rocks predominate, while below, ultramafic rocks such as peridotite predominate. Within the framework of lithospheric dynamics and mantle convection, the Moho (Mohorovičić discontinuity) is not only a material composition interface but also a transition interface for mechanical properties: the mantle material below the interface exhibits high-viscosity fluid behavior on a timescale of millions of years, with deformation mechanisms primarily consisting of diffusion creep and dislocation creep, and its thermal convection can be described by the incompressible Navier-Stokes equations under the Boussinesq approximation; the lithospheric crustal material above the interface exhibits elastic or elastoplastic solid behavior on the same tectonic timescale, capable of accumulating elastic strain energy under geological stress and releasing it through brittle fracture or ductile shear. In terrestrial lithospheric regions, the depth of the Moho typically ranges from 35 km to 45 km, a range that covers typical values ​​for normal terrestrial crustal thickness and is suitable for describing the crustal structure of tectonically stable zones, orogenic fronts, or intracontinental basins. Using the Moho depth as the interface parameter of the fluid-structure interaction computational domain enables the simulation to be customized for specific geological tectonic backgrounds, such as the periphery of the Tibetan Plateau, the North China Craton, or the Kunlun-Qilian Mountain tectonic belt. This provides a numerical experimental platform that conforms to the actual geophysical layered structure for studying how mantle convection drives the deformation of the overlying lithosphere through the drag force of the Moho and how the deformation of the lithosphere feeds back to the mantle convection by changing the boundary geometry and thermal state.

[0095] Furthermore, the velocity continuity condition requires that the fluid velocity equals the time derivative of the solid displacement, and a normal stress equilibrium condition is also applied at the Mohorovičić discontinuity (Moho), meaning that the projections of the fluid stress tensor and the solid stress tensor onto the interface normal direction are equal, with the interface normal direction being taken as the unit normal vector of the interface. These two conditions together constitute the complete coupling constraint at the fluid-solid interface, ensuring the coordinated transmission of physical field quantities on both sides of the interface from both kinematic and dynamic perspectives. The velocity continuity condition requires that the flow velocity of the mantle material at the Moho and the velocity of the solid material at the bottom of the lithosphere be equal everywhere on the interface, where the velocity on the solid side is represented by the time derivative of its displacement. The physical essence of this kinematic constraint is that during the process of mantle convection driving lithospheric deformation, the flow at the top of the mantle directly acts on the bottom of the lithosphere through viscous drag. The bottom of the lithosphere moves at the same velocity as the mantle flow, thus transferring the momentum of mantle convection to the overlying solid plate. Conversely, tectonic deformation of the lithosphere, such as the uplift of orogenic belts or the subsidence of basins, causes changes in the geometry and motion of the Mohorovičić discontinuity (Moho), which are also fed back to the fluid domain in real time through this velocity continuity condition. This allows the mantle convection pattern to adjust with the evolution of the lithosphere, and the two always maintain kinematic coordination at the interface, avoiding numerical non-physical phenomena such as fluid-solid separation or penetration. The normal stress equilibrium condition requires that the projections of the fluid stress tensor and the solid stress tensor onto the interface in the normal direction be equal, that is, the action and reaction forces exerted by the media on both sides of the interface in the direction perpendicular to the Moho satisfy Newton's third law. In specific scenarios of dynamic topographic simulation, the normal direction of the Moho discontinuity is usually not parallel to any coordinate axis, but rather exhibits a spatially varying tilt due to changes in crustal thickness and tectonic deformation. Therefore, stress equilibrium conditions need to be calculated by projection based on the local normal vector of the interface. The stress on the fluid side includes pressure terms and viscous deviatoric stress terms, while the stress on the solid side is determined by elastic strain and possible plastic yielding. The balance between the two on the normal projection ensures that the normal force of the mantle on the bottom of the lithosphere can be accurately applied as the force boundary condition of the solid domain to the bottom of the lithosphere. This allows the buoyancy and dynamic effects of mantle convection to be realistically transmitted to the interior of the lithosphere, providing a mechanically consistent numerical calculation basis for studying major geological processes involving the transmission of interface normal forces, such as mantle plume upwelling, lithospheric delamination, or subducting slab dragging.

[0096] Specifically, in step S2, updating the viscosity parameter includes updating the viscosity in real time based on the current strain rate, temperature field, and pressure distribution, according to at least one of the following: viscoelastic-plastic constitutive relations, the Drucker-Prager yield criterion, and strain softening or strain hardening laws. In geodynamic simulations, the rheological behavior of lithosphere and mantle materials is the core factor controlling tectonic deformation and mantle convection patterns. Viscosity, as a physical parameter characterizing a material's resistance to flow deformation, is not a constant under real geological conditions but varies significantly with physical conditions such as temperature, pressure, and strain rate. In the high-temperature environment of the mantle, the creep strength of major rock-forming minerals such as olivine is strongly dependent on temperature; an increase of several hundred degrees Celsius can lead to a decrease in effective viscosity by several orders of magnitude. Within the lithosphere, with increasing depth and confining pressure, the deformation mechanism of rocks gradually transitions from brittle frictional sliding near the surface to intragranular dislocation creep and diffusion creep in the deep interior. This transition is accompanied by a gradual change in effective viscosity from the equivalent high viscosity of brittle fracture to the low viscosity of ductile flow. The operation of updating the viscosity parameter is precisely to reflect the spatiotemporal changes of the aforementioned rheological properties in real time within the numerical model. In practice, within each Newton iteration step or each time step, the second invariant of the strain rate tensor is calculated based on the currently solved velocity field. Combined with the temperature field and pressure distribution, the effective viscosity of each integration point or each element is calculated by calling the specified constitutive relation model. For regions employing viscoelastic-plastic constitutive relations, the effective viscosity is determined by the product of the elastic shear modulus and the relaxation time, which is itself a function of temperature and pressure. For regions employing the Drucker-Prager yield criterion, when the stress state reaches the yield surface, the equivalent plastic viscosity is calculated back using the plastic flow law, allowing the apparent viscosity of the material to adjust with the accumulation of plastic strain, thereby simulating the irreversible deformation behavior of rocks under tectonic stress. For fault zones with introduced strain softening or strain hardening laws, as the accumulated plastic strain increases, yield parameters such as cohesion and internal friction angle evolve according to preset decay or growth curves, correspondingly updating the ultimate strength and equivalent viscosity of the material. This causes the fault to undergo a softening stage with reduced strength after initial activation, thus spontaneously forming strain localization and continuous slip. This process plays a crucial mechanical control role in simulating the complete cycle of a fault from locking to unlocking and then to seismic slip. Through the aforementioned viscosity update mechanism, the entire simulation method can naturally transition between the friction intensity-dominated region of the brittle upper crust, the creep intensity-dominated region of the ductile lower crust, and the convective viscosity-dominated region of the mantle within a unified computational framework. This enables the dynamic topographic evolution model to possess the ability to describe a continuous medium from the surface to the deep mantle in terms of material mechanical behavior.

[0097] In addition, in step S2, the Newton iteration framework uses the incomplete Newton method, and the Jacobian matrix is ​​updated every 5 to 10 Newton iteration steps.

[0098] Newton's solution framework: An incomplete Newton method is used, updating the Jacobian matrix every 5-10 steps. Line search strategy: Determined by the Armijo criterion. Convergence criterion: and . Indicates the first Second Newton iteration Represents the residual.

[0099] in, X k+1 For the first k The unknown quantity after +1 Newton iteration update. X k For the current number k The unknown vector of the next iteration, Δ X k Solve the system of linear equations for the current iteration step (e.g.) J Δ X =- R The correction amount (Newton direction) obtained. α k The line search step size factor (typically between 0 and 1), determined by the Armijo criterion, is used to control the magnitude of the advance along the correction direction to ensure that the residual function decreases sufficiently. k This is for Newton's iteration counting. F ( X k ) is the current number k The global residual vector obtained by assembling the discrete equations of each physical field (such as fluid, solid, and thermodynamic equations) in each iteration. The norm of a vector (usually the 2 norm). This is the preset absolute residual convergence tolerance (used to determine the absolute precision of the equation equilibrium). This is the preset relative residual convergence tolerance (used to measure the magnitude of the correction relative to the current solution). Let the norm of the current correction vector be denoted as . The norm of the current unknown vector, together with the norm of the unknown vector, constitutes the convergence criterion: Newton iteration is considered convergent only when the residual norm is less than the absolute tolerance and the correction norm is less than the current solution norm multiplied by the relative tolerance. k For Newton's iteration counting, Generally, this refers to the preset tolerance threshold.

[0100] In the nonlinear solution of fluid-solid-thermal coupling, each complete Newton iteration step requires calculating the Jacobian matrix, which involves taking the partial derivatives of all residual terms of the multiphysics system with respect to all unknowns. This operation involves calculating the cross-derivatives between the Navier-Stokes equations in the fluid domain, the finite deformation equations in the solid domain, and the global thermodynamic equations. Given the coupling of multiphysics variables, the computational cost is far higher than that of constructing the Jacobian matrix for single-physics problems. Adopting an incomplete Newton method strategy, where the Jacobian matrix remains unchanged across multiple Newton iteration steps and is only recalculated every 5 to 10 steps, effectively reduces the number of Jacobian matrix assembly steps and the associated numerical integration and matrix assembly overhead without significantly affecting the nonlinear convergence speed. Between two Jacobian matrix updates, the Newton iteration still uses the previously calculated Jacobian matrix to generate the correction direction. Although this matrix is ​​not the exact derivative at the current iteration point, since the change in the Jacobian matrix is ​​relatively slow when the nonlinear iteration of the multiphysics system has entered the local convergence stage, this approximation can still provide a sufficiently accurate descent direction. The update interval is set to a range of 5 to 10 steps, determined through a convergence trade-off within the typical parameter range of geodynamic multiphysics simulations. Too short an interval offers limited improvement in computational efficiency, while too long an interval may lead to a significant decrease in the convergence speed of the Newton iteration due to outdated Jacobian matrix information, or even divergence. By introducing this incomplete Newton strategy, the computational cost of the entire long-term evolution simulation is effectively distributed and reduced, further ensuring the practical feasibility of large-scale three-dimensional dynamic topographic thermodynamic simulations under limited computational resources.

[0101] In addition, the Newton iteration step size is adjusted to determine the line search step size using the Armijo criterion. The Armijo criterion is a line search method that obtains sufficient descent of the objective function by gradually reducing the trial step size, so that the updated residual satisfies the preset sufficient descent condition.

[0102] In multiphysics nonlinear block coupling iterations, although the correction direction generated by the Newton method has a local second-order convergence rate, in the initial stage of iteration or when the system is highly nonlinear, the full-step Newton correction may cause the residual function value to increase instead of decrease, or even cause iteration divergence. This problem is particularly prominent in fluid-solid-thermal strongly coupled systems. The sensitivity of the mantle convection velocity field to the temperature field, the nonlinear dependence of solid displacement on the interface drag force, and the variation of viscosity parameters with strain rate and temperature may all cause complex non-convex shapes of the residual surface in local regions. The Armijo criterion addresses the above problem by introducing a backtracking mechanism that gradually reduces the trial step size. In each Newton iteration step, the new residual is first calculated by trial with the full step size. If the new residual satisfies the preset sufficient descent condition, that is, the reduction of the residual is not less than the line search step size multiplied by the projection of the residual gradient in the correction direction multiplied by a preset constant, then the step size is accepted; if it is not satisfied, the step size is reduced by a preset ratio and the trial is repeated until a step size that satisfies the descent condition is found or the step size is reduced to a preset lower limit value. This backtracking process ensures that each Newton iteration step produces a substantial reduction in residuals, preventing iterative oscillations or divergence caused by excessively large single corrections crossing local minima on the residual surface. In the application scenario of this invention, since the Jacobian matrix is ​​updated only once every multiple steps using the incomplete Newton method, the matrix information may lag behind the current iteration point, resulting in a decrease in the quality of the correction direction compared to the standard Newton method. In this case, the Armijo line search becomes more crucial. By adaptively scaling the correction step size, it compensates for the correction direction deviation that may be caused by the lag in Jacobian matrix information, ensuring the overall convergence robustness of the nonlinear coupled iteration.

[0103] Furthermore, in step S3, the adaptive determination of the step size of the next time step is achieved through three sequentially executed steps, which combines the physical numerical stability constraints with the time scale characteristics of the simulation process to form a step size control mechanism that adapts to the simulation requirements of geodynamics over a time span of millions of years.

[0104] like Figure 3 As shown, the time-step adaptive advancement control logic of the method of the present invention is illustrated. After the residual converges, the process enters the step size control stage, which includes calculating the candidate step size based on the numerical stability condition (S31), judging the order of magnitude interval and applying logarithmic time scale adjustment (S32), and updating the step size (S33). The next time step is prepared by using velocity extrapolation, forming a closed-loop iteration.

[0105] Specifically, step S3 adaptively determines the step size of the next time step, including: S31: Calculate candidate step sizes based on numerical stability conditions.

[0106] In the time-progression process of multiphysics coupled simulations, the step size is first strictly constrained by the physical propagation velocity. Various physical processes with different characteristic velocities exist in the geodynamic system: mantle convection typically occurs at speeds of several centimeters per year, corresponding to convection timescales of millions of years; while convective transport of the temperature field near the thermal boundary layer can generate significant temperature gradient changes in a relatively short time; heat source terms such as shear heat generation and radioactive decay have their own characteristic response times. Step S31 comprehensively considers the characteristic velocities of these physical processes and, based on numerical stability conditions, ensures that the distance the fluid travels within a time step does not exceed the characteristic size of the grid cell, and that the propagation distance of heat diffusion within a time step does not exceed the spatial variation scale of the temperature field, thus calculating a candidate step size that meets the numerical stability requirements. This candidate step size is the maximum allowable step size that ensures the numerical solution does not diverge under the current physical field conditions, providing an upper limit benchmark determined by physical constraints for subsequent step size adjustments.

[0107] S32: Determine the order of magnitude range of the current simulation time, apply logarithmic time scale adjustment to the candidate step size according to the order of magnitude range, and obtain the final step size of the next time step to be adopted.

[0108] A significant characteristic of geodynamic simulations is that their total simulation duration spans multiple orders of magnitude, from hundreds of years for earthquake cycles to tens of millions of years for orogenic belt evolution. The characteristic timescales of physical processes fundamentally change with different stages of tectonic evolution. If a uniform step size strategy is adopted throughout the simulation, the time step may be too small in the early stages, wasting computational resources, or too large in the later stages, failing to identify key geological events. The logarithmic timescale adjustment mechanism in step S32 is designed to address this problem: the total simulation duration is divided into multiple order-of-magnitude intervals. As the simulation progresses from one interval to the next, the candidate step size is adjusted according to a preset mapping relationship, ensuring that the logarithm of the step size is approximately linearly related to the simulation progress. This adjustment method allows the step size to grow reasonably as the simulation progresses. In the early stages, a smaller step size accurately identifies rapid evolutionary processes such as the initiation of mantle convection and initial topographic uplift; in the middle and later stages, a larger step size efficiently advances long-term tectonic evolution, thus achieving a dynamic balance between computational accuracy and efficiency that adapts to the characteristics of the geological timescale.

[0109] S33: Update the current time step with the final adopted next time step size, and use it as the time advance step size for the next time step cycle.

[0110] The final step size, determined after logarithmic time-scale adjustment, is assigned to the system's time step variable, replacing the step size value used in the previous time step, and becomes the time increment driving the advancement of all physical fields in the next time step. This update operation makes the entire simulation's time step sequence form a curve that adaptively adjusts with the simulation process. It strictly satisfies the numerical stability condition within each time step, and conforms to the inherent characteristic of geological processes spanning multiple orders of magnitude in time scale. This provides numerical feasibility and computational efficiency guarantees at the time integration level for simulating dynamic topographic evolution on a million-year scale.

[0111] Furthermore, step S31 calculates the candidate step size based on the numerical stability condition. This candidate step size is not the result of calculation under a single criterion, but a safe and reasonable upper limit of time step size determined after considering multiple control factors.

[0112] Specifically, step S31 calculates the candidate step size based on the numerical stability condition, including: S311: Obtain the preset maximum allowable step size, the current step size amplified by the growth factor, and the upper limit of the stable step size determined based on the Courant-Friedrich-Liouvi conditions.

[0113] These three factors constrain the time step from different perspectives. The preset maximum allowable step size is an upper limit set by the simulation designer before computation begins. Its purpose is to prevent insufficient time resolution due to excessively large step sizes in regions with slow physical process evolution or finer meshing, while also providing a controllable safety boundary for the entire simulation. The current step size amplified by a growth factor is obtained by multiplying the actual step size used in the previous time step by a preset growth factor greater than 1. This mechanism allows the step size to gradually increase during a continuous and stable solution phase, avoiding drastic jumps between adjacent time steps, thus maintaining the accuracy of the time discretization scheme and the convergence stability of the solver. The stable step size upper limit determined based on the Courant-Friedrich-Liouvi conditions is a step size constraint calculated from the physical essence of numerical solutions to partial differential equations, based on the maximum characteristic velocity of the current velocity field and the minimum characteristic size of the mesh element. This condition requires that the distance of fluid motion or thermal convection propagation within a time step does not exceed the scale of one mesh element; otherwise, the numerical scheme will lose stability and produce non-physical numerical oscillations.

[0114] S312: Take the minimum value among the preset maximum allowable step size, the amplified current step size, and the upper limit of the stable step size as the candidate step size.

[0115] Adaptive time step control: Long-term integration uses a logarithmic time scale, supporting simulations over millions of years. Represents a limiting factor, which can be taken as... The values ​​between.

[0116] Where, Δ t n+1 The candidate step size for the next time step is Δ, and min() is the minimum function used to synthesize multiple constraints to ensure numerical stability. t max The preset maximum allowable step size is a safety limit set by humans to prevent the step size from increasing indefinitely when the physical process evolves slowly, which would lead to insufficient time resolution. γ Δ is the step size growth factor (i.e., the limiting coefficient), with a value between 0 and 1. It is used to proportionally increase the step size of the previous time step, ensuring a smooth transition between adjacent time steps and preventing numerical oscillations caused by abrupt changes in step size; t n Δ is the actual step size used in the current time step (previous time step); t CFL The upper limit of the stable step size is calculated based on the Courant-Friedrich-Liouvi (CFL) numerical stability condition, which is the maximum allowable step size that satisfies the condition that the fluid motion within a time step does not exceed the size of one grid cell under the physical propagation and thermal diffusion characteristic velocities.

[0117] Taking the minimum of the above three factors as the candidate step size is essentially a conservative and robust step size determination strategy: when the physical process is intense and the flow rate is high, the step size constrained by the Courant-Friedrich-Liouvi conditions is usually the most stringent and plays a dominant role; when the physical process tends to be stable and the flow rate decreases, the current step size amplified by the growth factor or the preset maximum allowable step size becomes the actual constraint; taking the minimum of the three factors ensures that the candidate step size simultaneously satisfies multiple constraints of physical stability, numerical smoothness, and safety boundary at any time, providing a numerically safe and reliable, and physically reasonable and self-consistent benchmark step size for subsequent logarithmic time scale adjustment.

[0118] Furthermore, step S32, based on the fact that the candidate step size has met the numerical stability condition, introduces an adjustment mechanism related to the time scale of the simulation process, so that the final adopted time step size can adapt to the long-term integration requirements of geodynamic simulations spanning multiple orders of magnitude. This adjustment mechanism is achieved by first dividing the time step into orders of magnitude intervals, then determining the interval switching and applying a mapping adjustment.

[0119] Specifically, in step S32, the order-of-magnitude range of the current simulation time is determined, and a logarithmic time scale adjustment is applied to the candidate step size based on the order-of-magnitude range to obtain the final step size for the next time step, including: S321: Divide the total simulation duration into multiple order-of-magnitude intervals.

[0120] In long-term geodynamic evolution simulations, the typical total time span can extend from hundreds of years to tens of millions of years or even longer, with the characteristic timescales of physical processes fundamentally changing depending on the stage of tectonic evolution. For example, in the initial stage of mantle convection initiation, the thermal boundary layer is still in the adjustment period, the spatial gradient of the temperature field changes rapidly, and the establishment of the convection velocity field involves a relatively short timescale, usually hundreds to thousands of years. After entering the stable convection stage, the changes in the mantle convection pattern tend to slow down, and the characteristic timescale rises to hundreds of thousands to millions of years. In the lithospheric tectonic response stage, such as the isostatic adjustment of orogenic belt uplift and erosion or the subsidence and filling process of basins, the timescale may again be between hundreds of thousands and tens of millions of years. To address this multi-scale characteristic, during simulation initialization, the entire geological time span to be simulated is pre-divided into several orders of magnitude intervals on a logarithmic scale, for example, each interval spanning an order of magnitude corresponds to a specific step size mapping relationship. This division allows the step size adjustment strategy to distinguish between the need to capture rapid physical processes with fine resolution in the early stages of the simulation and the need to efficiently advance long-term tectonic evolution with large step sizes in the later stages of the simulation, providing a structured interval framework for subsequent adaptive adjustment.

[0121] S322: Determine the order of magnitude interval of the current simulation time. When the simulation process moves from the previous order of magnitude interval to the next order of magnitude interval, adjust the candidate step size according to the preset interval step size mapping relationship so that the logarithm of the adjusted step size is linearly related to the simulation process, and obtain the final step size of the next time step.

[0122] This step performs a judgment once during the determination of the step size at each time step. Specifically, it obtains the current simulation time value, identifies its order-of-magnitude interval, and compares this interval with the interval corresponding to the previous time step. If the current simulation time is still within the same order-of-magnitude interval as the previous time step, the candidate step size remains unchanged. If the simulation process has crossed the interval boundary and entered the next order-of-magnitude interval, the candidate step size is adjusted according to the step size mapping relationship corresponding to the new interval, so that the logarithm of the adjusted step size is approximately linearly related to the simulation progress. This logarithmic linear relationship means that as the simulation progresses from the early to the middle and late stages, the step size does not increase continuously without constraints, but rather increases gradually according to a logarithmic law. The time increment of each step matches the current geological timescale, which avoids missing key physical processes such as convection initiation or initial deformation due to excessively large step sizes in the early stages of the simulation, and also avoids wasting computational resources by continuing to advance with small step sizes in the later stages of the simulation. The step size obtained after the above logarithmic time scale adjustment is the final step size of the next time step, which will be used as the input value for updating the current time step in step S33 to complete the time advancement from the current time step to the next time step.

[0123] Furthermore, numerical solutions for multiphysics coupled simulations of geodynamics face a variety of instabilities from different sources, including error accumulation from long-term integration, algebraic ill-conditioning caused by fluid incompressibility constraints, and numerical oscillations induced by heat transport dominated by strong convection. A single stabilization method cannot simultaneously cover all of these problems. Therefore, this method integrates three complementary numerical stabilization mechanisms during the time-progression process.

[0124] During the time progression in steps S1 to S3, a multi-mechanism numerical stabilization process is applied. The multi-mechanism numerical stabilization process includes: discretizing the time derivative term using a backward second-order difference scheme in step S2; performing pressure regularization when solving the fluid domain equation in step S2; and applying artificial viscosity based on entropy viscosity when solving the thermodynamic equation in step S2.

[0125] In step S2, the time derivative term is discretized using a backward second-order difference scheme. During time progression, the choice of the time discretization scheme directly affects the accuracy and stability of long-time integration. The backward second-order difference scheme is an implicit multi-step method that constructs a second-order accurate approximation of the time derivative using the function values ​​at the current time, the previous time, and two time points prior. Compared to the implicit Euler scheme, it has higher truncation accuracy; compared to the Crank-Nicolson scheme, it exhibits better numerical damping characteristics when dealing with rigid systems, effectively suppressing high-frequency spurious oscillations caused by nonlinear coupling effects during long-time integration. For geodynamic simulations spanning millions of years, the numerical damping effect of this scheme helps filter out physically unreliable short-period fluctuations while preserving the true evolution signal over the tectonic timescale.

[0126] In step S2, pressure regularization is performed when solving the fluid domain equations. When the fluid domain adopts the incompressible Navier-Stokes equations under the Boussinesq approximation, the discretized linear system has a saddle-point structure. Pressure, acting as a Lagrange multiplier, forces the incompressible constraint, resulting in a zero-valued diagonal block in the coefficient matrix. This significantly increases the condition number of the entire system, severely reducing the convergence speed of the iterative solver. Pressure regularization scales the pressure variable using the ratio of fluid viscosity to characteristic length as a scaling factor, making the pressure-related sub-blocks numerically match the velocity-related sub-blocks. This improves the spectral properties of the Stokes subsystem, reduces the condition number, and allows the iterative solver to achieve convergence accuracy in fewer iterations.

[0127] Pressure regularization: Pressure is amplified using viscosity and mantle plume diameter L, scaling = Improve the condition number of the Stokes system.

[0128] In step S2, an artificial viscosity treatment based on entropy viscosity is applied when solving the thermodynamic equations. When thermal convection dominates, the temperature field may exhibit steep gradients or near-discontinuous characteristics in local regions. Standard Galerkin finite element discretization in such convection-dominated problems will produce non-physical numerical oscillations, leading to negative values ​​or spurious spatial fluctuations in the temperature field. These errors then propagate to the momentum equation through the buoyancy term, severely compromising the accuracy and stability of the entire coupled solution. The artificial viscosity treatment based on entropy viscosity automatically identifies convection-dominated regions by calculating the entropy residual of the temperature field. Using the local entropy yield as an indicator of oscillation intensity, artificial viscosity is adaptively introduced at the required locations, proportional to the mesh feature size and the local entropy residual. In regions dominated by diffusion or with gentle temperature changes, the artificial viscosity automatically decays to near zero, thereby effectively suppressing convective oscillations without introducing excessive numerical dissipation, maintaining the physical realism and numerical stability of the temperature field. The three stabilization methods mentioned above are applied to the three levels of time discretization, algebraic solution, and physical field discretization, respectively, and work together to ensure the overall numerical robustness of multi-physics coupling simulation with a time span of millions of years.

[0129] The discrete form of the temperature field time derivative of the backward second-order difference scheme involves the temperature fields at the current time, the previous time, and the previous two time points. In step S1, the temperature extrapolation adopts an extrapolation formula that matches the backward second-order difference scheme. The predicted temperature field of the current time step is constructed by linear combination of the temperature fields at the previous two known time points.

[0130] The backward second-order difference scheme is an implicit linear multistep method that constructs a second-order accurate approximation of the time derivative using the temperature values ​​of three time layers: the current time, the previous time, and the two time layers before that, rather than using only the first-order difference between two adjacent time layers. Therefore, when discretizing multi-timescale thermal evolution processes such as mantle convection and lithosphere heat conduction, it can simultaneously provide high truncation accuracy and good numerical damping characteristics.

[0131] Time discretization: using a backward second-order difference scheme (BDF2) + two-step extrapolation: .

[0132] / This represents the partial derivative of temperature with respect to time (i.e., the time derivative term). n Indicates the sequence number of the current time step. n -1 indicates the sequence number of the previous moment (the previous time step). n -2 represents the sequence number of the two previous time steps, therefore The temperature field at the current moment, The temperature field at the previous moment. This refers to the temperature field at the first two moments; here...k Represents the time step (usually denoted by Δ in the document). t The time interval (represented by the time increment between the current step and the previous step) is the time increment between the current step and the previous step.

[0133] The numerical damping of this scheme helps to attenuate high-frequency pseudo-oscillations caused by nonlinear coupling effects during long-term integration, while preserving low-frequency physical signals that match the construction evolution timescale. Correspondingly, the temperature extrapolation in step S1 uses an extrapolation formula of the same order as the backward second-order difference scheme, calculating the predicted temperature field for the current time step through a linear combination of the temperature fields from the previous two known time points. Both schemes maintain consistency in data dependency and truncation order. The discrete scheme uses temperature values ​​from three time levels to construct a derivative approximation, while the extrapolation formula uses temperature values ​​from two known time points to provide an initial estimate for the third time point through linear extrapolation. This matching design ensures that the predicted temperature field has time accuracy compatible with the discrete scheme before entering the Newton iteration solution, effectively reducing the residual level in the initial stage of nonlinear iteration and reducing the number of iterations required to reach convergence.

[0134] Pressure regularization improves the condition number of the Stokes system by scaling the pressure variable using the ratio of fluid viscosity to characteristic length as a scaling factor, with the scaling factor on the order of 10. 17 The characteristic length is the characteristic diameter of the mantle plume.

[0135] Pressure regularization improves the condition number of the Stokes system by scaling the pressure variable using the ratio of fluid viscosity to characteristic length as a scaling factor. This scaling factor is on the order of 10^17, and the characteristic length is taken as the characteristic diameter of the mantle plume. In each Newton iteration step of a geodynamic multiphysics simulation, the algebraic subsystem formed after discretizing the fluid domain exhibits a typical saddle-point structure. The submatrices between velocities originate from the discretization of the elliptic operator of the viscosity term, possessing good numerical properties. However, pressure, as a Lagrange multiplier, is forced to satisfy incompressible constraints, and its corresponding diagonal sub-block is a zero matrix. This significantly increases the condition number of the entire Stokes subsystem, causing the iterative solver to converge slowly or even stall when solving such systems. The principle of pressure regularization is to introduce a scaling factor with physical dimensions to make the pressure variable dimensionless or normalize its magnitude, so that the pressure-related sub-blocks match the velocity-related sub-blocks in numerical magnitude. This improves the spectral distribution and condition number of the coefficient matrix, reducing the difficulty of iterative solutions. The scaling factor is taken as the ratio of fluid viscosity to characteristic length, with a value on the order of 10 to the power of 17. This value is not arbitrarily chosen, but is naturally determined by typical physical parameters of the mantle. The effective viscosity of mantle rocks is on the order of 10¹⁹ to 10²¹ Pascals per second, while the characteristic diameter of mantle plumes is on the order of tens to hundreds of kilometers. Dividing the two naturally falls into 10. 17The characteristic length of the mantle plume is near the order of magnitude. Using the characteristic diameter of the mantle plume as the characteristic length has clear physical significance: the mantle plume is a vertically rising flow structure driven by thermal buoyancy in the mantle convection system. Its diameter represents a typical spatial scale for thermal boundary layer instability and upflow focusing, and it is also the region where the mantle convection velocity and temperature fields change most dramatically. Pressure regularization at this scale can most effectively improve the numerical properties of the Stokes subsystem associated with the mantle plume. After pressure regularization, the condition number of the Stokes subsystem can be reduced by several orders of magnitude, and the number of iterations required by the iterative solver under the same convergence tolerance is significantly reduced, providing an algebraic solution-level efficiency guarantee for the engineering implementation of large-scale three-dimensional mantle convection simulations.

[0136] Artificial viscosity treatment based on entropy viscosity includes: constructing entropy viscosity on the cell using the entropy residual of the temperature field and the cell characteristic size, and taking the smaller value between the local maximum permissible viscosity and this entropy viscosity as the applied artificial viscosity to suppress numerical oscillations in the region dominated by thermal convection.

[0137] Entropy viscosity method: Artificial adhesiveness: This suppresses numerical oscillations caused by strong convection. Among them, v Indicates viscosity, Indicates the current unit K The entropy viscosity on, where the subscript h The superscript indicates that the grid scale is relevant. E Representing entropy as the metric, vertical lines K Limited to the unit K Take the value above; c R It is an empirical constant (usually taken as a dimensionless constant of the order of O(1)). h K For unit K Feature dimensions (such as mesh cell diameter or side length); R E ( T ( ) represents the temperature field T The entropy residual is used to measure the degree to which the numerical solution deviates from the physical true solution within a local grid cell by calculating the transient change of entropy and the balance deviation between convection transport and diffusion. For unit K On Norm (i.e., the maximum absolute value within the given region); E ( T Let be the entropy function corresponding to the temperature field; Here, Ω represents the entropy calculated at the reference temperature; Ω represents the entire computational domain. Represents the global domain (That is, the maximum absolute value within the entire region). In the artificial viscosity formula... middle, For the current unit K The preset local maximum permissible viscosity, The entropy viscosity is calculated above. This is the artificial viscosity that is ultimately applied to the thermodynamic equation, taking the smaller of the two values.

[0138] Artificial viscosity based on entropy viscosity suppresses non-physical oscillations in convection-dominant regions by adaptively introducing local numerical dissipation when solving thermodynamic equations. In mantle thermal convection simulations, when convective transport effects are much stronger than thermal conduction, the temperature field forms extremely steep gradients near the thermal boundary layer and rising or falling streams. Standard finite element discretization in such regions produces spurious numerical fluctuations, manifesting as temperature overshoot, negative temperature values, or high-frequency spatial oscillations. These non-physical signals, coupled back to the momentum equation through the buoyancy term, lead to a decrease in the accuracy of the entire coupled solution or even divergence. The entropy viscosity method, based on the second law of thermodynamics, treats the temperature field as a state variable of a thermodynamic system, defines its entropy function, and measures the degree to which the numerical solution deviates from the physical true solution within a local grid cell by calculating the transient change of entropy and the balance residual between convective transport and diffusion. This entropy residual approaches zero in diffusion-dominant regions with gentle temperature changes, but increases significantly in convection-dominant regions with sharp temperature gradients and insufficient grid resolution to capture the convection front, thus providing a rigorously physical adaptive indicator for applying artificial viscosity. For each cell, the entropy viscosity is constructed by multiplying the modulus of the entropy residual by the square of the cell's characteristic size and dividing by the magnitude of the global entropy change. To avoid over-application of artificial viscosity in local regions and masking the actual physical dissipation, the actual applied artificial viscosity is taken as the smaller of the preset local maximum allowable viscosity for that cell and the calculated entropy viscosity. This provides just enough minimum dissipation to suppress oscillations in areas of intense convection, while automatically degenerating to zero additional viscosity in areas of diffusion dominance. This treatment preserves the true heat transport characteristics of mantle thermal convection to the greatest extent possible while ensuring numerical stability, providing a physically self-consistent numerical stabilization guarantee for the simulation of thermodynamic evolution on a million-year scale.

[0139] like Figure 4 As shown, Figure 4The process starts with the "initial temperature field". After loading the initial temperature field, the process proceeds to "temperature extrapolation", which uses the temperature field information of the previous time step or the two previous time steps to make a preliminary prediction of the temperature field at the current time step. After the temperature extrapolation is completed, the process enters "solving the Stokes equation", that is, using the incompressible Navier-Stokes equation under the Boussinesq approximation to solve the momentum conservation and mass conservation problems of the fluid domain, namely the mantle / lower crust, to obtain the velocity field and the pressure field. Then the "velocity coupling condition" is executed, under the arbitrary Lagrangian-Eulerian framework, the fluid velocity and the deformation rate of the solid domain, namely the lithosphere / upper crust, are matched and constrained on the fluid-solid interface, namely the Moho surface, to ensure continuous interface velocity and force transfer. Next, "solving the solid equation" is performed to solve the finite deformation elastic or viscoelastic-plastic equation of the lithosphere under the action of gravity and fluid drag force, to obtain the solid displacement field. After solving the solid equation, the "displacement coupling condition" is executed to map the solid deformation to the fluid grid boundary; then "updating the grid" is performed. After completing the grid update, the process enters the "residual < tol?" judgment: if the residual does not meet the tolerance, "line search" is executed, then "reducing the incremental step size" is executed, and the process returns to "solving the Stokes equation" to continue the iteration; if the residual meets the tolerance, "velocity extrapolation" is executed, then "grid velocity" is calculated, then "solving the temperature equation" is performed, then "updating the viscosity" is performed, then "updating the CFL time step" is performed, and finally "temperature extrapolation" is executed to enter the next time step cycle, until the simulation ends, and the dynamic topography evolution result is output.

[0140] Wherein, according to Figure 4 the sequence shown, "solving the solid equation" follows the "velocity coupling condition"; "displacement coupling condition" and "updating the grid" follow "solving the solid equation"; "updating the grid" is located before the "residual < tol?" judgment; "solving the temperature equation" and "updating the viscosity" are located after the residual converges, and "solving the temperature equation" is after "grid velocity" and before "updating the viscosity"; "updating the CFL time step" follows "updating the viscosity"; "temperature extrapolation" follows "updating the CFL time step", and the coupled solution of the next time step is entered from the "temperature extrapolation". All the above processes are based on Figure 4 the standard, and expressions inconsistent with Figure 4 such as "velocity coupling condition → updating viscosity", "solving solid equation → solving temperature equation", and "performing displacement coupling and grid updating after convergence" are no longer adopted.

[0141] Accordingly, a second aspect of the present invention provides an electronic device, including: at least one processor; and a memory connected to the at least one processor; wherein the memory stores instructions executable by the at least one processor, the instructions being executed by the at least one processor to cause the at least one processor to perform the above-described dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework.

[0142] Accordingly, a third aspect of the present invention provides a computer-readable storage medium having computer instructions stored thereon, which, when executed by a processor, implement the above-described dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework.

[0143] The embodiments of this invention aim to protect a dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework, which has the following effects: 1. By identifying the block sparse structure of multi-physics coupled algebraic equation systems, the unknowns are divided into a first variable group containing fluid-solid interface displacement and temperature and a second variable group containing fluid velocity, pressure and other solid displacements according to the type of physical field. Geometric multigrid and algebraic multigrid preprocessing are applied to the sub-matrix blocks corresponding to different physical fields, which significantly accelerates the convergence speed of large-scale nonlinear systems and achieves higher computational efficiency than general solvers. This solves the problem of low computational efficiency caused by the failure to fully utilize the system matrix structure characteristics when dealing with fluid-solid-thermal three-way coupling problems in existing technologies. 2. By constructing the fluid dynamics equations, solid mechanics equations, and thermodynamics equations into a single-block coupled variational model with fluid velocity, pressure, solid displacement, and temperature as unknowns within a unified arbitrary Lagrange-Euler framework, the velocity continuity and stress equilibrium conditions at the fluid-solid interface are automatically satisfied in this model. This enables a closed-loop coupled simulation of thermal driving, force transmission, and deformation feedback during mantle convection-driven lithospheric deformation, fundamentally overcoming the shortcomings of existing technologies that cannot fully characterize the thermo-mechanical-fluid feedback mechanism due to separating fluid, solid, and thermodynamic processes. 3. By combining candidate step size calculation based on numerical stability conditions with logarithmic time scale adjustment during the time progression, the step size is adaptively adjusted according to the order of magnitude range of the simulation process. At the same time, backward second-order difference scheme time discretization, pressure regularization processing and artificial viscosity processing based on entropy viscosity are applied to ensure the numerical stability and accuracy of the integration process with a time span of millions of years. This solves the technical problem that existing technologies cannot balance long-term integration efficiency and numerical stability of multi-physics coupling.

[0144] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product embodied on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code.

[0145] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0146] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0147] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0148] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention and not to limit it. Although the present invention has been described in detail with reference to the above embodiments, those skilled in the art should understand that modifications or equivalent substitutions can still be made to the specific implementation of the present invention. Any modifications or equivalent substitutions that do not depart from the spirit and scope of the present invention should be covered within the scope of protection of the claims of the present invention.

Claims

1. A dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework, characterized in that, Includes the following steps: Step S1: Load the initial temperature field of the target geological area. At the beginning of each time step, perform temperature extrapolation using the temperature fields of the previous two known times to obtain the predicted temperature field of the current time step. The temperature extrapolation is an operation that uses the temperature field at the known time as input and constructs an approximate value of the temperature field at the current time through time extrapolation. Step S2: Using the predicted temperature field as the initial thermodynamic value, under the same Newton iteration framework, sequentially solve the momentum and mass conservation equations in the fluid domain, update the viscosity parameters, solve the finite deformation equations in the solid domain, and solve the global thermodynamic equations. Apply the velocity continuity condition at the fluid-solid interface to form a nonlinear block coupling iteration that solves all physical field variables jointly. Step S3: Determine whether the residual of the current iteration has converged. If it has not converged, adjust the Newton iteration step size and return to step S2 to continue iterating. If it has converged, update the mesh and field quantity according to the converged solution, adaptively determine the step size of the next time step, perform velocity extrapolation and enter the next time step loop until the simulation ends, and output the dynamic terrain evolution results.

2. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 1, characterized in that, The joint solution of all physical field variables described in step S2 includes: S21: Identify the block sparse structure of the algebraic equation system in the nonlinear block coupling iteration, and divide the unknowns into the first variable group and the second variable group according to the physical field type; S22: Construct the Schur complement matrix for the divided variable groups to form a reduced system with the first variable group as the unknown; S23: Solve the reduced system to obtain the correction amount for the first set of variables; S24: Substitute the correction amount of the first variable group back into the second variable group to calculate the correction amount of the second variable group; S25: Update all physical field variables together with the correction values ​​of the first variable group and the second variable group to complete this joint solution.

3. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 2, characterized in that, Step S21, which involves identifying the block sparse structure of the algebraic equation system in the nonlinear block-coupled iteration, includes: S211: The algebraic equations are divided into blocks according to the four physical field variables of fluid velocity, pressure, solid displacement and temperature to obtain a block sparse matrix; S212: Identify the symmetric positive definiteness of the fluid velocity submatrix and the block diagonality of the solid displacement submatrix in the block sparse matrix, wherein the symmetric positive definiteness is used to guide the selection of a geometric multigrid preprocessor and the block diagonality is used to guide the selection of an algebraic multigrid preprocessor.

4. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 2, characterized in that, The method of classifying by physical field type in step S21 is as follows: fluid-solid interface displacement and temperature are used as the first variable group, and fluid velocity, pressure and solid displacement other than interface displacement are used as the second variable group. This allows the key physical quantities controlling fluid-solid coupling and thermo-mechanical coupling to be assigned to the first variable group, and the subordinate physical quantities that can be driven by the first variable group through coupling relationship to be assigned to the second variable group.

5. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 2, characterized in that, In step S221, the submatrix of the second variable group is constructed using a hybrid multigrid preprocessing method to approximate its inverse. S221a: Apply geometric multigrid preprocessing to the sub-matrix block corresponding to fluid velocity and pressure to obtain the approximate inverse of the sub-matrix block; S221b: Apply algebraic multigrid preprocessing to the sub-matrix block corresponding to the solid displacement to obtain the approximate inverse of the sub-matrix block; S221c: Assemble the approximate inverses obtained in the two steps according to the block diagonal structure to form the overall approximate inverse of the second variable group submatrix.

6. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 1, characterized in that, The adaptive determination of the step size for the next time step in step S3 includes: S31: Calculate candidate step sizes based on numerical stability conditions; S32: Determine the order of magnitude range of the current simulation time, apply logarithmic time scale adjustment to the candidate step size according to the order of magnitude range, and obtain the final next time step size to be adopted; S33: Update the current time step with the final adopted next time step size, and use it as the time advance step size for the next time step cycle.

7. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 6, characterized in that, Step S32, which involves determining the order-of-magnitude range of the current simulation time and applying a logarithmic time scale adjustment to the candidate step size based on the order-of-magnitude range to obtain the final next time step size, includes: S321: Divide the total simulation duration into multiple order-of-magnitude intervals; S322: Determine the order of magnitude interval of the current simulation time. When the simulation process moves from the previous order of magnitude interval to the next order of magnitude interval, adjust the candidate step size according to the preset interval step size mapping relationship so that the logarithm of the adjusted step size is linearly related to the simulation process, and obtain the final next time step size.

8. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 1, characterized in that, In step S2, the momentum and mass conservation equations for the fluid domain are adopted using the incompressible Navier-Stokes equations approximated by Boussinesq. The finite deformation equations for the solid domain are adopted using hyperelastic or viscoelastic-plastic constitutive equations described by the Piola-Kirchhoff stress tensor. In step S3, the mesh motion in the updated mesh and field quantities is controlled by a biharmonic equation, which is a mesh motion equation with the composite of the Laplace operator as the controlling operator.

9. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to claim 8, characterized in that, The fluid domain momentum and mass conservation equations, the solid domain finite deformation equations, and the global thermodynamic equations are constructed in a unified variational weak form under the reference configuration. The unified variational weak form uses fluid velocity, pressure, solid displacement, and temperature as unknowns and includes fluid momentum terms, fluid convection terms, fluid stress terms, pressure divergence terms, solid inertia terms, solid stress terms, thermal transient and convection terms, thermal diffusion terms, as well as buoyancy terms, volume force terms, and internal heat source terms. The velocity continuity and stress equilibrium conditions at the fluid-solid interface are satisfied in the unified variational weak form.

10. The dynamic terrain thermodynamic simulation method based on a fluid-structure interaction framework according to any one of claims 1-9, characterized in that, During the time progression in steps S1 to S3, a multi-mechanism numerical stabilization process is applied. The multi-mechanism numerical stabilization process includes: discretizing the time derivative term using a backward second-order difference scheme in step S2; performing pressure regularization when solving the fluid domain equation in step S2; and applying artificial viscosity based on entropy viscosity when solving the thermodynamic equation in step S2.

Images