Finite element discretization analysis method for long strip foundation pits based on spatiotemporal mechanical coupling

CN122572041APending Publication Date: 2026-08-14ANHUI MINGSHENG ELECTRIC POWER DESIGN CO LTD
View PDF 0 Cites 0 Cited by

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2026-05-26
Publication Date
2026-08-14

AI Technical Summary

Technical Problem

[0008]为了克服现有技术的上述缺陷,本发明的实施例提供基于时空力学耦合的长条形基坑有限元离散分析方法,将土体蠕变时间效应与基坑纵向空间约束效应在统一的三维有限元框架内强耦合求解,以解决现有方法长期变形预测失真的技术问题

Benefits of technology

1.本发明通过设定基坑纵向端部的法向位移约束条件以引入空间约束效应,同时引入土体蠕变参数以引入时间效应,并在计算各时间增量步的总应变增量、更新有效应力及蠕变内变量的过程中将两者纳入统一的求解流程,实现了三维空间约束与土体蠕变时间效应的强耦合求解,解决了现有技术中空间效应与时间效应割裂分析导致长期变形预测失真的技术问题。

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN122572041A_ABST
    Figure CN122572041A_ABST
Patent Text Reader

Abstract

This invention discloses a finite element discrete analysis method for elongated foundation pits based on spatiotemporal mechanical coupling, belonging to the field of numerical analysis technology in geotechnical engineering. The method includes the following steps: acquiring geometric data of the foundation pit, soil creep parameters, and construction step sequence time data; establishing normal displacement constraints at the longitudinal ends of the foundation pit; determining the excavation unloading range based on the construction step sequence time data and calculating the time-varying unloading vector; calculating the total strain increment for each time increment based on the constraints and the time-varying unloading vector, updating the corresponding effective stress and creep internal variables, and updating the consistent tangent stiffness matrix accordingly, iteratively solving the displacement and stress fields; and outputting the spatiotemporal distribution results of foundation pit deformation and internal forces of the support structure based on the displacement and stress fields for each time increment. This invention solves the technical problem of distortion in long-term deformation prediction caused by the separation of spatial and temporal effects in existing methods, significantly improving the accuracy of full-lifecycle deformation prediction for elongated foundation pits in soft soil areas.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This invention relates to the field of numerical analysis technology in geotechnical engineering, and more specifically, to a finite element discrete analysis method for long strip foundation pits based on spatiotemporal mechanical coupling. Background Technology

[0002] With the large-scale development of urban underground space, long and narrow foundation pit projects such as subway stations, underground commercial streets, and integrated utility tunnels are becoming increasingly common. A typical characteristic of these foundation pits is that their length is much greater than their width, with a length-to-width ratio often exceeding 5, and sometimes even reaching over 10. The excavation process involves soil unloading, leading to lateral displacement of the retaining structure, bottom heave, and settlement of the surrounding strata. Extensive engineering measurements and theoretical studies have shown that the deformation behavior of long and narrow foundation pits is significantly controlled by two effects: firstly, a spatial effect, where the constraint of the end walls at both ends of the pit causes deformation to be distributed longitudinally with a "larger deformation in the middle and smaller deformation at both ends" pattern; secondly, a temporal effect, where the creep and consolidation of the soft clay after unloading cause deformation to continue to increase even after excavation is completed.

[0003] To assess and predict the deformation and safety status of elongated foundation pits, the following numerical analysis methods are currently mainly used: The plane strain finite element method. This method simplifies the foundation pit into a typical cross-section for two-dimensional analysis, resulting in high computational efficiency and widespread application. However, it implicitly assumes an infinitely long longitudinal section of the foundation pit and completely ignores the constraints of the end walls and longitudinal structures. This leads to a systematic overestimation of the deformation and internal forces in the mid-span region, while the calculated values ​​near the ends are severely distorted, failing to reflect the true three-dimensional spatial effects.

[0004] The conventional three-dimensional finite element method (FEM) is used to construct a three-dimensional numerical model to account for the spatial geometry and constraints of the excavation pit, effectively capturing longitudinal deformation differences. However, existing conventional three-dimensional analyses often employ time-independent elastoplastic models for soil constitutive models, such as the Mohr-Coulomb model and hardened soil models. These models tend to stabilize after each excavation step, failing to describe the continuous deformation of soft soil caused by creep and consolidation during long-term static conditions. This severely limits the predictive ability of conventional three-dimensional methods for long-term post-construction deformation, underestimating the subsequent load increase on the support structure and the development of surface settlement, posing a threat to the long-term safety of the project.

[0005] Finite element methods considering time effects. To compensate for the above shortcomings, some researchers have introduced Biot's consolidation theory or empirical creep models. However, existing foundation pit analyses considering time effects are mostly limited to two-dimensional plane strain frameworks, or although they are three-dimensional models, they only consider pore pressure dissipation and ignore eccentric creep, or they adopt simple linear superposition or sequential weak coupling strategies in the coupling treatment of spatial and time effects. Since the creep rate of long strip foundation pits is highly dependent on the eccentric stress level determined by spatial constraints at that location, and creep development, in turn, changes the spatial stress distribution, there is a strong nonlinear interaction between the two. Weakly coupled or uncoupled methods cannot truly reproduce this physical process, resulting in significant deviations between the calculated results and the actual situation in both the time domain and spatial distribution.

[0006] In summary, current technologies lack a finite element discrete analysis method capable of fully simulating the longitudinal constraint effect of elongated foundation pits in three-dimensional space, while simultaneously coupling soil mass ratio-related mechanical behaviors (creep and / or consolidation) with three-dimensional stress redistribution. This gap leaves reliable theoretical tools and technical support unavailable for the full life-cycle deformation prediction and refined construction control of elongated foundation pits in soft soil areas.

[0007] Therefore, it is urgent to propose a finite element discrete analysis method for long strip foundation pits based on spatiotemporal mechanical coupling, so as to accurately solve the nonlinear coupling problem between three-dimensional spatial constraints and soil time effects, thereby improving the accuracy of foundation pit deformation and internal force prediction, and providing a basis for engineering design, construction optimization and long-term safety assessment. Summary of the Invention

[0008] To overcome the aforementioned deficiencies of the prior art, embodiments of the present invention provide a finite element discrete analysis method for long strip foundation pits based on spatiotemporal mechanical coupling. This method strongly couples the soil creep time effect with the longitudinal spatial constraint effect of the foundation pit within a unified three-dimensional finite element framework, thereby solving the technical problem of long-term deformation prediction distortion in existing methods.

[0009] To achieve the above objectives, the present invention provides the following technical solution: The finite element discrete analysis method for long strip foundation pits based on spatiotemporal mechanical coupling includes the following steps: acquiring the geometric data of the foundation pit, soil creep parameters, and construction step sequence time data; setting normal displacement constraints at the longitudinal ends of the foundation pit based on the geometric data; determining the excavation unloading range based on the construction step sequence time data and calculating the time-varying unloading vector; calculating the total strain increment for each time increment based on the constraints and the time-varying unloading vector, and updating the corresponding effective stress and creep internal variables; updating the consistent tangent stiffness matrix based on the effective stress and creep internal variables, and iteratively solving the displacement field and stress field for each time increment; and outputting the spatiotemporal distribution results of foundation pit deformation and internal forces of the support structure based on the displacement field and stress field for each time increment.

[0010] In a preferred embodiment, setting the normal displacement constraint conditions at the longitudinal ends of the foundation pit includes: extracting the set of mesh nodes located at the outer boundaries of both longitudinal ends of the foundation pit; and setting the normal displacement degree of freedom of each node in the set of mesh nodes to zero.

[0011] In a preferred embodiment, determining the excavation unloading range based on the construction step sequence time data includes: parsing the construction step sequence time data and extracting the longitudinal segment boundary coordinates and excavation depth elevation corresponding to the current time increment step; traversing the finite element discrete mesh pre-established from the foundation pit geometry data and extracting the elements whose spatial positions fall within the range of the longitudinal segment boundary coordinates and excavation depth elevation as the target element set; marking the elements in the target element set as removed, and using the spatial region corresponding to the target element set as the excavation unloading range of the current time increment step.

[0012] In a preferred embodiment, the calculation of the time-varying unloading vector includes: extracting the element stress of the target element set before it is marked as removed; calculating the initial equivalent release node force generated by the target element set on the remaining effective grid nodes; extracting the creep intrinsic variables of the remaining effective grid elements at the excavation boundary; calculating the creep state evolution increment of the current time increment step; using the creep state evolution increment as an endogenous driving parameter to control the release process of the initial equivalent release node force; and autonomously allocating the load release amount of the current time increment step by the internal rheological mechanism of the soil to generate the corresponding time-varying unloading vector.

[0013] In a preferred embodiment, the step of calculating the total strain increment for each time increment step based on constraints and time-varying unloading vectors includes: constructing the system equilibrium equation corresponding to the current time increment step based on the time-varying unloading vector and normal displacement constraints; using the displacement field obtained from the previous time increment step as the initial displacement guess, and solving the system equilibrium equation using Newton-Raphson iteration to obtain the nodal displacement increment for the current time increment step; and converting the nodal displacement increment into the total strain increment for the current time increment step through geometric equations.

[0014] In a preferred embodiment, updating the corresponding effective stress and creep internal variables includes: calculating the test effective stress based on the total strain increment according to the elastic constitutive relation; calculating the corresponding creep strain increment based on the test effective stress and soil creep parameters; correcting the test effective stress based on the creep strain increment to obtain the effective stress corresponding to the current time increment step; and updating the creep internal variables corresponding to the current time increment step based on the creep strain increment.

[0015] In a preferred embodiment, updating the consistent tangent stiffness matrix based on effective stress and creep intrinsic variables includes: extracting the additional stiffness contribution generated by the nonlinear evolution of the time-varying unloading vector with creep intrinsic variables, and constructing an additional boundary coupling stiffness term; constructing a basic consistent tangent modulus matrix based on the effective stress and creep intrinsic variables of the current time increment step, combined with the creep evolution rule; and fusing and assembling the additional boundary coupling stiffness term with the basic consistent tangent modulus matrix to update the consistent tangent stiffness matrix of the current time increment step.

[0016] In a preferred embodiment, the construction of the additional boundary coupling stiffness term includes: calculating the first partial derivative of the time-varying unloading vector with respect to the creep internal variable; calculating the second partial derivative of the creep internal variable with respect to the total strain increment; combining the first partial derivative, the second partial derivative, and the finite element geometric equation, analytically obtaining the non-zero gradient matrix of the external load vector with respect to the nodal displacement increment using the chain rule, and using the non-zero gradient matrix as the additional boundary coupling stiffness term.

[0017] In a preferred embodiment, the iterative solution of the displacement field and stress field for each time increment step includes: updating the time-varying unloading vector according to the creep internal variable corresponding to the displacement field of the current iteration step within each nonlinear iteration step; calculating the difference between the updated time-varying unloading vector and the internal nodal force vector obtained from the effective stress to obtain the global unbalanced residual force vector; constructing a linearized incremental equation from the uniform tangent stiffness matrix and the global unbalanced residual force vector, solving the iterative displacement correction vector and updating the nodal displacements; iterating until the global unbalanced residual force vector satisfies the preset convergence tolerance, and outputting the displacement field and stress field of the current time increment step.

[0018] In a preferred embodiment, the step of outputting the spatiotemporal distribution results of the excavation pit deformation and the internal forces of the support structure based on the displacement field and stress field of each time increment step includes: extracting the displacement sequence of the target observation node from the displacement field to obtain the excavation pit deformation of the current time increment step; extracting the integral point stress of the support structure unit from the stress field, and obtaining the internal forces of the support structure of the current time increment step through cross-sectional integration calculation; and assembling the excavation pit deformation and the internal forces of the support structure of each time increment step sequentially according to the calculation time axis to generate and output the distribution results characterizing the spatiotemporal evolution characteristics of the entire excavation process.

[0019] The technical effects and advantages of this invention's finite element discretization analysis method for elongated foundation pits based on spatiotemporal mechanical coupling are as follows: 1. This invention introduces a spatial constraint effect by setting a normal displacement constraint condition at the longitudinal end of the foundation pit, and simultaneously introduces soil creep parameters to introduce a time effect. In the process of calculating the total strain increment, updating the effective stress and creep internal variables at each time increment step, both are incorporated into a unified solution process, realizing a strong coupling solution of three-dimensional spatial constraints and soil creep time effect. This solves the technical problem of distortion in long-term deformation prediction caused by the separate analysis of spatial and time effects in the prior art.

[0020] 2. This invention improves the accuracy of predicting the full life cycle deformation of long strip foundation pits in soft soil areas by updating the consistent tangential stiffness matrix based on the effective stress and creep internal variables in each time increment step and iteratively solving the displacement field and stress field, so that the evolution direction of the system stiffness matrix is ​​consistent with the evolution direction of the soil creep state. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the finite element discretization analysis method for long strip foundation pits based on spatiotemporal mechanical coupling of the present invention. Figure 2 This is a structural schematic diagram illustrating the three-dimensional finite element mesh discretization and longitudinal end normal constraint in an embodiment of the present invention. Figure 3 This is a curve showing the distribution of lateral displacement of the retaining structure along the longitudinal direction of the foundation pit, output by an embodiment of the present invention. Figure 4 This is a curve showing the evolution of the support axial force over time, output by an embodiment of the present invention. Figure 5 This is a comparison curve of the nonlinear iterative convergence process of the method of this invention and the existing method in a typical excavation step; Figure 6 This is a bar chart comparing the prediction accuracy of key deformation indicators of the method of this invention with existing methods and field measured values. Detailed Implementation

[0022] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0023] Example 1, Figure 1 This invention presents a finite element discretization analysis method for elongated foundation pits based on spatiotemporal mechanical coupling, comprising the following steps: S1, acquire the geometric data of the foundation pit, soil creep parameters and construction sequence time data; In this embodiment, acquiring the geometric data of the foundation pit, soil creep parameters, and construction sequence time data specifically includes: Obtain the geometric data of the foundation pit: Extract the three-dimensional spatial dimensions and geological structure information of the foundation pit, including but not limited to the total length, total width, maximum design excavation depth, and spatial distribution elevation and thickness of each soil layer; simultaneously obtain the spatial layout parameters of the support structure, such as the cross-sectional thickness and embedment depth of the diaphragm wall or piles, and the planar spacing and vertical elevation of each internal support. Generate an initial finite element discrete mesh model based on the geometric data, and establish the node topology relationship between soil solid elements and support structure elements.

[0024] Obtaining Soil Creep Parameters: A nonlinear constitutive model and corresponding physical control parameters are determined to characterize the rheological properties of the soil surrounding the excavation pit. As a preferred embodiment, the constitutive model may be an empirical time-hardening creep model, a Burgers viscoelastic rheological model, or a unified hardening (UH) creep model based on critical states. The corresponding soil creep parameters include, but are not limited to: initial elastic modulus, Poisson's ratio, material viscosity coefficient, reference creep time, and secondary consolidation coefficient. These soil creep parameters can be obtained through inversion from indoor triaxial rheological tests or actual field monitoring data to ensure the physical authenticity of the spatiotemporal analysis.

[0025] Obtain construction step sequence time data: Analyze the construction organization design scheme and extract work condition division data with absolute or relative time axis attributes. Specifically, this includes: multiple excavation stage sequences, the start time, end time, and current excavation depth and elevation corresponding to each excavation stage; in particular, it is also necessary to extract the duration of excavation operations (to characterize the unloading rate) and the downtime exposure or subbase pouring waiting time between two adjacent excavation stages (to characterize the pure creep evolution duration of the soil under constant external force boundaries). The construction step sequence time data serves as the control basis for subsequent finite element time increment step division, providing time integration boundary conditions for triggering the endogenous creep driving mechanism.

[0026] S2, set the normal displacement constraint conditions at the longitudinal end of the foundation pit based on geometric data; In this embodiment, the setting of normal displacement constraints at the longitudinal end of the foundation pit based on geometric data specifically includes: S201. Based on the geometric data obtained in step S1, such as the total length, total width, excavation depth, and diaphragm wall thickness of the foundation pit, a finite element discrete mesh model is generated, and a system global coordinate system is established in computer memory. In this embodiment, the longitudinal direction of the foundation pit is the X-axis, the transverse direction of the short side of the foundation pit is the Y-axis, and the depth direction of the foundation pit is the Z-axis. According to the start and end points of the design length of the foundation pit in the geometric data, the computational domain geometric boundary surfaces at both ends of the longitudinal direction of the foundation pit are determined. The outer boundaries at both ends of the longitudinal direction include the end nodes of the foundation pit retaining structure (such as a diaphragm wall) in the longitudinal direction, and the longitudinal boundary nodes of the soil extending outward from both ends of the foundation pit. Specifically, in the three-dimensional finite element mesh, all mesh nodes whose coordinate values ​​are equal to or exceed the boundary surfaces at the endpoints of the design length of the foundation pit are selected along the longitudinal direction of the foundation pit to form the set of mesh nodes. For long strip foundation pits, since the length is much greater than the width and depth, the end constraints are only applied to the two ends of the longitudinal direction, while the transverse boundary of the foundation pit remains free or adopts other corresponding constraint conditions.

[0027] S202, for each node in the extracted mesh node set, in the system's global coordinate system, its displacement degree of freedom along the longitudinal direction (X-axis) of the foundation pit is set to zero, that is:

[0028] in, Indicates the node along the longitudinal direction of the foundation pit (i.e. The normal translational displacement component (in the axial direction).

[0029] This operation is achieved by applying displacement constraints numerically using the finite element method. Specifically, this corresponds to the assembly and solution stage of the discrete control equations of the system, where the diagonal elements of the system's consistent tangent stiffness matrix corresponding to the degrees of freedom of the corresponding nodes are multiplied by a large number (i.e., the penalty function method), and the corresponding components of the global unbalanced residual force vector on the right-hand side of the system are set to zero; or the matrix is ​​directly reduced in order using the row and column division method.

[0030] The significance of employing the aforementioned normal displacement constraint condition in this invention lies in the fact that for long, narrow foundation pits (such as subway stations and elongated integrated utility tunnels), there are usually unexcavated original soil or rigid end well retaining structures at both ends. These longitudinal ends exert a strong normal deformation constraint on the soil inside the foundation pit. By forcibly constraining the normal displacement of the mesh nodes at the outer boundaries of both ends to zero, this "spatial constraint effect" can be realistically and accurately introduced into the three-dimensional finite element calculation. This spatial constraint will directly affect the deformation transmission and stress redistribution in the middle of the foundation pit through the subsequent assembly of the stiffness matrix, thereby avoiding the deficiency of the traditional plane strain assumption (2D model) which cannot consider the end support effect and thus leads to conservative calculation results.

[0031] Figure 2This diagram illustrates the three-dimensional finite element mesh discretization and longitudinal end normal constraint of the elongated foundation pit in this embodiment. In the diagram, the X-axis represents the longitudinal direction of the foundation pit, the Y-axis represents the transverse direction, and the Z-axis represents the depth direction. The blue areas at both ends of the foundation pit represent the end boundaries where normal displacement constraints are applied. These constraints cover the end nodes of the foundation pit retaining structure and the far-end boundary nodes of the soil in the longitudinally extended area. Red arrows indicate the location and direction of the applied normal displacement constraints. The three-dimensional finite element mesh domain (represented by mesh filling in the diagram) covers the excavation area of ​​the foundation pit and its surrounding influence zone, providing a spatial discretization basis for the numerical calculations in subsequent steps S3 to S6.

[0032] S3, determine the excavation unloading range based on the construction step sequence time data, and calculate the time-varying unloading vector; In this embodiment, determining the excavation unloading range based on construction step sequence time data includes: S301, based on the construction step sequence time data obtained in step S1, analyze the longitudinal segmentation range and vertical layering depth corresponding to each excavation step recorded therein. Let the excavation step sequence number corresponding to the current time increment step be... Extracting the first step from the construction sequence time data Vertical segment boundary coordinates of the step and excavation depth and elevation ,in and These are the start and end coordinates of the current excavation step along the longitudinal direction of the foundation pit. This is the starting elevation of the current excavation face (usually the bottom elevation of the previous excavation step). This represents the design bottom elevation of the current excavation step. The coordinates and elevations mentioned above are derived from the spatial range information recorded in the construction step sequence time data in step S1.

[0033] S302, in the finite element discrete mesh pre-established from the foundation pit geometry data in step S1, traverse all soil solid elements. For each element, extract its centroid space coordinates. Determine whether the coordinates simultaneously satisfy the following two conditions: Longitudinal range conditions:

[0034] Depth range conditions:

[0035] Extract all units that simultaneously meet the above two conditions to form the target unit set for the current time increment step. .

[0036] In the spatial range determination, the coordinates of the longitudinal segment boundary are... and The excavation depth and elevation are directly given by the construction sequence time data. and The maximum design excavation depth of the foundation pit was obtained by combining the construction sequence time data with S1. And the layered excavation scheme is determined. Therefore, the extraction of the target unit set is the result of the combined effect of construction step time data and foundation pit geometric data - that is, time data determines "when to excavate", and geometric data determines "where to excavate".

[0037] S303, set the target unit All elements are marked as removed. In the underlying numerical implementation of the finite element method, this operation is performed by multiplying the stiffness matrix of the corresponding element by a very small attenuation factor (e.g., ...). This is accomplished by setting its mass matrix and load vector to zero, ensuring that the removed elements no longer contribute to the stiffness and load of the remaining finite element system. The spatial region corresponding to the target element set is used as the excavation unloading range for the current time increment step.

[0038] S304, in assembling the target unit Before marking a element as removed, the element stress of each element under the current stress state is extracted. The element stress is the effective stress field after initialization in step S1 or convergence of the previous time increment step. Based on the finite element virtual work principle, the initial equivalent release nodal force generated by the target element set on the remaining effective mesh nodes is calculated. :

[0039] in, The geometric strain matrix, For unit volume, For unit The initial effective stress tensor before removal This indicates that the summation is performed by iterating through all cells in the target cell set.

[0040] S305, extract the creep intrinsic variables of the remaining effective mesh elements at the excavation boundary. The excavation boundary refers to the target element set. The interface between the remaining effective mesh elements and the interface. For each remaining effective mesh element located on this interface, extract its creep intrinsic variable. Based on the extracted creep intrinsic variables, the creep state evolution increment at the current time increment step is calculated. :

[0041] in, This represents the real-time creep internal variable within each nonlinear iteration step at the current time increment step (this value is dynamically updated with the iterative displacement field). These are the historical creep internal variables after the previous time increment step converged.

[0042] Furthermore, the creep state evolution increment As an endogenous driving parameter, it controls the initial equivalent release node force. The release process is autonomously allocated by the soil's internal rheological mechanism to determine the load release amount for the current time increment step, generating a corresponding time-varying unloading vector. :

[0043] in, This is a release process control function with the creep state evolution increment as the independent variable. In this embodiment, The preferred form is the normalized form:

[0044] In the formula, The norm of the creep state evolution increment; This is a characteristic creep reference value (determined by the scalar value of the creep intrinsic variable corresponding to the reference creep time in the soil creep parameters obtained in step S1). This function ensures that: when the creep rate of the soil at the excavation boundary is large... When the creep is large, the load release is accelerated; when the creep tends to stabilize, The load release is slowed down accordingly, thus achieving adaptive synchronization between the load release process and the soil creep evolution.

[0045] This invention increments the creep state evolution. As an endogenous driving parameter, it replaces the pre-defined explicit time decay function (such as the exponential decay function or piecewise linear decay function) in traditional methods. In traditional methods, the load release process is controlled by an externally given time function, which is independent of the actual creep state of the soil itself. This leads to the dual inclusion of external load decay and internal creep relaxation in the time integral, causing serious distortion in deformation calculations. This invention transfers the control of the unloading process from the external time function to the creep state evolution increment of the soil itself, so that the driving force of load release originates from the internal rheological mechanism of the soil rather than an externally imposed mathematical setting. This completely eliminates this inherent technical bias of "spatiotemporal disconnect and dual calculation" from the physical underlying mechanism.

[0046] S4, based on the constraints and time-varying unloading vector, calculates the total strain increment for each time increment step and updates the corresponding effective stress and creep internal variables; In this embodiment, the step of calculating the total strain increment for each time increment step based on constraints and time-varying unloading vectors, and updating the corresponding effective stress and creep internal variables, includes the following steps: S401, based on the longitudinal end normal displacement constraint condition of the foundation pit set in step S2, and the time-varying unloading vector obtained in step S3. Construct the equilibrium equations of the finite element system:

[0047] in, The overall stiffness matrix of the system. Let be the nodal displacement increment vector to be solved. The nodal force vectors within the system are obtained by integrating the effective stresses at the integration points of each element corresponding to the current displacement field using the geometric strain matrix.

[0048] In the formula, This indicates that the summation is performed by iterating through all valid grid cells. The strain-displacement matrix, composed of the partial derivatives of the shape functions of the finite element with respect to spatial coordinates, is a standard operator in the finite element method, used to represent discrete nodal displacement increments. Converted to continuous total strain increments at integration points , The current effective stress tensor of the element. This is the integral domain of the element. The physical meaning of the internal nodal force vector is the equivalent resistance generated by the internal stress state of the soil on the node, and the difference between it and the external load vector constitutes the unbalanced residual force of the system.

[0049] The normal displacement constraint is applied to the system equilibrium equations using either the row-and-column method or the multiplication method: the diagonal elements of the stiffness matrix corresponding to the degree of freedom of the constrained node are multiplied by a large number (e.g., 10151015), and the corresponding components of the right-hand load vector are set to zero, ensuring that the solution automatically satisfies the constraint condition that the longitudinal end normal displacement is zero. The time-varying unloading vector... As an external load term, it directly enters the right-hand side of the system equilibrium equations.

[0050] S402, increment the previous time step (the nth time step) Step 1) Solve for the displacement field As the initial displacement guess for the current time increment step, the Newton-Raphson iterative method is used to solve the system equilibrium equations. In each iteration step... Within, perform the following sub-steps: S402-1, Guess the displacement value based on the current iteration step. The effective stress at each element integration point is updated using geometric equations and constitutive relations, and according to the above... The calculation formula assembles the internal node force vector of the current iteration step. .

[0051] S402-2, Calculate the global unbalanced residual force vector of the system :

[0052] S402-3, Judgment Does the norm satisfy the preset convergence tolerance? In this embodiment, the convergence tolerance... Pick to The specific value is determined based on the engineering precision requirements. If (in If the initial unbalanced residual force vector is given, then the iteration converges, and proceed to step S403.

[0053] S402-4, if the convergence condition is not met, then the current uniform tangent stiffness matrix is ​​used. Construct a linearized incremental equation and solve for the iterative displacement correction vector. :

[0054] Update node displacement:

[0055] Return to substep S402-1 to proceed to the next iteration.

[0056] After iterative convergence, the node displacement increment at the current time increment step is obtained. :

[0057] in, This represents the node displacement vector at the current time increment step when the iteration converges. This is the node displacement vector at the convergence of the previous time increment step.

[0058] S403, based on the finite element geometric equations, obtains the nodal displacement increments. Transformed into the total strain increment at the integration points of each grid element :

[0059] S404, based on the calculated total strain increment First, assume that the strain increment is entirely elastic strain, and calculate the effective stress of the test according to the elastic constitutive relation. :

[0060] in, The effective stress tensor that converged in the previous time increment step. Here is the elastic stiffness matrix. This indicates the double dot product operation.

[0061] S405, based on the obtained experimental effective stress Based on the soil creep parameters obtained in step S1, the creep strain increment corresponding to the current time increment step is calculated. In this embodiment, the creep model adopted is the Singh-Mitchell creep model, and the creep strain increment is determined by the following formula:

[0062] in, The time step is the increment of the current time. The creep coefficient is 1. The stress level index, For the deviatoric stress ratio level, The time decay exponent, Let be the partial derivative of the creep potential function with respect to the effective stress. For reference time, This is the cumulative time since the start of the current excavation step.

[0063] S406, using the return mapping algorithm, the creep strain increment obtained in step (5) is used to correct the tested effective stress, and the true effective stress corresponding to the current time increment step is obtained. :

[0064] S407, based on the obtained creep strain increment Update the creep internal variables corresponding to the current time increment step. :

[0065] in, These are the creep internal variables from the previous time increment step. This represents the increment of the creep internal variable caused by the creep strain increment. In this embodiment, the creep internal variable is taken as the equivalent creep strain. :

[0066] Updated creep internal variables The creep state evolution increments stored at the integration points of each grid cell will be used for the next time increment step. The calculation of the uniform tangent stiffness matrix is ​​performed, and the update of the uniform tangent stiffness matrix is ​​performed in subsequent steps.

[0067] This embodiment constructs a rigorous bidirectional data transmission stream from "global displacement field solution" to "elastic probing and creep pullback at local integration points," and then back to "internal nodal forces of the global system" through inverse integration mapping. By introducing strain decomposition and stress pullback mechanisms at the Gaussian integration point level, not only are time-dependent inelastic rheological components of the material accurately isolated, ensuring the physical authenticity of constitutive updates under complex stress paths, but more importantly, the true effective stress calculated in this step is... It became the sole source of updating the internal node forces (internal forces) of the system, while the creep internal variables were updated synchronously. This becomes the intrinsic driver for refreshing the time-varying unloading vector (external force) in step S3. This makes both the "internal force" and the "external force" essentially strongly coupled with the current deformation state of the system within the same iteration step, laying an absolutely rigorous legal and mathematical foundation for constructing the "boundary equivalent coupling stiffness" through the chain partial derivative rule and achieving the second-order accelerated convergence of the strongly nonlinear equation in subsequent steps.

[0068] S5, update the consistent tangent stiffness matrix based on the effective stress and creep internal variables, and iteratively solve the displacement field and stress field for each time increment step; In this embodiment, the step of updating the consistent tangent stiffness matrix based on the effective stress and creep internal variables, and iteratively solving the displacement field and stress field at each time increment step, includes the following steps: S501, Effective stress based on the current time increment step output in step S4. and creep internal variables Based on the creep evolution law, the first-order partial derivatives of the creep constitutive relation with respect to the total strain are calculated at the integration points of each grid element. These partial derivatives reflect the real-time influence of the soil creep state on the material's tangential stiffness. The first-order partial derivatives are then compared with the elastic stiffness matrix. By combining the results, we can construct the fundamental consistent tangent modulus matrix at each integration point. .

[0069] S502, Extract the time-varying unloading vector generated in step S3. Relative to creep internal variables The partial derivative, i.e., the first partial derivative. This partial derivative characterizes the sensitivity of the soil to the creep state caused by the external load, and its value is determined by the release process control function defined in step S3. The specific form is determined. Simultaneously, creep intrinsic variables are extracted. Relative to the total strain increment The partial derivative, i.e., the second partial derivative. The partial derivative is obtained by uniformly linearizing the return mapping algorithm of elastic prediction-creep correction in step S4, and reflects the local constitutive relationship of creep internal variables as strain evolves.

[0070] The first partial derivative, the second partial derivative, and the finite element geometric equation ( By combining these methods, the non-zero gradient matrix of the external load vector with respect to the nodal displacements can be analytically obtained using the chain rule:

[0071] This non-zero gradient matrix represents the additional boundary coupling stiffness term. It is important to note that in conventional finite element analysis, external loads are typically treated as known quantities independent of nodal displacements, and their partial derivatives with respect to displacement are always zero. However, in this invention, due to the time-varying unloading vector... Defined as creep internal variable in step S3 Since the creep internal variable is itself a function of the displacement field, the above partial derivatives are no longer zero. The additional boundary coupling stiffness term is the quantitative expression of this non-zero coupling relationship. It precisely transforms the physical coupling between the external boundary unloading evolution and the internal soil creep into a calculable additional stiffness contribution at the numerical implementation level.

[0072] S503, the basic consistent tangent modulus matrix obtained in step S501. Following the standard finite element assembly process (i.e., combining the strain-displacement matrix) Perform spatial integration Generate the basic global stiffness matrix and combine it with the additional boundary coupling stiffness term obtained in step S502. Perform system-level fusion to obtain the complete consistent tangent stiffness matrix for the current time increment step. :

[0073] in, The fundamental global stiffness matrix is ​​obtained by assembling the fundamental consistent tangent modulus matrix using standard finite element methods. This is the contribution term of the additional boundary coupling stiffness term after assembly with global degrees of freedom. In the Newton-Raphson iteration, this complete and consistent tangent stiffness matrix simultaneously absorbs the dual unbalanced residual forces generated by the internal soil creep relaxation and the endogenous unloading evolution of the external boundary, providing a tangent direction consistent with the actual physical evolution for the iterative solution and ensuring the quadratic convergence of the algorithm.

[0074] S504, within each nonlinear iteration step of the current time increment step, execute the following strongly coupled iterative process: S504-1, based on the displacement field of the current iteration step Corresponding creep internal variables Recalculate the time-varying unloading vector according to the method in step S3. This step breaks the implicit assumption in the conventional Newton-Raphson iteration that the external load remains constant within the incremental step, allowing the boundary load to be updated synchronously with the soil creep state.

[0075] S504-2, the effective stress field corresponding to the current displacement field. via strain-displacement matrix Integral assembly of internal nodal force vectors The difference between the updated time-varying unloading vector and the internal nodal force vector is calculated to obtain the global unbalanced residual force vector of the system.

[0076] This residual force vector contains the dual contributions of internal stress imbalance force and external load endogenous evolution imbalance force.

[0077] S504-3, based on current effective stress and creep internal variables Update the consistent tangent stiffness matrix according to steps S501 to S503. Following step S402-4, a linearized incremental equation is constructed using the uniform tangent stiffness matrix and the global unbalanced residual force vector, and the iterative displacement correction vector is solved. and update node displacements. .

[0078] S504-4, Determine the global unbalanced residual force vector Does the norm satisfy the preset convergence tolerance? In this embodiment, the convergence tolerance... Pick to If satisfied, output the displacement field at the current time increment step. With stress field If the condition is not met, return to step S504-1 to perform the calculation for the next iteration until convergence.

[0079] It should be noted that this embodiment, starting from the underlying partial derivative logic of the Newton-Raphson iteration, discovers a non-zero term that has long been neglected in traditional finite element analysis—the partial derivative of the external load vector with respect to the nodal displacement. In conventional analysis, this term is always zero because the external load is independent of structural deformation. However, under the endogenous driving mechanism established in step S3 of this invention, the time-varying unloading vector is defined as a function of the creep internal variable, and the creep internal variable is associated with the total strain (and thus with the nodal displacement) through the return mapping algorithm in step S4, making the above partial derivative no longer zero. This step uses the chain rule to analytically express this non-zero term as an additional boundary coupling stiffness term and achieves system-level fusion with the consistent tangent modulus of the foundation, so that the consistent tangent stiffness matrix fully reflects the dual contribution of "material constitutive stiffness + boundary coupling stiffness". In each nonlinear iteration step, the external load is dynamically refreshed with the update of the creep internal variable and is performed synchronously with the residual calculation of the internal nodal force vector, realizing the coordinated convergence of "internal creep relaxation" and "external load evolution". This technical solution mathematically guarantees the quadratic convergence rate of the strongly nonlinear equations and physically eliminates the deformation prediction distortion caused by the spatiotemporal effect fragmentation in the traditional sequential weak coupling method. It represents a fundamental reconstruction of the finite element analysis method for foundation pits from the underlying numerical mechanism level.

[0080] S6 outputs the spatiotemporal distribution results of foundation pit deformation and internal forces of the support structure based on the displacement field and stress field of each time increment step.

[0081] In this embodiment, the step of outputting the spatiotemporal distribution results of the foundation pit deformation and the internal forces of the support structure based on the displacement field and stress field at each time increment step specifically includes the following steps: S601, Displacement field from all time increment steps output in step S5 ( , Within the total time increment steps, the displacement sequence of the target observation nodes is extracted. The target observation nodes include: nodes of the retaining structure (such as a diaphragm wall) at each longitudinal section, nodes arranged longitudinally on the bottom surface of the pit, and nodes arranged laterally and longitudinally on the ground surface around the pit.

[0082] For the lateral displacement of the retaining structure, the horizontal displacement components of each observation node along the depth direction (Z-axis direction) of the foundation pit are extracted at each time increment step, generating the distribution curves of the lateral displacement of the retaining structure along the depth direction and longitudinal direction of the foundation pit. This curve intuitively reflects the distribution pattern of the lateral displacement of the retaining structure along the depth at different locations in the longitudinal direction of the foundation pit (mid-span, quarter-span, end) and its evolution over time.

[0083] For the pit bottom heave, the vertical displacement components of each observation node on the pit bottom surface are extracted at each time increment step to generate the distribution curve of the pit bottom heave along the longitudinal direction of the pit. This curve reflects the spatial non-uniform distribution characteristics of the pit bottom rebound deformation under the combined action of excavation unloading and soil creep.

[0084] For surface settlement, the vertical displacement components of the surface observation nodes around the foundation pit are extracted at each time increment step to generate the distribution curves of the surface settlement trough along the transverse and longitudinal directions.

[0085] S602, Stress field from all time increment steps output in step S5 In the process, the stress components of the support structure unit (including diaphragm wall unit and support unit) at each integration point are extracted.

[0086] For the axial force of the supports, the axial stress components of each support element at each time increment step are extracted, and the axial force value of each support is calculated by cross-sectional integration:

[0087] in, For the cross-sectional area of ​​the support, To support axial stress, the evolution curves of the axial force of each support along the longitudinal direction of the foundation pit and over time are generated.

[0088] For the bending moment of the enclosure structure, the bending normal stress components of the diaphragm wall element at each integration point are extracted, and the bending moment value at each section is obtained by calculating the stress integral of the section:

[0089] in, This refers to the cross-sectional area per unit width of the diaphragm wall. For bending normal stress, The distance from the integration point to the neutral axis of the cross section is given. The curves of the bending moment of the retaining structure along the depth of the excavation pit, longitudinally, and over time are generated.

[0090] S603, the excavation deformation data (lateral displacement of retaining structure, heave of pit bottom, and ground settlement) and the internal force data of the support structure (axial force of support and bending moment of retaining structure) obtained in each time increment step of steps S601 and S602 are sequentially assembled according to the calculation time axis. The sequential assembly refers to: assembling the data of each time increment step... ( The calculation results are arranged in chronological order to form a continuous spatiotemporal evolution data sequence covering the entire time range from the first excavation step to the last settling period.

[0091] Based on the assembled data sequence, the distribution results characterizing the spatiotemporal evolution of the entire excavation process are generated and output, including but not limited to: spatiotemporal distribution cloud maps of the lateral displacement of the retaining structure along the depth and longitudinal direction, time history curve clusters of pit bottom uplift along the longitudinal direction, evolution curve clusters of axial forces of each support over time, envelope diagram of the bending moment of the retaining structure along the depth and its evolution over time.

[0092] This embodiment transforms the complete calculation process of the strongly coupled spatiotemporal mechanical finite element analysis constructed in steps S1 to S5 into spatiotemporal distribution results of foundation pit deformation and internal forces of the support structure, which can be directly applied in engineering, through a three-step progressive data processing chain of "displacement field extraction → stress field extraction → time axis sequence assembly". Unlike conventional finite element analysis, which only outputs the deformation or internal force envelope value of a single time step, the spatiotemporal distribution curves and evolution curves output in this step completely record the evolution trajectory of the foundation pit's mechanical behavior from the beginning of excavation to long-term static placement. In particular, since the previous steps have strongly coupled the spatial constraint effect and the soil creep time effect within a unified framework, the results output in this step can simultaneously reflect the superposition effect of "non-uniform deformation space caused by end constraints" and "continuous growth of deformation time caused by soil creep", as well as the nonlinear interaction between the two formed by endogenous driving and equivalent boundary coupling stiffness. This provides accurate numerical basis for the full life cycle deformation prediction, long-term safety assessment of support structures, and optimization of construction procedures for long strip foundation pits in soft soil areas, which cannot be provided by existing methods.

[0093] Figure 3 The diagram shows the distribution curves of the lateral displacement of the retaining structure along the longitudinal direction of the foundation pit, as output in this embodiment. The two curves correspond to the lateral displacement distribution immediately after excavation and after 30 days of settling, respectively. The "groove-shaped" distribution of lateral displacement, which is largest in the mid-span region and decreases towards both ends, reflects the spatial constraint effect generated by the longitudinal end normal displacement constraint in step S2. The overall increase in lateral displacement after 30 days of settling reflects the time effect generated by soil creep relaxation in step S4. The increase in lateral displacement in the mid-span region is significantly greater than that in the end region, demonstrating the inhibition of creep development by spatial constraints—the spatial effect and the time effect exhibit a strong coupling characteristic here.

[0094] Figure 4 The curve showing the evolution of the support axial force over time is illustrated in this embodiment. Traditional methods (gray dashed line) ignore the creep time effect, resulting in a constant support axial force after excavation. The method of this invention (red solid line) accurately captures the physical process of soil creep relaxation leading to continuous lateral displacement of the retaining structure, passive compression of the support, and a continuous increase in axial force over time with a gradually decreasing rate of increase, consistent with the measured patterns in foundation pit engineering in soft soil areas.

[0095] Example 2: To further verify the effectiveness and significant technological advancement of the method of the present invention, the following comparative analysis results between the method of the present invention and existing traditional algorithms are presented, combining numerical convergence experiments and measured data from a long strip foundation pit project: 1. Verification of algorithm's underlying convergence For the additional boundary coupling stiffness term constructed in step S5, this embodiment extracts the nonlinear iterative process data of a typical strong creep-unloading coupled excavation step. Table 1 shows a comparison of the convergence process of the algorithm of this invention and the traditional sequential weak coupling algorithm (i.e., assuming that the partial derivative of the external load with respect to displacement is zero).

[0096] Table 1. Comparison of the convergence process of the algorithm of this invention and the traditional algorithm in a typical excavation step.

[0097] As shown in Table 1, when facing strong spatiotemporal coupling conditions, due to the introduction of an additional boundary coupling stiffness term into the system equilibrium equation in this invention, the residual force norm of the system exhibits a quadratic exponential decrease in the fourth iteration step and satisfies the convergence tolerance ( Traditional methods, by neglecting the non-zero partial derivatives of external load evolution with respect to nodal displacements, cause the iteration direction to deviate from the true tangent, eventually diverging after several oscillations. This rigorously demonstrates, from a numerical experimental perspective, the significant advancements of the algorithm in computational stability and convergence rate of this invention.

[0098] Figure 5 The convergence process data in Table 1 are presented as curves. The horizontal axis represents the Newton-Raphson iteration step sequence, and the vertical axis represents the norm of the global unbalanced residual force vector. In the figure, the blue solid line represents the convergence curve of the method of this invention, the red dashed line represents the convergence curve of the traditional method, and the gray dashed line marks the position of the preset convergence tolerance. Figure 5 As can be seen, the method of this invention satisfies the convergence tolerance in the fourth iteration, and the residual force norm decreases rapidly at a quadratic rate; the traditional method oscillates and gradually diverges after the third step. The root cause of this significant difference is that this invention constructs an equivalent coupling stiffness at the boundary and incorporates a consistent tangent stiffness matrix, ensuring that the iteration direction is consistent with the actual physical tangent; the traditional method, by assuming that the partial derivative of the external load with respect to displacement is zero, lacks the guidance of this coupling stiffness term, causing the iteration to deviate from the correct direction and diverge.

[0099] 2. Verification of Engineering Prediction Accuracy Based on the spatiotemporal distribution results output in step S6, this embodiment extracts the key deformation indicators of a long strip foundation pit (length-to-width ratio greater than 5) in a soft soil area after excavation to the base and 60 days of cessation of work. The results are compared with the field measured data and the prediction results of traditional numerical methods, as shown in Table 2.

[0100] Table 2 Comparison of Prediction Accuracy of Key Deformation Indicators during the Stagnation Period (60 days) of Long Strip Foundation Pit

[0101] As shown in Table 2, the traditional plane strain method suffers from extremely high end displacement error (up to +282.9%) due to neglecting end spatial constraints; while the traditional three-dimensional finite element method considers spatial effects, it underestimates the overall long-term deformation prediction by approximately 25% to 30% because it ignores the creep time effect. In contrast, the spatiotemporal distribution results output by the method of this invention are in high agreement with the measured values ​​on site (the relative errors of all key deformation indicators are strictly controlled within 6%), which precisely demonstrates the extremely high accuracy and engineering practical value of this invention in predicting the deformation of long strip foundation pits in soft soil throughout their entire life cycle. Figure 6 A bar chart comparing the prediction accuracy of key deformation indicators of the method of this invention with existing methods and field measured values ​​is provided.

[0102] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0103] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented using software, the above embodiments can be implemented, in whole or in part, in the form of a computer program product.

[0104] Those skilled in the art will recognize that the modules and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented in electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this application.

[0105] In addition, the functional modules in the various embodiments of this application can be integrated into one processing module, or each module can exist physically separately, or two or more modules can be integrated into one module.

[0106] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

[0107] In conclusion, the above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A finite element discretization analysis method for elongated foundation pits based on spatiotemporal mechanical coupling, characterized in that, Includes the following steps: Obtain geometric data of the foundation pit, soil creep parameters, and construction sequence time data; The normal displacement constraint conditions at the longitudinal end of the foundation pit are set based on geometric data; The excavation unloading range is determined based on the construction sequence time data, and the time-varying unloading vector is calculated. Based on the constraints and time-varying unloading vector, the total strain increment for each time increment step is calculated, and the corresponding effective stress and creep internal variables are updated. The consistent tangent stiffness matrix is ​​updated based on the effective stress and creep internal variables, and the displacement and stress fields are iteratively solved for each time increment step. Based on the displacement and stress fields at each time increment step, the spatiotemporal distribution results of the foundation pit deformation and the internal forces of the support structure are output.

2. The method according to claim 1, characterized in that, The normal displacement constraint conditions at the longitudinal ends of the foundation pit include: Extract the set of mesh nodes located at the outer boundaries of both longitudinal ends of the foundation pit; Set the normal displacement degree of freedom of each node in the mesh node set to zero.

3. The method according to claim 2, characterized in that, The determination of the excavation unloading range based on construction sequence time data includes: Analyze the construction step sequence time data and extract the longitudinal segment boundary coordinates and excavation depth elevation corresponding to the current time increment step; Traverse the finite element discrete mesh pre-established from the foundation pit geometry data, and extract the elements whose spatial location falls within the range of the longitudinal segment boundary coordinates and the excavation depth elevation as the target element set; Mark the elements in the target element set as removed, and use the spatial region corresponding to the target element set as the excavation unloading range for the current time increment step.

4. The method according to claim 3, characterized in that, The calculation of the time-varying unloading vector includes: Extract the element stress of the target element set before it is marked as removed, and calculate the initial equivalent release nodal force of the target element set on the remaining effective mesh nodes; Extract the creep internal variables of the remaining effective grid cells at the excavation boundary and calculate the creep state evolution increment of the current time increment step; Using the creep state evolution increment as an endogenous driving parameter, the release process of the initial equivalent release node force is controlled. The load release amount of the current time increment step is autonomously allocated by the soil's internal rheological mechanism, and the corresponding time-varying unloading vector is generated.

5. The method according to claim 4, characterized in that, The calculation of the total strain increment for each time increment step based on constraints and time-varying unloading vectors includes: Based on the constraints of time-varying unloading vector and normal displacement, the system equilibrium equation corresponding to the current time increment step is constructed. Using the displacement field obtained from the previous time increment step as the initial displacement guess, the Newton-Raphson iteration is used to solve the system equilibrium equations to obtain the nodal displacement increments for the current time increment step. The nodal displacement increments are converted into the total strain increments of the current time increment step using geometric equations.

6. The method according to claim 5, characterized in that, The updated effective stress and creep intrinsic variables include: Based on the total strain increment, the effective stress of the test is calculated according to the elastic constitutive relation; Calculate the corresponding creep strain increment based on the tested effective stress and soil creep parameters; The effective stress corresponding to the current time increment step is obtained by correcting the creep strain increment and probing the effective stress. Update the creep intrinsic variable corresponding to the current time increment step based on the creep strain increment.

7. The method according to claim 6, characterized in that, The step of updating the consistent tangent stiffness matrix based on effective stress and creep intrinsic variables includes: Extract the additional stiffness contribution generated by the nonlinear evolution of the time-varying unloading vector with creep internal variables, and construct an additional boundary coupling stiffness term; Based on the effective stress and creep internal variables of the current time increment step, a basic consistent tangent modulus matrix is ​​constructed in combination with the creep evolution rule; The additional boundary coupling stiffness term is fused and assembled with the foundation consistent tangent modulus matrix to update the consistent tangent stiffness matrix for the current time increment step.

8. The method according to claim 7, characterized in that, The construction of the additional boundary coupling stiffness term includes: Calculate the first partial derivative of the time-varying unloading vector with respect to the creep internal variable; Calculate the second partial derivative of the creep internal variable with respect to the total strain increment; By combining the first and second partial derivatives and the finite element geometric equations, the non-zero gradient matrix of the external load vector with respect to the nodal displacement increment is obtained analytically using the chain rule, and the non-zero gradient matrix is ​​used as the additional boundary coupling stiffness term.

9. The method according to claim 8, characterized in that, The iterative solution of the displacement and stress fields at each time increment step includes: Within each nonlinear iteration step, the time-varying unloading vector is updated according to the creep internal variable corresponding to the displacement field of the current iteration step. The difference between the updated time-varying unloading vector and the internal nodal force vector obtained from the effective stress is calculated to obtain the global unbalanced residual force vector; A linearized incremental equation is constructed using the uniform tangent stiffness matrix and the global unbalanced residual force vector. The iterative displacement correction vector is then solved to update the nodal displacements. Iterate until the global unbalanced residual force vector satisfies the preset convergence tolerance, and output the displacement field and stress field of the current time increment step.

10. The method according to claim 9, characterized in that, The displacement and stress fields based on each time increment step output the spatiotemporal distribution results of the foundation pit deformation and the internal forces of the support structure, including: The displacement sequence of the target observation node is extracted from the displacement field to obtain the pit deformation at the current time increment step; Extract the integral point stress of the support structure element from the stress field, and obtain the internal force of the support structure at the current time increment step through cross-sectional integration calculation; The excavation pit deformation and the internal forces of the support structure at each time increment step are sequentially assembled according to the calculation time axis to generate and output the distribution results that characterize the spatiotemporal evolution of the entire excavation process.