Numerical simulation method for time-space evolution law of negative friction of high fill pile foundation

Abstract

The present application relates to geotechnical engineering numerical calculation technical field, disclose a kind of high fill pile foundation negative friction time-space evolution law numerical simulation method, the method is first built including no-thickness contact surface element pile soil flow coupling finite element model, during calculation, real-time monitoring load change rate and mapping is procedure state weight factor, to this quantification identification fast loading and static consolidation state;Using the factor simultaneously executes the resistance correction to contact surface local permeability coefficient, and executes the damping update to consolidation strengthening factor, to simulate procedure rate to interface fluid-structure coupling effect and the dynamic regulation and control of shear strength evolution;Based on the updated state variable, construct consistent tangent stiffness matrix for nonlinear iterative solution, finally by tracking the depth of neutral point of shear stress sign inversion, integral obtains negative friction distribution.The present application realizes the reproduction of construction whole process fluid-structure coupling mechanism, improves the calculation precision of high fill engineering pile foundation long-term settlement and stress analysis.

Description

Technical Field

[0001] This invention relates to the field of numerical calculation technology in geotechnical engineering, specifically to a numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations. Background Technology

[0002] In high embankment projects, due to the self-weight consolidation of the backfill soil and the action of the superstructure load, the settlement of the soil along the pile is greater than the settlement of the pile itself, thus generating downward negative skin friction on the pile foundation surface. This poses a challenge to the bearing capacity and structural safety of the pile foundation. Numerical simulation, as an important means of analyzing this complex geotechnical behavior, is widely used to predict the stress state and long-term evolution of the pile-soil interface.

[0003] However, existing numerical simulation methods for pile-soil contact surfaces typically treat permeability characteristics and mechanical strength parameters as relatively independent variables when dealing with fluid-structure interaction problems, or only establish a simple coupling relationship based on changes in void ratio, lacking an endogenous response mechanism to the key dynamic factor of construction loading rate. This limitation makes it difficult for existing models to automatically identify and distinguish between two different physical conditions—rapid filling (approximate undrained shear) and static intermittent (drained consolidation)—within a unified computational framework. Specifically, conventional methods neglect the transient impediment effect of soil skeleton deformation on seepage channels under rapid loading conditions, as well as the time lag effect of effective stress intensity increase relative to the increase in total load. This lack of a physical mechanism makes the simulation results unable to accurately reproduce the accumulation of excess pore water pressure at the contact surface, resulting in significant deviations in the calculation of the spatiotemporal distribution of pile side shear strength and negative skin friction.

[0004] Furthermore, the numerical analysis of high embankment projects spans a vast timescale, from short-term intense construction loading to long-term slow consolidation. Traditional finite element method (FEM) solutions often employ fixed time steps or step-size control methods based solely on mathematical truncation errors, lacking an adaptive and collaborative control mechanism directly related to the physical process states (such as loading or consolidation). This leads to the solver failing to converge when faced with transient responses caused by abrupt changes in the process due to excessively large time steps. Conversely, when dealing with long-term consolidation problems, the inability to effectively scale up the time step results in low computational efficiency, making it difficult to balance numerical stability and computational timeliness in full-lifecycle simulation. Summary of the Invention

[0005] To address the shortcomings of existing technologies, this invention provides a numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations. This method solves the problems of existing technologies being unable to automatically respond to the dynamic control of the pile-soil interface fluid-structure coupling mechanism and shear strength evolution under a unified computational framework. This results in the difficulty in accurately reproducing the excess pore water pressure accumulation effect and strength lag characteristics during the rapid filling period, as well as the difficulty in balancing computational accuracy and solution efficiency in the full life cycle simulation of high-fill piles.

[0006] To achieve the above objectives, the present invention provides the following technical solution: This invention provides a numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations, the method comprising: First, a pile-soil flow coupled finite element model containing elements with no thickness contact surface is constructed, and the state variables are initialized at the integration points of the elements with no thickness contact surface. At the start of the time increment step of the numerical calculation, the process state weight factor is calculated based on the global loading data, and the process state weight factor is substituted into the pile-soil-flow coupled finite element model. Subsequently, the local permeability coefficient of the contact surface of the pile-soil fluid coupling finite element model and the consolidation strengthening factor in the state variables were simultaneously corrected using the process state weight factor to simulate the regulatory effect of the process on fluid-structure coupling and strength growth. Then, based on the updated state variables, the dynamic friction coefficient is synthesized and a consistent tangent stiffness matrix is ​​constructed. Nonlinear iterative solution is performed on the pile-soil-flow coupled finite element model to obtain converged shear stress field data. Finally, the depth location where the shear stress sign reverses is traced based on the shear stress field data as the neutral point, and the negative skin friction distribution is calculated by integration.

[0007] In the numerical simulation method described above, the state variables include cumulative plastic shear work, consolidation strengthening factor, shear damage factor, and the process state weight factor from the previous increment step. During initialization, the cumulative plastic shear work is assigned a value of zero, the consolidation strengthening factor and shear damage factor are assigned unit values, and the process state weight factor from the previous increment step is assigned a value of zero. The thickness-free contact surface element is topologically composed of master-face nodes belonging to the pile foundation surface and slave-face nodes belonging to the soil surface, to simulate the discontinuous deformation of the pile-soil interface.

[0008] Furthermore, regarding the calculation of the process state weighting factor, this method calculates the ratio of the effective load increment to the time increment within the current increment step, and divides this ratio by a preset reference pressure to obtain the dimensionless load change rate. Subsequently, a nonlinear smooth step function model containing a natural exponential term is constructed. This model includes shape control parameters, a process judgment threshold, and a transition zone bandwidth parameter. This nonlinear smooth step function model is used to map the load change rate to the process state weighting factor. When the load change rate is greater than the process judgment threshold, the process state weighting factor approaches a set upper limit value (e.g., 1), representing the undrained shear-dominated state; when the load change rate is less than the process judgment threshold, the process state weighting factor approaches a set lower limit value (e.g., 0), representing the drained consolidation-dominated state.

[0009] In terms of fluid field correction, the process of correcting the local permeability coefficient using the process state weighting factor follows a multiplicative coupling rule. The local permeability coefficient is obtained by multiplying the initial permeability coefficient, the void ratio correction term, and the process retardation correction term. The void ratio correction term is defined as an exponential function of the volumetric strain of the contact surface element without thickness, used to characterize the influence of soil skeleton geometric deformation on the pore size of the seepage channels. The process retardation correction term is defined as a linear function of the process state weighting factor and the preset process retardation coefficient; when the process state weighting factor approaches the set upper limit value, the value of the process retardation correction term decreases, forcibly reducing the local permeability coefficient, thereby simulating the accumulation effect of transient excess pore water pressure at the contact surface under rapid loading conditions.

[0010] In terms of solid field correction, the consolidation strengthening factor is updated with damping using a process state weighting factor. This step first reads the excess pore water pressure values ​​and initial peak excess pore water pressure of adjacent soil elements in the pile-soil coupled finite element model to calculate the local degree of consolidation. The theoretical target strengthening value is set as a power function of this local degree of consolidation, including the maximum consolidation gain coefficient and the consolidation hardening exponent. The update amplitude is determined based on the process state weighting factor: when the process state weighting factor approaches the set upper limit, the change in the consolidation strengthening factor is suppressed, maintaining its value from the previous moment, simulating the phenomenon that strength cannot increase under undrained conditions; when the process state weighting factor approaches the set lower limit, the consolidation strengthening factor is controlled to approach the theoretical target strengthening value, simulating the strength increase brought about by drained consolidation.

[0011] Furthermore, the method also includes a step of updating the shear damage factor based on accumulated plastic shear work. The result is calculated by adding the dot product of the shear stress vector and the plastic relative displacement increment vector within the current time increment step to the accumulated plastic shear work. The shear damage factor at the current moment is calculated using a negative exponential decay function that includes parameters for residual damage ratio and damage rate. This causes the shear damage factor to nonlinearly decay from a unit value to the residual damage ratio as the accumulated plastic shear work increases, thus reflecting the shear softening characteristics under large deformation.

[0012] In terms of the solution mechanism, a dynamic friction coefficient is synthesized based on the updated state variables, and the consistent tangent stiffness matrix is ​​updated accordingly. The dynamic friction coefficient is defined as the product of the initial static friction coefficient, the consolidation strengthening factor, and the shear damage factor. Based on the yield criterion, if the state is plastic, a consistent tangent stiffness matrix including a plasticity correction term is constructed. This matrix is ​​obtained by subtracting the plasticity projection component from the elastic stiffness matrix. The dynamic hardening modulus, as the denominator of the plasticity projection component, is obtained by weighted summation of the partial derivatives of the consolidation strengthening factor and the shear damage factor with respect to the plastic displacement. This ensures the consistency between the tangent stiffness matrix and the constitutive integration algorithm, guaranteeing the convergence rate of the nonlinear iteration.

[0013] Furthermore, this method also performs adaptive time-step collaborative control and long-term calculation. The process state weighting factor is fed back to the time-step controller of the global solver of the pile-soil-flow coupled finite element model. A linear growth model for the next increment time span is constructed based on the process state weighting factor: when the process state weighting factor approaches a set upper limit, the growth of the time increment step is limited to accurately capture transient responses; when the process state weighting factor approaches a set lower limit, the time increment step is allowed to be amplified to accelerate the long-term consolidation calculation process. During this period, the method activates a collaborative calculation mechanism for secondary consolidation rheology and consolidation damping: the additional viscous strain increment generated by secondary consolidation is calculated based on the logarithmic time law to simulate the long-term creep of the soil skeleton under constant effective stress; simultaneously, a hydraulic and mechanical coupled damping term is superimposed during the effective stress update process. This damping term is defined as a damping stress proportional to the time derivative of pore water pressure and the time increment step, used to suppress numerical oscillations under large time steps.

[0014] In the post-processing stage, the depth location where the shear stress sign reverses is used as the neutral point. Specifically, this involves extracting the shear stress components distributed along the pile axis in the pile-soil-flow coupled finite element model, retrieving the sequence of two adjacent element nodes where the shear stress component sign changes, and using a linear interpolation algorithm to calculate the depth corresponding to the zero shear stress component, which is recorded as the instantaneous neutral point depth at the current moment. When integrating to calculate the negative skin friction distribution, all thickness-free contact surface elements located above the instantaneous neutral point depth are traversed. For each element, the product of the pile perimeter, the element axial length, and the absolute value of the axial shear stress at the element integration point is calculated, and all products are summed to obtain the total downward load borne by the pile.

[0015] This invention provides a numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations. It has the following beneficial effects: 1. This invention establishes a process state weight factor mapping mechanism based on load change rate, and uses a single factor to synchronously correct the local permeability coefficient and consolidation strengthening factor of the contact surface. This fluid-solid bidirectional collaborative correction strategy can automatically identify and distinguish between rapid loading and static consolidation conditions, solving the problem that traditional methods are difficult to take into account the dynamic influence of different construction rates on the seepage channel and shear strength evolution of the pile-soil interface in a unified model, and improving the physical reality of multi-physics coupling simulation.

[0016] 2. This invention constructs a process retardation correction model for local permeability and a damping update model for consolidation strength. In the undrained dominance stage of rapid loading, by forcibly reducing the permeability coefficient and suppressing the strength growth, the accumulation effect of transient excess pore water pressure and the lag phenomenon of the growth of effective stress intensity of soil are reproduced. This enables a more accurate calculation of the spatiotemporal distribution characteristics of pile side negative skin friction during the settlement of high embankments, avoiding calculation deviations caused by ignoring the accumulation of excess pore pressure.

[0017] 3. This invention utilizes process state weighting factors to adaptively and collaboratively control the global time step, and combines it with a nonlinear iterative algorithm based on a consistent tangent stiffness matrix. This mechanism automatically densifies the time step to ensure computational convergence when capturing transient load responses, and automatically amplifies the time step to improve computational efficiency during the long-term consolidation stage. This enables efficient and stable solutions to the evolution law of negative skin resistance of high-fill pile foundations from the construction period to the operation period. Attached Figure Description

[0018] Figure 1 This is an overall flowchart of the numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations provided in an embodiment of the present invention; Figure 2 This is a flowchart illustrating the calculation and mapping logic of the process state weight factor in this embodiment of the invention. Figure 3 This is a flowchart illustrating the constitutive integral and parameter update logic based on the process state weight factor in this embodiment of the invention. Figure 4 This is a flowchart of the post-processing for neutral point identification and negative friction integral calculation in an embodiment of the present invention. Detailed Implementation

[0019] The technical solutions in 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.

[0020] See attached document Figure 1 This invention provides a numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations. This method runs on a computer device containing a processor and memory, and is implemented through finite element analysis software and a user-defined subroutine interface. The method includes the following steps: First, a pile-soil fluid-flow coupled finite element model incorporating the geological characteristics of high embankments is constructed. Based on engineering survey data and design drawings, this step establishes a three-dimensional geometric model including the undisturbed soil layer, pile foundation structure, and layered fill soil. During discretization, pore pressure and displacement coupled elements are used for the soil and fill to calculate the fluid-structure interaction effect during consolidation; linear elastic solid elements are used for the pile foundation. Specifically, a thickness-free contact surface element is set between the pile foundation and the soil, and a set of state variables is initialized at the integration points of this thickness-free contact surface element. These state variables include at least the cumulative plastic shear work. Consolidation strengthening factor Shear damage factor and the process state weighting factor of the previous incremental step. Among them, cumulative plastic shear work The initial value is set to zero, and the consolidation strengthening factor is... With shear damage factor The initial value is set to a unit value.

[0021] Next, at the start of each time increment step in the numerical computation, real-time identification and weight mapping of the process status are performed. The processor accesses the globally loaded data and calculates the rate of change of the fill load within the current increment step. The load change rate is defined. This is the ratio of the increment of vertical stress applied to the fill per unit time to the reference pressure. Based on this load variation rate... The process state weighting factor at the current moment is calculated using a smooth step function. Process status weighting factor This is used to quantify the degree to which the current calculation step is in a rapid loading or intermittent consolidation process. The calculation logic of this smooth step function is as follows: construct about A nonlinear function that includes shape control parameters that control the steepness of the transition curve. Preset process judgment thresholds and transition zone bandwidth parameters This function guarantees It changes continuously between 0 and 1, when Significantly greater than hour, Approaching 1 represents the undrained shear-dominated state; when Significantly smaller than hour, A value close to 0 indicates a state dominated by drainage consolidation.

[0022] Subsequently, based on the calculated process state weighting factor The constitutive integral update, which involves bidirectional coupling of shear and seepage, is performed at the integration point of the contact surface element without thickness. This process includes dynamic correction of the permeability coefficient and damped update of the constitutive state variables. The processor first corrects the local permeability coefficient near the interface. The correction logic sets the local permeability coefficient. With initial permeability coefficient Porosity correction term And related to process hindrance correction items. Among them, the porosity correction item. It is about volumetric strain The exponential function; the process stagnation correction term is a weighting factor related to the process state. and the preset process resistance coefficient The linear interpolation function. When When the value approaches 1, this logic forces a decrease in the local permeability coefficient, thereby simulating the accumulation effect of transient excess pore water pressure in the solver. Simultaneously, the processor uses process state weighting factors... consolidation strengthening factor Perform damping update. Calculate the theoretical target enhancement value. This value is set as the local degree of consolidation. The power function; then according to Determine the update magnitude for the current step, when When it is large, The update is suppressed, maintaining the value from the previous moment, reflecting the lag in intensity growth under rapid loading; when When smaller, Approaching Shear damage factor Then, based on the cumulative plastic shear work of the current step Updated according to the negative exponential decay law.

[0023] After updating the state variables, the dynamic friction coefficient is synthesized and the uniform tangent stiffness matrix is ​​updated. The processor calculates the dynamic friction coefficient at the current time and depth based on the updated state variables. This coefficient is defined as the initial static friction coefficient. Consolidation strengthening factor With shear damage factor The product of the three. Based on this dynamic friction coefficient. Given the current effective normal stress, calculate the ultimate shear stress and perform yield testing. If the state is slipping, apply a back-bias correction to the stress. During this process, calculate the partial derivative of the shear stress with respect to the shear displacement to obtain the consistent tangent stiffness matrix. Since... and under its control and All of them are continuously differentiable functions, which ensures the continuity of the stiffness matrix, thereby ensuring the numerical convergence of the nonlinear iteration.

[0024] Simultaneously, adaptive time-step cooperative control is executed. The processor assigns process state weighting factors. The time step controller is fed back to the global solver. When As the value approaches 1, the solver limits the growth of the time increment step, using a small step size to capture the transient response; when As the value approaches 0, the solver allows for scaling up the time increment step to accelerate the calculation process of the consolidation period.

[0025] Furthermore, to accurately simulate the spatiotemporal evolution characteristics of high-fill projects during the long-term post-construction operation phase, this invention introduces a long-term secondary consolidation rheology and consolidation damping co-calculation mechanism. When the process state weighting factor... When the system remains in a near-zero drained consolidation state for an extended period, the processor calculates the additional viscous strain increment caused by the secondary consolidation based on the logarithmic time law. The calculation formula is as follows: ; in, The secondary consolidation coefficient is determined by long-term indoor consolidation tests. The initial void ratio; This is the current cumulative consolidation time; This is the time increment for the current computation step, and this value is dynamically amplified by an adaptive algorithm during long-term computation. This is the equivalent projection vector. This formula quantifies the creep deformation of the soil skeleton under constant effective stress over time after the excess pore water pressure has dissipated, thus perfecting the long-term evolution law of negative skin friction in the time dimension.

[0026] Meanwhile, to address the potential numerical oscillations caused by large time steps in long-term calculations, the processor introduces a hydraulically and mechanically coupled damping term into the consolidation equations. This consolidation damping is defined as a numerical viscous force capable of suppressing non-physical high-frequency pore pressure fluctuations. When updating the effective stress, the processor superimposes a damping stress proportional to the rate of change of pore pressure. The calculation formula is as follows: ; in, is the consolidation damping coefficient, which is a function of the mesh feature size and the permeability coefficient; This is the time increment for the current computation step; The time derivative of pore water pressure; It is a second-order unit tensor.

[0027] The physical significance of this damping mechanism lies in the fact that, during long-term computation, when the numerical solver attempts to make a large time jump (i.e., When the value is large, This generates a small reverse correction stress, simulating the viscous resistance of pore fluid in the microscopic channels of the soil. This not only smooths the pore pressure dissipation curve over a long period, preventing computational divergence, but also more realistically reflects the damping effect brought about by fluid viscosity during soil consolidation, ensuring the numerical stability and physical realism of the spatiotemporal evolution simulation throughout its entire life cycle.

[0028] Finally, after the global solution converges, the result data is extracted and the spatiotemporal evolution law of negative skin friction is reconstructed. The processor outputs a four-dimensional data matrix containing time, depth, pile side shear stress, and pore water pressure. By tracking the depth position where the sign of the pile side shear stress reverses at each moment, the dynamic drift trajectory of the neutral point is determined; by integrating the dynamic friction coefficient evolution data at different depths, a negative skin friction distribution cloud map reflecting the influence of the filling process is generated.

[0029] When constructing the multiphysics coupling model, a geometric entity is first generated in a three-dimensional Cartesian coordinate system based on the geological stratification data and pile foundation design parameters from the engineering survey report. This geometric entity includes the underlying undisturbed soil-bedrock composite, the pile foundation entity penetrating the soil layers, and the fill body above the undisturbed soil surface. To accurately simulate the loading history during the construction of high fills, the fill body is discretized into several independent horizontal fill layers during the geometric modeling stage. The thickness of each fill layer corresponds to the single fill height or compaction layer thickness in the actual project, so that the corresponding unit set can be activated according to the time step sequence in subsequent calculations.

[0030] In the mesh discretization step, specific element types are assigned to the physical properties of different media to construct a computational domain capable of responding to the shear-seepage coupling mechanism. For undisturbed soil layers and layered fill bodies, pore pressure-displacement coupling elements are selected for discretization. Each node of this type of element has translational degrees of freedom in the X, Y, and Z directions, as well as one pore water pressure degree of freedom, thus enabling the simulation of the interaction between soil skeleton deformation and pore fluid seepage. To capture the drastic stress gradient changes in the soil around the pile, a radially gradient meshing strategy is adopted with the pile axis as the center. That is, a high-density refined mesh is set in the region immediately adjacent to the pile body, and the mesh size gradually increases geometrically with the radial distance to balance computational accuracy and efficiency. For the pile foundation entity, since its permeability is much lower than that of the soil and it mainly bears structural loads, linear elastic solid elements are selected for discretization. The nodes of this linear elastic solid element only contain three translational degrees of freedom and do not include pore water pressure degrees of freedom.

[0031] Specifically, to achieve refined simulation of the complex mechanical behavior of the pile-soil interface, this invention constructs a thickness-free contact surface element between the nodes on the outer surface of the pile foundation and the nodes on the inner surface of the adjacent soil. Topologically, this thickness-free contact surface element consists of a series of node pairs, each pair containing a master node belonging to the pile foundation surface and a slave node belonging to the soil surface. In the initial undeformed state, the master and slave nodes coincide in spatial geometry but belong to different node numbering sequences in the finite metadata structure. This thickness-free contact surface element defines a local coordinate system, encompassing compression and separation behaviors along the normal direction and relative slip behaviors along the tangential direction.

[0032] Furthermore, considering the hydraulic boundary conditions in high-fill environments, differentiated values ​​were assigned to the model boundaries. For the bottom surface of the model, displacement degrees of freedom in all directions were constrained, and pore water pressure degrees of freedom were closed to simulate impermeable bedrock; for the sides of the model, normal displacement degrees of freedom were constrained and set as impermeable boundaries; for the top surface of the fill and undisturbed soil, a free drainage boundary with zero pore water pressure was set. More importantly, at the pile-soil interface, the integration points of the contact elements without thickness were assigned dynamically adjustable permeability properties. These permeability properties allow pore fluid to cross the interface tangentially or normally under specific hydraulic gradient conditions, thus making the interface not only a boundary for mechanical transmission but also a channel or barrier for fluid seepage. This provides a mesh topology basis for subsequent calculations that change the interface permeability based on the process state.

[0033] After mesh discretization and boundary condition setting, to support nonlinear constitutive evolution calculations based on process fingerprinting, an independent state variable storage space needs to be allocated at the integration point location of each thicknessless contact surface element. This state variable storage space is address-mapped through the user subroutine interface of the finite element solver, forming a state variable array that corresponds one-to-one with the integration point. During the preprocessing stage before the calculation begins, the processor performs physical definition and numerical initialization for each component in this array to establish the initial benchmark for subsequent incremental step iterations.

[0034] First, define and initialize the cumulative plastic shear work. It is used to record the total energy dissipated by irreversible plastic slip at the pile-soil interface during historical loading. At the initial moment, since no shear displacement has occurred, the accumulated plastic shear work at the integration points of all thickness-free contact elements is recorded. The value is assigned to zero. This variable serves as the core parameter driving the evolution of shear damage in subsequent calculations, and its value monotonically increases with the product of the plastic shear displacement generated in each increment step and the current shear stress.

[0035] Secondly, define and initialize the consolidation strengthening factor. This variable characterizes the effect of excess pore water pressure dissipation in the soil on the enhancing effect of interfacial shear strength. A consolidation strengthening factor is defined. The initial value of is 1, representing the baseline strength state before excess pore water pressure dissipation occurs. During the calculation, this variable will be dynamically updated by reading the pore water pressure degrees of freedom values ​​of soil elements at the same location or adjacent soil elements, combined with the degree of consolidation function, to reflect the enhanced interlocking effect of soil particles caused by the increase in effective stress.

[0036] Next, define and initialize the shear damage factor. This variable describes the structural wear and strength degradation at the pile-soil interface caused by large deformation shear under high confining pressure. A shear damage factor is set. The initial value is 1, representing that the interface is in a intact, undamaged state. With the accumulation of plastic shear work... As the value increases, this variable will gradually decrease according to a preset negative exponential law, simulating the degradation process of interface roughness.

[0037] At the same time, define and initialize the process state weighting factor of the previous incremental step. This variable stores the process state determination result at the end of the previous incremental step to support the time continuity smoothing in the constitutive algorithm. Its initial value is set to zero, representing a default initial state of rest or complete drainage. Furthermore, to calculate the local consolidation degree, the initial peak excess pore water pressure also needs to be defined. The variable is initialized to zero, and the maximum excess pore water pressure generated at the current integration point is captured and locked at the instant the filling load is applied, serving as the denominator for subsequent calculations of the degree of consolidation normalization. Through the precise definition and initialization of the above state variable array, a microscopic data structure of the pile-soil interface capable of remembering the loading history and evolution path is constructed.

[0038] At the start of each incremental step in the numerical simulation, the processor first executes a global monitoring command. This command accesses the global loading state data of the current model through the finite element solver's data exchange interface. Specifically, the processor identifies and extracts the vertical filling load values ​​acting on the top surface of the soil. To accurately capture the dynamic characteristics of the construction process, this monitoring process distinguishes the time increment of the current incremental step. And the increment of external loads applied during that time span.

[0039] The processor calculates the stress increase rate per unit time. First, it reads the total filling load value at the end of the current increment step and the total filling load value at the beginning of the current increment step; the difference between the two is the effective load increment for the current increment step. This effective load increment represents the additional vertical stress applied to the foundation surface during the current simulation period due to the laying of new fill layers or mechanical compaction. Subsequently, the effective load increment... Divide by the current time increment Thus, the stress change rate with physical dimensions is obtained.

[0040] To eliminate the dimensional effects of varying project scales and geological depths and establish a unified process discrimination standard, the processor further performs dimensionless processing on the aforementioned stress variation rate. The load variation rate is defined. Its calculation logic is to divide the calculated stress change rate by the preset reference pressure. The reference pressure As a constant normalization factor, it is used to convert stress rate, which has physical units, into a purely numerical index. Through the above calculation steps, the processor outputs in real time a load change rate reflecting the current construction intensity. When the filling construction is in the rapid loading phase, due to... shorter and Larger, calculated The values ​​increase significantly; however, during the intermittent filling period or the simple consolidation and load-bearing period, since the external load no longer increases, Approaching zero, leading to The value approaches zero. This monitoring result... It is directly used as the input parameter for the subsequent calculation of the process status weight factor.

[0041] See attached document Figure 2 In obtaining the load change rate The processor then performs a nonlinear mapping operation to transform it into a normalized process state weight factor. This mapping process does not employ traditional binary logic judgments, but instead constructs a continuously differentiable smooth step function model. This smooth step function model aims to establish a quantitative correlation between the load application rate and the soil drainage conditions, while eliminating the iterative non-convergence phenomenon caused by abrupt state changes in numerical calculations.

[0042] Specifically, process status weighting factor It is calculated using a fractional expression that includes a natural index term, and the formula is as follows: ; in, Represented by the natural constant Operations on exponential functions with base 0; The process judgment threshold is a dimensionless scalar whose physical meaning corresponds to the critical loading rate at which the soil changes from a drained state to an undrained state. The specific value is determined by the characteristic ratio of the soil's permeability coefficient to the drainage path length. This is the transition zone bandwidth parameter, used to set the range of the mixed calculation interval near the threshold, ensuring that the mapping curve transitions smoothly on both sides of the threshold rather than jumping vertically; This is a shape control parameter used to adjust the slope of the tangent line in the transition zone of the S-curve. The larger the value of this parameter, the more sensitive the state switching; the smaller the value, the wider the transition zone.

[0043] Based on the above mapping algorithm, when the calculated load change rate Significantly higher than the set process judgment threshold When the denominator of the exponential term approaches 0, the process state weighting factor... Asymptotically approaching a value of 1 indicates that the current computational step is in a rapid loading process, and the physical field response is dominated by undrained shear mechanisms. Conversely, when... Significantly below the process judgment threshold When the exponential term approaches infinity, the process state weight factor becomes... The asymptotic approach to 0 indicates that the current period is an intermittent period or an extremely slow loading period, and the physical field response is dominated by the drainage consolidation mechanism.

[0044] In particular, when Falling into Centered on, with a width of When within the neighborhood, the process state weight factor Output an intermediate value between 0 and 1. This mathematical construction ensures the weighting factor of the process state. Regarding the rate of change of load The first derivative is continuous and bounded throughout the entire domain. This continuity ensures the smooth assembly of the full stiffness matrix in the finite element solver, allowing subsequent calculations based on... The smoothness of the modified constitutive model in the tangential direction avoids residual force oscillations caused by non-physical jumps in state variables, thus ensuring that the Newton-Raphson iterative algorithm can converge stably along a smooth path.

[0045] After calculating the process state weight factor for the current increment step Afterward, the processor does not immediately proceed to the next incremental step of stiffness matrix assembly, but instead first executes adaptive cooperative control of the time step. This step aims to resolve the contradiction between the need for extremely small time steps in the fast loading process to capture transient responses, and the need for larger time steps in the long-term consolidation process to improve computational efficiency. The processor will then calculate... The numerical feedback is sent to the global time stepping controller of the finite element solver, which determines the time span of the next incremental step through a step size scaling algorithm based on process weights. .

[0046] Specifically, the time span of the next incremental step Based on the current time span Weighting factors related to process status The linear growth model is constructed, and its calculation formula is as follows: ; in, The maximum acceleration factor for the time step is a dimensionless constant greater than zero, used to set the upper limit of the allowable time step growth rate under the fully drained consolidation state. This characterizes the degree to which the current numerical model is in a consolidation-dominated state.

[0047] Based on the above control logic, when the process state weight factor When the value approaches 1, indicating that the numerical model is in a high-intensity filling loading stage, the growth term in the formula approaches zero. At this point, the time span of the next incremental step... The constraint is set to prevent growth or maintain only the current magnitude, thus forcing the solver to integrate using fine time slices. This constraint ensures that the transient accumulation peak of pore water pressure and the minute slip at the pile-soil interface can be accurately captured during rapid loading, avoiding the accumulation of numerical integration errors caused by excessively large step sizes.

[0048] Conversely, when the process status weight factor When the value approaches zero, that is, when the system enters a long period of intermittent consolidation or settling, the growth term in the formula reaches... At this point, the solver is allowed to scale the time step exponentially. This strategy enables the simulation of long-duration physical evolution processes, such as the dissipation of excess pore water pressure, to be completed with fewer computational steps, thus significantly reducing computational resource consumption. Furthermore, the processor incorporates a mutation truncation mechanism that detects... If the change between two adjacent increment steps exceeds a preset safety threshold, the next increment step will be forcibly reset to the preset minimum time step size. This is to prevent convergence failure due to excessive step size during process switching.

[0049] During the constitutive integration of the pile-soil interface contact element, the processor calculates the dynamic friction coefficient at the current moment and location by synthesizing the basic physical properties and evolving state variables. This step aims to overcome the limitation of the friction coefficient being a constant in the traditional Coulomb friction model and establish a dynamic strength model that can simultaneously reflect the time effect (consolidation) and the deformation effect (shear softening).

[0050] Specifically, the processor calculates the dynamic friction coefficient based on the multiplication coupling rule. The calculation formula is as follows: ; in, The initial static friction coefficient is determined by laboratory direct shear tests. It represents the basic tribological property of the pile side surface and soil when there is no damage and no consolidation reinforcement. It is used as a dimensionless reference constant input to the model. The consolidation strengthening factor is a dimensionless variable greater than or equal to 1, used to characterize the contribution of soil excess pore water pressure dissipation and soil particle rearrangement to the interfacial shear strength over time. The shear damage factor is a dimensionless variable less than or equal to 1, used to characterize the shear strength attenuation effect caused by interfacial roughness wear or cementation failure as the cumulative plastic slip increases.

[0051] In this synthesis logic, the consolidation strengthening factor With shear damage factor This constitutes a competitive coupling mechanism. When the fill material is statically consolidated, Numerical growth was dominant, leading to The increase simulates the time-dependent growth of the bearing capacity of the pile foundation; when the pile foundation experiences significant relative slippage... The numerical decay is dominant, leading to The decrease simulates the softening phenomenon of the interface. The processor independently performs the multiplication operation at each integration point of each increment step, thereby forming a non-uniformly distributed friction coefficient field that evolves over time along the entire length of the pile.

[0052] Calculated dynamic friction coefficient This is then substituted into the limit state equation of Coulomb's friction law. The processor utilizes the normal effective stress at the current integration point. Multiplying this by the dynamic friction coefficient yields the interfacial ultimate shear stress at the current moment. The ultimate shear stress at this interface The threshold value used for subsequent elastic yielding determines whether the interface is in a bonded or slipping state. Through the above synthesis logic, the constitutive model achieves a quantitative mapping from the evolution of microscopic state variables to the macroscopic interface mechanical response.

[0053] When calculating the damage component in the dynamic friction coefficient of the pile-soil interface, the processor executes a shear damage accumulation and update algorithm based on the principle of energy dissipation. This algorithm uses the energy of irreversible plastic deformation as an internal variable driving the damage evolution, thereby quantifying the physical wear process of interface roughness.

[0054] First, the processor calculates the increment of plastic shear work generated within the current time increment step. This calculation process extracts the shear stress vector at the integration point of the current thicknessless contact element. and the corresponding plastic relative displacement increment vector The shear stress vector is multiplied by the plastic relative displacement increment vector to obtain the scalar plastic work increment. This increment is then added to the accumulated plastic shear work stored in the state variable array. In this process, real-time tracking of historical energy dissipation during loading is achieved.

[0055] Next, based on the updated cumulative plastic shear work The processor uses a negative exponential decay function to calculate the shear damage factor at the current moment. The mathematical expression for this function is constructed as follows: ; in, The residual damage ratio is a dimensionless constant between 0 and 1, representing the final residual ratio of the friction coefficient relative to the initial value after the interface undergoes infinite shear deformation. This is a damage rate parameter, with the dimension of the reciprocal of energy, used to control the rate at which the damage factor decays with increasing plastic work. Represented by the natural constant Operations on exponential functions with base 0.

[0056] According to this evolutionary model, in the initial unloaded state, the accumulated plastic shear work... The calculated shear damage factor is 0, with the exponent term being 1. A value of 0 indicates that the interface is in its initial, structurally intact state. As relative slippage between the pile and soil occurs, the accumulated plastic shear work... The shear damage factor increases monotonically, with the exponential term gradually approaching 0, leading to... It decays non-linearly from 1 to .

[0057] This calculation process is executed in each iteration step, and the processor calculates... The value is used in real time for the synthesis calculation of the dynamic friction coefficient. When the numerical model simulates the large deformation and settlement caused by high embankments, if significant plastic slippage occurs at a certain location on the pile side, the value at that location... The value will drop rapidly, thereby reducing the shear strength at that point in subsequent calculations, thus simulating the mechanical behavior of the pile-soil interface softening from peak strength to residual strength.

[0058] When calculating the strengthening component in the dynamic friction coefficient of the pile-soil interface, the processor executes a consolidation strengthening update algorithm based on the process state retardation effect. This algorithm aims to solve the problem that traditional models rely solely on pore water pressure values ​​while neglecting the influence of loading rate on the time lag of soil particle rearrangement. By introducing a process state weighting factor, it achieves dynamic control of the strength growth process.

[0059] First, the processor determines the current local consolidation state. The processor accesses the integration point of the soil element topologically adjacent to the element with no thickness contact surface, reads the current pore water pressure value, and calculates the current excess pore water pressure by combining this value with the hydrostatic pressure at that location. Simultaneously, the processor accesses the initial peak excess pore water pressure experienced at that location during the loading history, stored in the state variable array. Based on these two variables, the local degree of consolidation is calculated. It is defined as 1 minus the ratio of the current excess pore water pressure to the initial peak excess pore water pressure.

[0060] Next, based on the calculated local consolidation degree The processor calculates the theoretical target enhancement value without considering the time lag effect. This target value represents the maximum potential reinforcement ratio that the soil structure can achieve after consolidation and reorganization at the current effective stress level. The calculation is constructed using a power function, and its formula is as follows: ; in, The maximum consolidation gain coefficient is a dimensionless constant greater than zero, representing the value when the excess pore water pressure is completely dissipated (i.e., When the value is equal to 1, the maximum percentage increase in the interfacial friction coefficient relative to the initial value; It is the consolidation hardening index, used to control the nonlinear rate of strength growth relative to the change in degree of consolidation.

[0061] Subsequently, to reflect the inhibitory effect of the rapid loading process on the strength growth of the soil structure, the processor utilizes the process state weighting factor of the current increment step. Damping control is applied to the renewal process of the strengthening factor. The current consolidation strengthening factor... Based on the consolidation strengthening factor of the previous incremental step Strengthening value towards theoretical target The incremental update is performed using the following formula: ; Under this damped update mechanism, when the process state weight factor When the value approaches 1, corresponding to the rapid filling and loading process, the update term in the formula approaches 0. This makes the consolidation strengthening factor... The values ​​forced to remain at the previous moment numerically simulate the freezing effect, reflecting that under undrained shear conditions, although pore water pressure may fluctuate slightly due to shear dilatation, the soil skeleton has not yet completed its strengthening through creep. Conversely, when... When the value approaches 0, corresponding to the intermittent consolidation process, the damping term in the formula disappears. able to quickly approach The time-dependent recovery of pile side friction during the static placement period was simulated.

[0062] Finally, the processor calculates the consolidation strengthening factor. Perform a monotonicity constraint check. Ensure the current step... The value must be no less than the value of the previous increment and must always be greater than or equal to 0, to conform to the irreversible physical law of soil consolidation and hardening. (Verified) The values ​​are updated in the state variable array and used as multipliers in the subsequent synthesis calculation of the dynamic friction coefficient.

[0063] During the finite element calculation of the seepage field, the processor, at the beginning of each increment step, calculates the local permeability coefficient at the integration point of the element with no thickness contact surface. Real-time dynamic correction is performed. This correction strategy aims to establish a dynamic hydraulic boundary model that can simultaneously respond to the microscopic volumetric deformation of the soil and the macroscopic construction sequence rhythm.

[0064] Specifically, the processor calculates the local permeability coefficient at the current moment based on the multiplication coupling rule. The calculation formula is constructed as follows: ; in, The initial permeability coefficient represents the baseline water conductivity of the soil at the contact surface under initial unloaded conditions. This is a porosity correction term, which relates to volumetric strain. The exponential function is used to characterize the physical effect of soil skeleton geometric deformation on the pore size of seepage channels; The process state weighting factor is calculated for the current incremental step; The preset process inhibition coefficient is a dimensionless constant between 0 and 1 (e.g., a value of 0.95 to 0.99), used to quantify the degree to which rapid construction inhibits seepage discharge.

[0065] Among them, the porosity correction term The specific mathematical expression is as follows: ; in, The volume sensitivity index is a constant greater than zero; For volumetric strain, a positive value is defined for shear expansion (volume increase) and a negative value for shear compression (volume decrease). When shear expansion occurs at the interface, A value greater than 1 indicates expansion of seepage channels due to loosening of soil particles; when shear compression occurs at the interface, A value less than 1 simulates the seepage barrier effect caused by pore closure.

[0066] More importantly, in the formula This term constitutes the process stagnation correction term. The processor uses this term to map the time-dimensional process characteristics to the spatial-dimensional hydraulic parameters. When the fill is in the rapid loading phase, the process state weighting factor... It approaches the value of 1. At this point, the value of the hindrance correction term decreases significantly (approaching the value of 1). Forced reduction of local permeability coefficient This numerical processing simulates the transient lock-in effect in a very short loading time, where the fluid cannot be discharged in time, thus forcing the finite element solver to calculate the accumulation of high-amplitude transient superpore water pressure inside the thin contact surface element.

[0067] Conversely, when the process status weight factor When the value approaches 0, it corresponds to the numerical model being in the intermittent consolidation or static simulation stage. At this point, the resistance correction term reverts to 1, and the local permeability coefficient... No longer inhibited by process hindrance, but only by the porosity correction term. Control. This allows the accumulated excess pore water pressure to dissipate gradually according to the actual physical permeability of the soil, thereby accurately simulating the smooth transition from undrained shear to drained consolidation. Through the above-mentioned bidirectional coupling update rule, this invention achieves adaptive simulation of complex seepage and stress coupling responses at different construction rates under a single set of model parameters.

[0068] See attached document Figure 3 In each global iteration step of the finite element solver, the processor executes a local stress update algorithm for the integration point of each thicknessless contact element. This algorithm employs an implicit backward Euler integration scheme combined with a regression mapping strategy, aiming to update the stress state and internal state variables of the previous time step to the true stress state at the current moment based on the displacement increment of the current step, and ensuring that the updated stress state satisfies the dynamically evolving yield criterion.

[0069] Specifically, the constitutive integration process begins with the elastic prediction step. The processor assumes that all deformations within the current increment step are purely elastic and calculates the trial shear stress vector accordingly. The calculation formula is as follows: ; in, This is the shear stress vector at the end of the previous convergence increment step; The tangential elastic stiffness coefficient is the element with zero thickness in the contact surface. This parameter is determined by the shear modulus of the interface soil and the numerical penalty stiffness factor. This is the tangential relative displacement increment vector calculated in the current iteration step.

[0070] Subsequently, the processor constructs a dynamic yield surface function to perform plasticity determination. At this point, the processor calls the dynamic friction coefficient that was calculated in the preceding steps. (This coefficient already includes the factors from...) Regulated consolidation strengthening factor and shear damage factor (Combined effects) and the current effective normal stress Dynamic yield function Defined as the difference between the modulus of the trial shear stress and the current shear strength, i.e.: ; in, Represents the Euclidean norm operation of vectors.

[0071] If the calculated dynamic yield function If the value is less than or equal to zero, it indicates that the current increment step is indeed within the elastic range, and the processor directly transmits the probe shear stress vector. Confirmed as the true shear stress vector at the current moment And it does not produce new plastic slip.

[0072] If the calculated dynamic yield function A value greater than zero indicates that the tested stress point is outside the yield surface, violating physical permissibility, and plastic correction must be performed. The processor performs a radial regression operation to pull the stress point back to the dynamic yield surface. The corrected true shear stress vector. The calculation formula is as follows: ; This correction step ensures that the updated shear stress direction is consistent with the test stress direction, but its modulus is strictly limited to the current ultimate shear strength envelope, which is jointly determined by the process state and historical damage. Simultaneously, the processor calculates the plastic slip increment generated during this process and feeds it back into the shear damage accumulation model for damage factor updates in the next iteration, thus achieving a strongly coupled iterative solution between stress updates and state variable evolution.

[0073] To ensure that the Newton-Raphson iterative algorithm maintains a quadratic convergence rate during the nonlinear finite element solution process, the processor must calculate the tangent stiffness matrix, which is mathematically identical to that of the implicit integral algorithm, after completing the local stress update. This matrix is ​​the linearized derivative of the local constitutive integral operator, used to provide the global solver with accurate gradient information about the variation of the residual force vector with displacement increments, rather than directly using the material's physical elastoplastic matrix.

[0074] Specifically, the processor first determines whether the current integration point is in an elastic or plastic state based on the yield determination result in the local constitutive integration process. If the yield function value is less than or equal to zero, it is determined to be in an elastic state, and the processor directly outputs the preset elastic stiffness matrix as a uniform tangential stiffness matrix. In this case, the matrix is ​​represented by diagonal elements that are all tangential elastic stiffness coefficients. The diagonal matrix represents the linear proportional relationship between stress increment and displacement increment within the elastic range.

[0075] If the yield function value is greater than zero, the condition is determined to be plastic. The processor needs to calculate the uniform tangent stiffness matrix, which includes the plasticity correction term. The matrix is ​​constructed following the linearization rules based on the radial regression algorithm, and its calculation formula is as follows: ; in, This represents the second-order unit tensor, used to maintain the basic elastic response structure; It is the tangential elastic stiffness coefficient; It is the unit normal vector in the direction of plastic flow, and its direction is consistent with the test shear stress vector; This represents the tensor product operation. This constitutes a projection operator used to extract the stiffness component along the flow direction; The dynamic hardening modulus is used to quantify the rate of change of the shear strength of the contact surface with the development of plastic slip.

[0076] In this calculation logic, the dynamic hardening modulus The specific value is obtained by the processor by calculating the partial derivative of the dynamic friction coefficient evolution equation with respect to plastic displacement. The processor also calculates the consolidation strengthening factor. and shear damage factor The result is obtained by weighting the two factors relative to the rate of change of plastic displacement and the current effective normal stress. When the interface is in the consolidation-enhancing dominant phase, A positive value makes the stiffness matrix slightly smaller than the elastic stiffness; when the interface is in the shear damage-dominated stage, Negative values ​​result in a significant decrease in the stiffness matrix. The processor will calculate the... By assembling the global stiffness matrix into the overall finite element equation, the iterative oscillations caused by using the elastic stiffness matrix instead of the plastic stiffness matrix in traditional methods are eliminated, ensuring numerical convergence stability when simulating strong nonlinear shear softening and rapid loading process switching.

[0077] See attached document Figure 4 To visually represent the evolution of the pile-soil interaction mechanism during the consolidation of high embankments, the processor executes a post-processing algorithm based on shear stress sign reversal to automatically identify and track the location of the boundary point (i.e., the neutral point) between positive and negative skin friction on the pile side surface, and quantify the downward load effect of negative skin friction on the pile foundation.

[0078] Specifically, at each output time step, the processor iterates through the integration point data of all thickness-free contact surface elements on the pile-soil interface. The processor first extracts the shear stress components distributed along the pile axis. According to the notation convention of the present invention, when the soil settlement is greater than the pile settlement, a downward shear stress is generated at the interface, defined as negative skin friction (…). Conversely, when the pile settlement is greater than the soil settlement, an upward shear stress is generated, defined as positive skin friction (SFR). ).

[0079] The processor then performs precise location of the neutral point depth. Because the finite element mesh is discrete, the neutral point is not typically located exactly at the node position. The processor then retrieves the shear stress components. The sequence of two adjacent unit nodes that have changed. Let the depth be... The nodal shear stress at the location is Depth is The nodal shear stress at the location is And satisfy At this point, the processor uses a linear interpolation algorithm to calculate the instantaneous neutral point depth corresponding to the zero shear stress. The calculation formula is as follows: ; in, This indicates the absolute value operation. This interpolation process ensures the continuity of the neutral point position and eliminates the influence of grid size on the accuracy of boundary point identification.

[0080] Next, based on the determined neutral point depth, the processor calculates the total downward load borne by the pile. This calculation process only applies to areas located above the neutral point (i.e., at depth). Numerical integration is performed on the region of [missing information]. The processor maps the negative skin friction distribution onto the pile perimeter, and the integral summation formula is constructed as follows: ; in, This represents the total number of thickness-free contact surface elements located above the neutral point; This is the calculated diameter of the pile foundation; For the first The length of a thicknessless contact surface unit along the pile axis; This represents the axial shear stress at the integration point of the element. Represents pi; Indicates to From 1 to Summing up all the terms.

[0081] Finally, the processor calculates the instantaneous neutral point depth. Total pull-down load It is stored as a time-varying sequence, and combined with the aforementioned process state weighting factor. Perform correlation analysis. The system automatically generates a neutral point migration curve, which reveals the relationship between the rapid filling phase and the neutral point migration curve. During the intermittent consolidation period, the neutral point temporarily shifts upward due to the accumulation of excess pore pressure inhibiting soil settlement. ), with consolidation strengthening factor The dynamic law of rapid downward displacement of the neutral point due to increased and effective stress recovery and significant soil settlement is demonstrated. This spatiotemporal data reconstruction process transforms complex constitutive integral results into macroscopic evaluation indicators required for engineering design.

Claims

1. A numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations, characterized in that, include: A pile-soil flow coupled finite element model containing elements with no thickness contact surface is constructed, and state variables are initialized at the integration points of the elements with no thickness contact surface. At the start of the time increment step of the numerical calculation, the process state weight factor is calculated based on the global loading data, and the process state weight factor is substituted into the pile-soil-flow coupled finite element model. The local permeability coefficient of the contact surface of the pile-soil fluid coupling finite element model and the consolidation strengthening factor in the state variables are simultaneously corrected using the process state weight factor to simulate the regulatory effect of the process on fluid-structure coupling and strength growth. Based on the updated state variables, a dynamic friction coefficient is synthesized and a consistent tangent stiffness matrix is ​​constructed. The pile-soil flow coupled finite element model is then subjected to nonlinear iterative solution to obtain converged shear stress field data. Based on the shear stress field data, the depth location where the shear stress sign reverses is traced as the neutral point, and the negative friction distribution is calculated by integration.

2. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The state variables include cumulative plastic shear work, consolidation strengthening factor, shear damage factor, and process state weighting factor of the previous increment step. During the initialization phase, the cumulative plastic shear work is set to zero, the consolidation strengthening factor and the shear damage factor are set to unit values, and the process state weight factor of the previous incremental step is set to zero. The thickness-free contact surface unit is topologically composed of master face nodes belonging to the pile foundation surface and slave face nodes belonging to the soil surface.

3. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The specific steps for calculating the process state weighting factor include: Calculate the ratio of the effective load increment to the time increment within the current increment step, and divide the ratio by the preset reference pressure to obtain the dimensionless load change rate. A nonlinear smooth step function model containing a natural exponent term is constructed. The nonlinear smooth step function model includes shape control parameters, process judgment threshold, and transition region bandwidth parameters. The load change rate is mapped to the process state weight factor using the nonlinear smooth step function model, so that when the load change rate is greater than the process judgment threshold, the process state weight factor approaches 1, representing the undrained shear-dominated state. When the load change rate is less than the process determination threshold, the process state weight factor approaches 0, representing the dominant state of drainage consolidation.

4. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The steps for correcting the local permeability coefficient using the process state weighting factor specifically include: The local permeability coefficient at the current moment is calculated based on the multiplicative coupling rule. The local permeability coefficient is equal to the product of the initial permeability coefficient, the porosity correction term, and the process hindrance correction term. The void ratio correction term is defined as an exponential function of the volumetric strain of the thicknessless contact surface unit, used to characterize the influence of soil skeleton geometric deformation on the pore size of the seepage channel; The process stagnation correction term is defined as a linear function of the process state weighting factor and the preset process stagnation coefficient, such that when the process state weighting factor approaches 1, the value of the process stagnation correction term decreases, forcibly reducing the local permeability coefficient, in order to simulate transient superpore water pressure accumulation.

5. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The specific steps of performing damping updates on the consolidation strengthening factor using the process state weighting factor include: Read the excess pore water pressure values ​​and initial peak excess pore water pressure values ​​of adjacent soil elements in the pile-soil flow coupled finite element model, and calculate the local degree of consolidation. The theoretical target strengthening value is set as a power function of the local consolidation degree, which includes the maximum consolidation gain coefficient and the consolidation hardening exponent. The update magnitude is determined based on the process state weight factor. When the process state weight factor approaches 1, the change of the consolidation strengthening factor is suppressed, so that the consolidation strengthening factor remains at the value of the previous moment. When the process state weight factor approaches zero, the consolidation strengthening factor is controlled to approach the theoretical target strengthening value.

6. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 2, characterized in that, The numerical simulation method further includes the step of updating the shear damage factor based on the accumulated plastic shear work: Calculate the dot product of the shear stress vector and the plastic relative displacement increment vector within the current time increment step, and sum the results to the accumulated plastic shear work; The shear damage factor at the current moment is calculated using a negative exponential decay function, which includes residual damage ratio and damage rate parameters, such that the shear damage factor decays nonlinearly from a unit value to the residual damage ratio as the cumulative plastic shear work increases.

7. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The steps of synthesizing the dynamic friction coefficient and updating the uniform tangent stiffness matrix based on the updated state variables specifically include: The dynamic friction coefficient is defined as the product of the initial static friction coefficient, the consolidation strengthening factor, and the shear damage factor. Based on the yield criterion, if the state is plastic, a uniform tangent stiffness matrix containing a plasticity correction term is constructed. The uniform tangent stiffness matrix is ​​obtained by subtracting the plasticity projection component from the elastic stiffness matrix. The dynamic hardening modulus is calculated as the denominator of the plastic projection component. The dynamic hardening modulus is obtained by weighted summation of the partial derivatives of the consolidation strengthening factor and the shear damage factor with respect to the plastic displacement.

8. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The numerical simulation method further includes the steps of performing adaptive time-step cooperative control and long-term time-history calculation: The process state weighting factor is fed back to the time step controller of the global solver of the pile-soil-flow coupled finite element model; Based on the process state weight factor, a linear growth model for the next increment step time span is constructed. When the process state weight factor approaches 1, the growth of the time increment step is restricted to capture transient responses; when the process state weight factor approaches 0, the time increment step is allowed to be amplified to accelerate the consolidation calculation process. The co-calculation mechanism for secondary consolidation rheology and consolidation damping specifically includes: The additional viscous strain increment generated by secondary consolidation is calculated based on the logarithmic time rule and used to simulate the long-term creep of the soil skeleton under constant effective stress. A hydraulic and mechanical coupling damping term is superimposed during the effective stress update process. The hydraulic and mechanical coupling damping term is defined as a damping stress that is proportional to the time derivative of the pore water pressure and the time increment step, and is used to suppress numerical oscillations.

9. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 1, characterized in that, The step of using the depth location where the sign of shear stress reverses as the neutral point specifically includes: Extract the shear stress components distributed along the axial direction of the pile in the pile-soil flow coupled finite element model; Retrieve the sequence of two adjacent element nodes where the sign of the shear stress component changes; The depth corresponding to the zero shear stress component is calculated using a linear interpolation algorithm, and the depth is recorded as the instantaneous neutral point depth at the current moment.

10. The numerical simulation method for the spatiotemporal evolution of negative skin friction in high-fill pile foundations according to claim 9, characterized in that, The specific steps for integral calculation of negative skin friction distribution include: Traverse all thicknessless contact surface units located above the instantaneous neutral point depth; For each of the aforementioned contact surface elements without thickness, calculate the product of the pile foundation perimeter, the element axial length, and the absolute value of the axial shear stress at the element integration point. The total downward load borne by the pile is obtained by summing the products of all the contact surface units with no thickness located above the instantaneous neutral point depth.

Images