A numerical method for calculating the long-term settlement of three-dimensional subgrade under numerous traffic loads

Settlement problems caused by various traffic loads significantly impact infrastructure sustainability. Accurately predicting the long-term settlement of subgrade is crucial for improving infrastructure resilience. This study determined the element strain of soil using an explicit strain model, based on the bounding surface model and associated flow rules. Using the method of generalized shear strain equivalent, a cyclic cumulative factor is proposed for determining the magnitude of cumulative strain increment. The direction of plastic strain increment was determined using the bounding surface model combined with associated flow rules. A cumulative strain calculation method integrating the bounding surface and explicit strain models under general stress states was established. This computational method, combined with equivalent finite elements, can be used for calculating the long-term settlement of three-dimensional subgrade under numerous cyclic loads at general stress states through a user-defined subroutine provided by the commercial finite element modeling (FEM) software. The reliability of the method was verified through unit tests.


Introduction
Traffic loads, characterized by their low-frequency, long-term nature and numerous cyclic actions, cause settlement problems in natural soft subgrade.Under long-term cyclic loading, soil deformation gradually accumulates with an increase in the number of cycles.Consequently, developing calculation methods for estimating subgrade deformation under traffic cyclic loading is necessary to predict settlement of transportation infrastructure during operational lifespans.
Current models for simulating the cumulative strain of soil under cyclic loading are primarily categorized into two types.One involves establishing a complex constitutive model based on the elastic-plastic theory, while the other involves establishing an explicit strain model based on experimental or measured data.Constitutive models based on the elastic-plastic theory, such as the subloading yield surface [1], kinematic hardening [2] and bounding surface [3][4][5] models, have rigorous frameworks, reflecting the stress-strain relationship of the soil in each cycle under dynamic loading more realistically.These models typically use a small-step numerical integration method in the calculation process, leading to high computational cost and time.By contrast, explicit strain models avoid the complex iterative calculations of cyclic principal structures.In addition, they consider the effects of static and dynamic stress levels, number of cycles, and other soil factors.For example, Monismith et al. [6] established an exponential model of strain variation with the number of cycles.Huang and Yao [7] established an explicit strain model of axial cumulative plastic strain and pore pressure based on the normalized static and dynamic stress levels.
Ma et al. [8] established an equivalent finite element method.This method determines the strain components of soil elements based on the explicit model and the principal stress direction, realizing the combination of the finite element method and explicit model.This method, established based on the assumptions of elastic mechanics and isotropy, enables calculating the long-term settlement of subgrade in a two-dimensional plane strain model.However, according to the elastic-plastic theory, the direction of the plastic strain occurring in the soil is related to the plastic potential surface and flow rules.In addition, this method can only compute the element strain under two-dimensional conditions and cannot be applied to a three-dimensional model to calculate complex engineering situations.
This study presents a novel approach for calculating the cumulative strain components of soil under the generalized stress state of the explicit and bounding surface models.This calculation method can be combined with the basic principle of equivalent finite element method.By writing a custom constitutive program in commercial finite element software, it can be used to calculate the long-term settlement of threedimensional subgrade under numerous traffic loads in complex engineering situations.

Computational principles of equivalent finite element methods
Ma et al. [8] proposed the equivalent finite element method based on the principle of the initial strain nonlinear analysis method, adapting the initial strain method to realize finite element calculations using explicit strain models.Assuming that the initial strain generated by the soil element under traffic loads is {ߝ }, and the elastic stiffness matrix is ‫,]ܦ[‬ the corresponding initial stress {ߪ } can be obtained.Moreover, the equivalent load vector on the nodes of the element strains generated by the initial strain can be obtained ‫ܨ{‬ ఌ బ }, expressed as follows: According to the principle of superposition, the initial stresses introduced into the original equilibrium conditions cause node displacements within the cells.The resulting node displacement vectors are then solved using basic finite element equations.
A comparison between Equations (1) and (2) reveals that the selection of an arbitrary elastic stiffness matrix ‫]ܦ[‬ has the same effect on the equivalent load vectors ‫ܨ{‬ ఌ బ } and total stiffness matrix ‫]ܭ[‬ generated by the equivalent nodes on the unit.Therefore, any elastic stiffness matrix ‫]ܦ[‬ can be selected without affecting the results of the calculation of the nodal displacement vectors ‫ݑ{‬ } within the cell, as long as it falls within the range of numerical convergence.

Fundamental assumption
Without considering the creep properties of soil, the total strain of the soil under cyclic loading can be divided into recoverable elastic and irrecoverable plastic strain increments, as in the following equation: IOP Publishing doi:10.1088/1755-1315/1332/1/012010 The cumulative plastic strain of soil can be divided into two parts: plastic shear and plastic volumetric strains, expressed as follows:

Boundary surface function
Based on the bounding surface model established by Dafalias and Herrmann [3], the bounding surface function of the model is divided into three parts, consisting of two elliptical segments and one hyperbolic segment.Notably, for the stress-strain response of subgrade soil under cyclic traffic loading, the stress paths are contained within the first segment of the elliptical bounding surface, taking the form of an elliptical yield surface based on the modified Cambridge model: where ‫ݏ‬ ‾ and ‫‬ ‾ are the bias stress tensor and average stress corresponding to the virtual stress points mapped by the real stress points on the bounding surface, respectively, ‫ܯ‬ is the critical state stress ratio in the modified Cambridge model and is a material parameter, and ‫‬ is the preconsolidation pressure on the bounding surface, which is the intersection with the p-axis and controls the size of the bounding surface.

Mapping rule
The mapping rule uses the radius mapping rule, with the center of projection serving as the origin in the pq stress space.The stress state at the projection point can be expressed as follows: where ܾ is the distance between the real and virtual stress points on the boundary surface.

Flow rule
Using the associated flow rule to characterize the cyclic cumulative plastic strain, the plastic flow direction is determined by the normal direction at the virtual stress point on the boundary surface.The plastic strain increment can be expressed as follows: where Λ is the loading factor, ߪ is the real stress, and ߪ ‾ is the virtual stress on the corresponding boundary surface using the radial mapping rule.

Strain increment
Suppose that the true stress state of the soil when a dynamic load is applied to the natural subgrade is ߪ .A virtual stress point ߪ ‾ on the boundary plane can be uniquely determined by the mapping rule.The direction of the plastic strain increment can be determined based on the normal direction to this point, and the loading factor Λ remains constant for any stress state.Thus, the relative magnitude of the plastic strain component depends on . For the real stress state ߪ and corresponding virtual stress state ߪ ‾ , based on the radial mapping rule, Substituting this into the boundary surface function, the virtual stress point can be identified as follows: ) , it follows that Consequently, ߪ ‾ = ܾߪ and பி பఙ ‾ ೕ can be determined.Subsequently, b is substituted into Equation (9).
Because p in the bounding surface function is a bounding surface model parameter, its value has no effect on the relative magnitude of each strain component.The above equation is primarily used to determine the direction of the increment of the plastic strain.Thus, it can be considered as constant.Therefore, its influence is disregarded, yielding where ߣ is the cyclic cumulative strain factor, characterizing the magnitude of the cumulative strain increment under each cycle.ߣ remains constant for a given stress state.
The generalized shear strain (or strain strength) reflects the deformation strength of the soil, and the incremental generalized shear strain can be expressed as follows: where ‫ܬ‬ ᇱ ଶ is the second invariant of the strain deviation.
The above equation reveals that for a given stress state ߪ , the generalized shear strain is solely dependent on the cyclic cumulative strain factor ߣ.Each strain component increment can be obtained and expressed by the following equation: where ߪ is the true stress state, ߣ is the cyclic cumulative strain factor, and M is the critical state stress ratio.

Explicit strain model for determining the generalized shear strain
In this study, the method of determining the generalized shear strain increment using an explicit model is illustrated based on the explicit model established by Huang and Yao [7].The expression of axial cumulative plastic strain is as follows: where ‫ݍ‬ ௨௧ is the undrained shear strength of saturated soft clay.Liu et al. [9] derived the following formula for the shear strength of ‫ܭ‬ -consolidated saturated soft clay based on the boundary surface and critical state theories: where ‫ܯ‬ is the critical state stress ratio, ‫‬ is the initial average effective circumferential pressure, ߣ and ߢ are the slopes of the normal consolidation and rebound lines in the space, respectively, and ߙ is the inclination angle of the yield surface in the ‫‬ − ‫ݍ‬ space, calculated using the following equation: where A is the correction factor.The model is a function of the number of cycles ܰ when the initial stress state, dynamic stress level, and relevant calculation parameters are known.Thus, the above equation is rewritten in incremental form as follows: This yields where ݀ߝ is the axial cumulative plastic strain, ߰൫ߪ , ܰ, ߦ൯ is the explicit model for the axial plastic strain increment established based on the results of the indoor test.This model is influenced by the soil dynamic load level, initial stress state, number of cycles, and other model parameters.The cyclic cumulative strain factor ߣ can be determined using the following equation: By substituting the cyclic cumulative strain factor ߣ into Equation ( 12), ቄ݀ߝ ቅ can be obtained.

Cumulative pore pressure dissipation corresponding to element strain
The results of indoor undrained triaxial tests reveal that the subgrade under cyclic loading generates cumulative pore pressures.The cumulative pore pressures under numerous cyclic dynamic loads can be expressed in the incremental form by the explicit strain model established by Yao et al. [10]: The volumetric strain ߝ ௩ due to the complete dissipation of the cumulative pore pressure considering long-term effects can be obtained from the e-lnp curve with the following expression: where ‫‬ is the preconsolidation pressure, and ݁ is the initial pore ratio.The dissipation of the cumulative pore pressure causes an increase in the effective mean body stress of the soil.The cumulative strain increment can be calculated using the following equation: where ‫ݑ‬ * is the cumulative pore pressure, obtained using the explicit strain model, and ܷ is the degree of consolidation.Considering long-term action, the cumulative pore pressure can be considered to be not fully dissipated.The degree of consolidation ܷ takes the value of 1.Therefore, the cumulative pore pressure dissipation caused by the element strain is given by From Equations ( 12) and ( 24), the plastic shear and plastic volumetric strain increments can be obtained, and the unit strain increment is obtained by adding the two parts of the strain increment.

Long-term settlement prediction of subgrade under numerous dynamic loads
The long-term settlement of the three-dimensional subgrade model can be predicted using the proposed soil element strain calculation method consider numerous cyclic dynamic loads combined with the basic principles of the equivalent finite element method and the UMAT subroutine in ABAQUS.Figure 1 shows the realization of the process.

Verification program
To verify the rationality of the proposed method and the correctness of the UMAT subroutine writing, undrained cyclic triaxial tests were simulated on soil units under deviatoric stress consolidation.This element verification did not consider the volumetric strain due to pore pressure dissipation, but only the cumulative shear strain of the soil under cyclic loading.The calculated results were compared with the explicit strain model predictions of Huang and Yao [7].A soil sample with dimensions of 40 mm × 40 mm × 40 mm was modeled using a first-order 8-node 3D solid unit (C3D8).Table 1 summarizes the simulation scheme for analyzing the deformation pattern of soil under different consolidation perimeter pressures and dynamic stress levels for deviatoric consolidation.

Simulation results
The results of the numerical simulation are compared with the results of the explicit strain model proposed by Huang and Yao [7], as shown in figure 2. The model parameters were kept consistent with those in the literature, with the specific parameter values shown in Table 2.The figure reveals that the results of the proposed strain calculation method are consistent with the predicted values using the explicit strain model under different consolidation confining pressures and dynamic stress levels.Notably, the rate of strain development was faster at smaller cycle counts.As the number of cycles increased, the rate of strain accumulation slowed down.Moreover, the simulation results under different consolidation pressures and dynamic stress ratios align with those obtained using the explicit strain model.The values were slightly smaller than those predicted using the explicit strain model, with a difference of approximately 5-10 %.Under high dynamic stresses, the difference further reduced to approximately 3%.This error arises because the proposed computational method is based on a joint computation of the explicit strain and bounding surface models.Nonetheless, the resulting error falls within acceptable limits.This indicates that the proposed method of calculating the soil element strain under numerous cyclic dynamic loads could better describe the cyclic cumulative deformation characteristics of soil samples under different consolidation circumferential pressures and dynamic stress ratios, exhibiting high numerical calculation accuracy.Therefore, applying the method in engineering cases can achieve better simulation results.

Results
The proposed long-term settlement prediction method for subgrade under cyclic dynamic loads is based on the general stress state.The overall subgrade settlement is obtained using model stress data combined with the explicit strain model, facilitating the analysis of the location and distribution of the occurrence of inhomogeneous settlement.The associated flow rules and boundary surface function consider the direction of the plastic strain increment of the soil, reflecting the plastic strain characteristics of the soil under cyclic loading.The cell validation results demonstrated that the finite element calculation results using the proposed method are consistent with the predictions made using the explicit strain model.The simulation results for different consolidation pressures and dynamic stress ratios align with those predicted using the explicit strain model.

7 Figure 1 .
Figure 1.Flowchart of the method for predicting the long-term settlement of subgrade under numerous cyclic dynamic loads.