The schematic diagram illustrating the calculation of soil consolidation under continuous drainage boundary conditions and cyclic loading is presented in Figure 1, where represents the cyclic load, which varies over time; H denotes the thickness of the soil layer; and z indicates the depth of the soil layer.

(1) The stress-strain relationship of the soil is assumed to be hyperbolic as shown in Equation 1 (Zhang and Sun, 2007):

Based on these basic assumptions, the governing equation for the consolidation of the hyperbolic model can be derived as shown in Equation 6: C v E 0 1 E 0 + n σ ′ 2 ∂ 2 u ∂ z 2 − 2 n E 0 + n σ ′ 3 ∂ σ ′ ∂ z ∂ u ∂ z = − E 0 E 0 + n σ ′ 2 ∂ σ ′ ∂ t (6)

It follows from the equality of the total stresses (Shi et al., 2001): C v ∂ 2 σ ′ ∂ z 2 − 2 n E 0 + n σ ′ ∂ σ ′ ∂ z 2 = ∂ σ ′ ∂ t (7)

Continuous drainage boundary conditions are introduced into Equation 7 to obtain the initial and boundary conditions as follows:

Initial condition as shown in Equation 8: t = 0 , u = p 0 , σ ′ = σ 0 ′ (8)

Top surface boundary conditions as shown in Equation 9: z = 0 , u = e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ , σ ′ = p t − e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ + σ 0 ′ (9)

Bottom boundary conditions: z = H , u = e − c t p 0 + ∫ 0 t e c τ p ′ τ d τ , σ ′ = p t − e − c t p 0 + ∫ 0 t e c τ p ′ τ d τ + σ 0 ′ (10) where b and c in Equation 10 are boundary permeability parameters which can be obtained through test fitting or engineering measurement. When the boundary condition is set to a completely permeable state, b and c approach infinity; conversely, when the boundary condition is impermeable, b and c equal zero.

Equation 7 is solved via variable substitution: ξ = 1 + e 0 σ ′ E 0 + n σ ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t (11)

Substituting ∂ σ ′ ∂ t , ∂ σ ′ ∂ z and ∂ 2 σ ′ ∂ z 2 into Equation 7 yields Equation 11 (Davis and Raymond, 1965): ∂ ξ ∂ t = C v ∂ 2 ξ ∂ z 2 − J t (12) where J t = E 0 1 + e 0 p ′ t E 0 + n p t + n σ 0 ′ 2

The initial and boundary conditions for Equation 12 are expressed as:

Initial condition as shown in Equation 13: ξ z , 0 = 1 + e 0 σ 0 ′ E 0 + n σ 0 ′ − 1 + e 0 σ 0 ′ + p 0 E 0 + n σ 0 ′ + p 0 (13)

Top surface boundary condition: ξ 0 , t = 1 + e 0 p t − e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ + σ 0 ′ E 0 + n p t − e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ + σ 0 ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t (14)

Bottom boundary condition: ξ H , t = 1 + e 0 p t − e − c t p 0 + ∫ 0 t e c τ p ′ τ d τ + σ 0 ′ E 0 + n p t − e − c t p 0 + ∫ 0 t e c τ p ′ τ d τ + σ 0 ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t (15)

The non-homogeneous boundary is homogenized as follows: ξ z , t = V z , t + χ t , z χ t , z = A t · z + B t (16)

Substituting Equations 14, 15 into Equation 16 yields: ξ z , t = V z , t + z H 1 + e 0 p t − e − c t ( p 0 + ∫ 0 t e c τ p ′ τ d τ ) + σ 0 ′ E 0 + n p t − e − c t ( p 0 + ∫ 0 t e c τ p ′ τ d τ ) + σ 0 ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t − 1 + e 0 p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ E 0 + n p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ + 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t + 1 + e 0 p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ E 0 + n p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t (17)

Next, we calculate the first-order derivative of with respect to t and the second-order derivative with respect to z in the above equation and substitute these results to Equation 12 to obtain: C v ∂ 2 V ∂ z 2 = ∂ V ∂ t + G t (18) where G t = z H E 0 1 + e 0 c e − c t p 0 + ∫ 0 t e c τ p ′ τ d τ E 0 + n p t − e − c t p 0 + ∫ 0 t e c τ p ′ τ d τ + σ 0 ′ 2 − E 0 1 + e 0 b e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ E 0 + n p t − e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ + σ 0 ′ 2 + E 0 1 + e 0 b e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ E 0 + n p t − e − b t p 0 + ∫ 0 t e b τ p ′ τ d τ + σ 0 ′ 2

At this point, the homogeneous boundary conditions are shown in Equation 19: z = 0 , V 0 , t = 0 z = H , V H , t = 0 t = 0 , V z , 0 = 0 (19)

The variable separation method is utilized to solve Equation 18, and the solution is expressed as: V z , t = ∑ m = 1 , 3 , 5 … ∞ f m t sin m π z H (20)

Substituting Equation 20 into Equation 18 and simplifying gives in Equation 21: f m ′ t + C v m π H 2 f m t = − 4 G t m π (21)

To solve the first-order nonhomogeneous differential equation of the above formula, we obtain: f m t = − 4 m π · e − C v m π H 2 t ∫ 0 t G x e C v m π H 2 x d x (22)

In accordance with the description in Chen (Chen, 2002), the method of separation of variables is employed. Substituting Equation 20 in Equation 22 yields the solution of Equation 18 as: V z , t = − ∑ m = 1 , 3 , 5 … ∞ 4 m π e − C v m π H 2 t ∫ 0 t G x e C v m π H 2 x d x sin m π z H (23)

Substituting Equation 23 into Equation 17 yields: ξ z , t = − ∑ m = 1 ∞ 4 m π · e − C v m π H 2 t ∫ 0 t G x e C v m π H 2 x d x · ⁡ sin m π z H + z H 1 + e 0 p t − e − c t ( p 0 + ∫ 0 t e c τ p ′ τ d τ ) + σ 0 ′ E 0 + n p t − e − c t ( p 0 + ∫ 0 t e c τ p ′ τ d τ ) + σ 0 ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t − 1 + e 0 p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ E 0 + n p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ + 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t + 1 + e 0 p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ E 0 + n p t − e − b t ( p 0 + ∫ 0 t e b τ p ′ τ d τ ) + σ 0 ′ − 1 + e 0 σ 0 ′ + p t E 0 + n σ 0 ′ + p t (24)

The analytical solution for the effective stress which is shown in Equation 25 is obtained by combining Equation 24 with Equation 11: σ ′ z , t = E 0 ξ z , t E 0 + n σ 0 ′ + p t + E 0 1 + e 0 σ 0 ′ + p t E 0 1 + e 0 − n ξ z , t E 0 + n σ 0 ′ + p t (25)

The increment of the foundation relative to the initial effective stress is expressed as Equation 26: Δ σ ′ z , t = σ ′ z , t − σ 0 ′ (26)

The analytical solution for the consolidation settlement is given in Equation 27: w t = ∫ 0 H σ ′ z , t E 0 + n σ ′ z , t dz − ∫ 0 H σ 0 ′ E 0 + n σ 0 ′ dz (27)

To gain deeper insight into the consolidation settlement behavior of the soil within this specific model, a comprehensive parametric analysis was meticulously conducted. The purpose of this investigation was to elucidate the intricate relationships between soil properties and the resultant settlement patterns, thereby enhancing the understanding of soil behaviors under varying conditions of stress and compaction. By systematically varying key soil parameters such as permeability, compressibility, and shear strength, we aimed to uncover the underlying mechanisms governing soil consolidation and settlement. Dynamic loads, which may vary over time due to factors such as traffic, wind, or seismic activity, impose additional stresses that can lead to increased soil deformation and potential failure. Figure 4 provides a detailed depiction of consolidation settlement variations over time under diverse boundary conditions, characterized by the parameters b and c. These parameters represent the boundary permeability quality of the soil and can be used to simulate the permeability behaviors of the corresponding consolidated soil in engineering practice. Specifically, larger b and c values corresponded to enhanced boundary permeability, facilitating the efficient drainage of pore water and thus accelerating the consolidation process. Conversely, smaller values denoted reduced permeability, which hindered drainage and subsequently slows down settlement. The figure clearly illustrated that, under sinusoidal loading, the amplitude of oscillation in the consolidation settlement curve is positively correlated with the model parameters b and c. This correlation highlighted the influence of boundary permeability on the soil’s dynamic response to cyclic loading. As permeability increased, the soil’s capacity to dissipate pore water pressure under load also increased, leading to a more pronounced oscillation in settlement. Furthermore, the settlement variation curve over time, similar to the effective stress, exhibited a trend of oscillatory incrementation. This oscillatory behavior was a manifestation of the soil’s viscoelastic behavior and its capacity to adapt to the applied loads over time. The cyclic nature of the load induced a corresponding cyclic response in settlement, which was superimposed on the overall trend of consolidation settlement. Moreover, as the permeability of the soil boundary increased, evidenced by larger values of parameters b and c, there was a corresponding escalation in the growth rate of settlement curve. Conversely, when the permeability of the boundary decreased, represented by smaller values of b and c, the settlement hysteresis effect became more pronounced. This effect arose from the restricted drainage conditions that led to a slower consolidation process and a more significant lag in settlement response to load changes, as observed in the soil behavior under cyclic loading. In addition, the similarity between the two consolidation settlement curves, b = 0.2 d − 1 , c = 0.3 d − 1 and b = 2 d − 1 , c = 3 d − 1 , was striking. This observation highlighted the superior convergence of the model’s computational outcomes, which demonstrated that it can be extensively applied to engineering practice.

When parameters b and c were assigned varying values, Figure 5 illustrates the resultant changes in the effective stress along the depth of the soil. As depicted in Figure 5, influenced by the permeability of the side boundaries, variation curve of the effective stress initially decreased and then increased with increasing soil depth. This behavior was a direct consequence of the interplay between fluid flow and stress distribution within the soil matrix. The fluid flow within the soil is influenced by its permeability, which is a measure of its ability to transmit fluids. The permeability is a key factor in determining how quickly water can drain from the soil, thereby affecting the rate at which effective stress is established. As water drains from the soil pores, the effective stress increases, leading to soil consolidation and settlement. When the permeability is high, indicated by larger values of b and c, the soil readily expelled pore water, resulting in a swift decrease in pore water pressure and a corresponding increase in effective stress near the surface. This initial decrease in the effective stress curve reflected the rapid dissipation of excess pore water pressure due to the high permeability conditions, which facilitated a more rapid consolidation process. As depth increased, the effective stress subsequently increased due to the overburden pressure exerted by the overlying soil layers. This resulted in a characteristic S-shaped curve, where the inflection point corresponded to the depth at which the influence of boundary permeability on stress distribution became less pronounced. In contrast, when the permeability was low, represented by smaller values of b and c, the pore water drainage was impeded, leading to a more gradual decrease in pore water pressure and a delayed increase in effective stress. The reduced permeability resulted in a more pronounced hysteresis effect in the subsidence response to load changes, as the soil’s ability to consolidate was constrained by the limited capacity to expel water from the soil matrix. This delayed consolidation was particularly evident in the shallower portions of the soil profile, where the effective stress curve exhibited a less steep initial decrease, followed by a more gradual increase with depth. As depicted in Figure 5, when b and c took different values, the distribution of effective stress with soil depth exhibited a non-monotonic pattern, which was in contrast to the scenario under constant loading conditions.

Figure 6 visually represents the impact of the initial void ratio e 0 , which describes the volume of voids relative to the volume of soil solids, on the soil consolidation settlement. The analysis revealed that variations in the initial void ratio e 0 have a negligible influence on the consolidation settlement. This minimal effect was primarily attributed to the slight changes observed in 1 + e 0 , which was related to the soil’s physical properties and had a minimal impact within the governing equation of soil behavior.

Figure 7 illustrates the impact of varying initial compression modulus on the consolidation settlement within the hyperbolic model. As the initial compression modulus gradually increased from 0.5 MPa to 5 MPa, there was a significant decrease in the rate of consolidation settlement of the soil. This trend was evident because the initial compression modulus that reflects the soil’s inherent stiffness was directly related to its resistance to compression under initial loading conditions. A higher initial compression modulus indicated a stiffer soil, which resulted in a lower settlement rate. The cumulative decrease in settlement became more pronounced in the later stages of consolidation. This observation was attributed to the increased stiffness of the soil at higher compression module, which limited further compression as consolidation progressed. Additionally, the oscillation caused by cyclic load was smaller. Conversely, a lower modulus suggested a softer, more compressible soil that settled more rapidly under the same load.

As depicted in Figure 8, the variation in the slope n of the compression curve corresponded to distinct changes in soil subsidence behavior. Here, n represents the rate of modulus growth during the loading process. A smaller value of n indicates a lower rate of modulus growth, while a larger value of n signifies a higher rate of modulus growth, which indicates that the soil’s stiffness increases rapidly under load. This rapid increase could lead to a more abrupt response to loading and a larger settlement rate during the early stages of loading. Therefore, as n increased, the rate of soil consolidation slowed down and the fluctuation of the curve became less pronounced.