An isochronous stress ratio logarithmic strain curve based clay creep model considering the effects of hardening and damage

Creep is one of the typical mechanical properties of clay, and studying the creep mechanical properties of clay is of great significance to construction projects in clay sites. This study conducted creep tests on Chengdu clay and found that the soil mass underwent elastic deformation, decay creep deformation, steady-state creep deformation, and accelerated creep deformation. The isochronous stress ratio-logarithmic strain curves and their mathematical models were proposed to thoroughly analyze clay creep mechanical properties. Creep automatic feature points, such as linear elastic extreme point, initial yield point, long-term strength point, and plastic point, were identified on the curve. Considering the hardening and damage effects during creep loading, linear elastic and viscoelastic elements considering the time-dependent damage, a viscoplastic element considering the load hardening effect, and viscoplastic and plastic elements considering the load damage effect were established based on the element model and the Riemann–Liouville fractional derivative. Based on the mechanical properties of the whole clay creep process, the creep mechanical feature points, and the established element model, a clay creep model was proposed considering the hardening and damage effects. The rationality and regularity of the creep model were verified using the creep test data. This research accurately revealed the creep mechanical properties of clay and facilitated soil deformation prediction, thus providing technical guidance and references for construction projects in clay sites.


Test methods
The creep tests were conducted on a GDS triaxial test system.The specimens were prepared as follows: (1) The specimens were 38 mm in diameter and 76 mm in height, and apply petroleum jelly to the gasket in contact with the sample; (2) The specimens were consolidated under 100 kPa confining pressure, which was loaded at a rate of 1 kPa/min, and the consolidation condition was that the pore water pressure dissipation of over 95%; (3) Creep tests were carried out using the stepwise loading method, and the deviator stress path was q = 20 kPa, 38 kPa, 55 kPa, 75 kPa, 90 kPa, 110 kPa, 128 kPa, and 145 kPa, with the deviator stress loading rate of 1 kPa/min.Figure 2 shows the Chengdu clay creep test.

Experimental result analysis
The creep test curve of Chengdu clay is shown in Fig. 3.It can be observed that when q ≤ 75 kPa, the clay creep strain is less than 1.2%.At this time, the strain includes the transient strain and the creep strain that gradually stabilizes, and the clay is in the decay creep stage.When q ≥ 145 kPa, the clay creep deformation increment accelerates until the failure of the specimen, at which time the clay is in the accelerated creep stage.When 75 kPa < q < 145 kPa, the clay gradually transitions from decay creep to accelerated creep, covering the steady-state creep stage.Therefore, the strain of Chengdu clay under creep test conditions mainly included transient strain and creep strain, and the creep process included decay creep, steady-state creep, and accelerated creep stages.

Stress-strain relationship under isochronous conditions Isochronous stress-strain curve
According to the creep test curves, the isochronous stress-strain curves are plotted in Fig. 4. It can be observed that when q ≤ 75 kPa, the isochronous stress-strain curves are roughly linear; When q > 75 kPa, the isochronous stress-strain curves gradually deflect to the strain axis, showing significant nonlinear characteristics; Meanwhile, the nonlinearity of the curves is more significant over time.The value at the inflection point of the isochronous stress-strain curve is generally taken as the long-term strength 19 .Significant inflection points can be observed on the isochronous stress-strain curves in Fig. 4 at q = 75 kPa and q = 110 kPa.Therefore, the long-term strength q L determined based on the isochronous stress-strain curve is 75 kPa to 110 kPa, which is 51.72% to 75.86% of the maximum loading stress in the creep tests, with a range of about 24%.Thus, the wide range of long-term strength determined based on the isochronous stress-strain curve makes it difficult to determine the long-term strength accurately.

Isochronous stress ratio-logarithmic strain curve
To accurately determine the creep mechanical characteristic values of clay, such as the long-term strength, the patterns of the isochronous stress-strain curve must be thoroughly analyzed.To this end, the stress and strain under isochronous conditions are transformed, as expressed in Eq. (1): where ε 1 and q are the axial strain and deviator stress on the isochronous stress-strain curve, respectively; q f is the creep failure load; q r and ε q are the stress ratio and logarithmic strain under isochronous conditions, respectively, with 0 < q r ≤ 1.
The isochronous stress-strain curve is converted into the isochronous stress ratio-logarithmic strain curve according to Eq. (1), and the converted curve is placed into the half-logarithmic coordinate system space, as shown in Fig. 5.It can be observed that the isochronous stress ratio-logarithmic strain curve has the following properties: (1) The curve has an inflection point C, and the ABC segment of the curve is concave, while the CDE segment is convex; (2) The upper boundary of the curve is q r = 1, and the lower boundary is q r = 0.With lim q r →1 ε q = 0 , ε 1 = +∞ , indicating that the soil mass reaches a state of complete plastic deformation or failure at this time; When lim q r →0 ε q = +∞,ε 1 =0 , indicating that the soil mass is in a fully elastic state at this time.
(1) q r = q q f ε q = − ln(ε 1 /100)  www.nature.com/scientificreports/According to the properties of the isochronous stress ratio-logarithmic strain curve, the tangent to the curve at point C intersects with the upper and lower boundaries at points M and N, respectively.The perpendiculars through points N and M intersect with the curve at points D and B, respectively, and the tangent at point D intersects segment MN through point F. According to the properties of the isochronous stress ratio-logarithmic strain curve and the mechanical properties of the whole creep process of the soil mass, points B, C, D, and F on the curve are defined as the creep mechanical feature points.Among them, point B is the linear elastic extreme point, point C is the initial yield point, point D is the complete yield point, and point F is the long-term strength point.Therefore, the AB segment on the isochronous stress ratio-logarithmic strain curve is a linear elastic stage describing the elastic deformation; The BC segment on the curve is the viscoelastic stage mainly representing the viscoelastic deformation in the decay creep stage; The CD segment on the curve is a viscoplastic stage primarily representing the viscoplastic deformation in the decay creep and steady-state creep stages; The DE segment on the curve is the plastic deformation stage mainly representing the deformation in the accelerated creep stage.
Isochronous stress ratio-logarithmic strain curve model.According to phenomenology, the isochronous stress ratio-logarithmic strain curve is highly similar to the soil-water characteristic curve.Therefore, the VG model of the soil-water characteristic curve is used as the mathematical model of the isochronous stress ratio-logarithmic strain curve, as expressed in Eq. ( 2): where k, n, and m are the model fitting parameters, m = 1-1/n.
To solve each characteristic point of the isochronous stress ratio-logarithmic strain curve with Eq. ( 2), the first-and second-order derivatives of Eq. ( 2) are derived, respectively, and the results are expressed in Eqs. ( 3) and ( 4).
At the inflection point C, we have q ′′ rC = 0 , and the coordinates of point C can be obtained from Eq. ( 4), as expressed in Eq. ( 5).
The slope K C of the tangent at point C can be derived from Eqs. ( 5) and (3), as expressed in Eq. ( 6).
The equation for the curve's tangent at point C can be derived from Eqs. ( 5) and ( 6), as expressed in Eq. (7).
The abscissa of point M and the coordinates of point D can be obtained by substituting q r = 1 into Eq.( 7), as expressed in Eq. (8). where . The abscissa of point N and the ordinates of point B can be obtained by substituting q r = 0 into Eq.( 7), as expressed in Eq. ( 9).

Chengdu clay creep test data analysis
The isochronous stress-strain curves in Fig. 4 are converted and analyzed according to Eqs. (1) to (11).The fitting results of k, n, and m are presented in Fig. 6.Limited by the length of the article, the isochronous stress ratiologarithmic strain curves and each creep mechanical feature points at t = 5 min and t = 1200 min are presented in Fig. 7.It can be observed from Fig. 6 that the parameter k increases linearly first with time and then increases nonlinearly; The parameter n decreases linearly at a larger rate first with time and then increases linearly at a lower rate.As shown in Fig. 7, the isochronous stress ratio-logarithmic strain test results have a fitting coefficient over 0.98 with their model, indicating a high degree of agreement between the test results and the model.Meanwhile, the calculation results of each mechanical feature point are clear.
According to Fig. 7, the elastic limit stress ratio q rB is between 0.1316 and 0.1322, and the deviator stress is about 19 kPa.This study selects the creep data under the loading level of q = 20 kPa to calculate the elastic modulus of the clay in the linear elastic stage, and the results are shown in Fig. 8.It can be observed that with the increase of time, the linear elastic modulus drastically decreases first and then slowly decreases until stabilization.In this study, the relationship between the linear elastic modulus and time fitted using Eq. ( 12).The fitting results are shown in Fig. 8.According to the fitting results in Fig. 8, the long-term elastic modulus of Chengdu clay E ∞ is 14.9578 MPa in the elastic deformation stage, E 0 = 39.5227MPa, and the parameter w is 0.005.( 11) where E 0 and E ∞ are the initial and long-term elastic modulus, respectively; w is the fitting parameter.

Riemann-Liouville fractional derivative
Fractional calculus is an important method to study viscoelastic materials.This study proposes to establish creep model elements using the Riemann-Liouville fractional derivative.According to the Riemann-Liouville theory 18,19 , we have the following.Let the function u be continuously integrable on (0, + ∞), and we have the following by finding its fractional integral over t ≥ 0 and Re(n) ≥ 0: where Ŵ(n) = ∞ 0 t n−1 e −t dt is the Gamma function, and α is the fractional order.The Riemann-Liouville fractional differential can be defined as follows: Let u ∈ C , γ > 0 , and m be the smallest integer greater than γ, we have the following by letting γ = β-α: If u(t) is integrable in the vicinity of t = 0, and 0 ≤ n ≤ 1, the Laplace transform of the fractional calculus is: where u(p) is the Laplace transform of u(t).

Elastic element considering time-dependent damage
The elastic element is used to describe the deformation pattern of the clay in the elastic stage.The linear elastic element is expressed in Eq. ( 16).
where σ is the loading deviator stress, σ = q, and ε e is the linear elastic strain.
According to Eqs. ( 12) and ( 16), the linear elastic element model considering the time-dependent damage can be expressed in Eq. (17).

Fractional viscoelastic element considering initial threshold
With σ e < σ ≤ σ y0 (σ e is the extreme value of linear elastic stress, and σ y0 is the initial yield strength), the clay has not yet yielded and is in the decay creep state.At this time, the viscoelastic element starts to work.The viscoelastic www.nature.com/scientificreports/element consists of a linear elastic element and a fractional viscous element in parallel, and its model is expressed in Eq. (18).
where σ ve is the stress of the viscoelastic element; σ e is the extreme value of the linear elastic stress, which can be determined by the feature point B on the isochronous stress ratio-logarithmic curve; ε ve is the strain of the viscoelastic element; E 1 and η 1 are the elastic modulus and viscosity coefficient of the viscoelastic element, respectively; α 1 is the fractional operator of the viscoelastic element.Through the Laplace transform, the viscoelastic element model can be obtained from Eq. ( 18), as expressed in Eq. (19). where

Viscoplastic elements considering initial yield strength hardening effect
With σ y0 < σ ≤ σ L (σ L is the long-term strength), the soil mass enters the yield stage.Yet, the soil mass deformation eventually stabilizes, and the structure rebalances.Therefore, the soil structure exhibits load hardening, and the soil mass is in the decay creep stage at this time.
In the initial yield-load hardening stage of the soil mass, the viscoplastic element considering the load hardening effect consists of a fractional viscous element and a plastic element in parallel, and its model is expressed in Eq. (20).
where σ vpH and ε veH are the stress and strain of the viscoplastic element considering the hardening effect, η 1 is the viscosity coefficient of the viscoplastic element, and α 1 is the fractional operator.σ y (t) is the yield strength of the soil mass after hardening over time.
where σ y0 is the initial yield strength, which can be determined by point C on the isochronous stress ratiologarithmic curve; H and λ are the material parameters, and 0 ≤ H < 1, 0 < λ < 1.The viscoplastic element model considering the load hardening effect can be derived from Eqs. ( 20) and ( 21), as expressed in Eq. (22).

Viscoplastic element considering soil mass strength damage effect
With σ L < σ ≤ σ a (σ a is the creep failure stress), the soil mass deformation after yielding increases at a steady rate, while the soil structure sustains damage until its failure.At this time, the soil mass is in a steady-state creep stage.
Kachanov's definition of the creep damage variable is expressed in Eq. ( 23).
where D is the damage variable of the nonlinear element, and A and β 1 are the material parameters.
The damage variable can be obtained by finding the integration of Eq. ( 23), as expressed in Eq. ( 24).
where t f and t L are the start times of accelerated creep and steady-state creep, respectively, which can be determined by points D and F on the isochronous stress ratio-logarithmic strain curve, with Based on the Lemaitre strain equivalence principle, the relationship between engineering stress and effective stress can be expressed in Eq. (25).
where σ and σ′ are the engineering stress and the effective stress, respectively.
In this study, the proposed steady-state creep-stage nonlinear damage element model is expressed in Eq. ( 26). ( 18) where ε vpD is the damaging strain and η 3 is the viscosity coefficient at the steady-state creep stage.
The nonlinear damage element model for the steady-state creep stage can be obtained from Eqs. ( 24) to ( 26), as expressed in Eq. (27).

Nonlinear accelerated creep element
With σ a < σ, the soil mass exhibits accelerated deformation until its failure and is in an accelerated creep stage.The accelerated creep pattern is described with a fractional viscous element, the model of which is expressed in Eq. (28).
where ε a , α 3 , and η 3 are the strain, fractional operator, and viscosity coefficient of the nonlinear accelerated creep element, respectively, with 1 < α 3 .σ a is the creep failure deviator stress, which can be determined by point D on the isochronous stress ratio-logarithmic strain curve.
Finally, the creep mechanics interval of each element model established on the isochronous stress ratiologarithmic strain curve is shown in Fig. 9.

Creep model
According to the creep properties of clay and the established element model, the clay creep model proposed in this study is shown in Fig. 10.The clay creep model is as follows.
(1) With 0 < σ ≤ σ e , only the elastic element undergoes deformation at this time, and the strain model is expressed in Eq. ( 29).
(4) With σ L < σ ≤ σ a , the viscoplastic element considering the damage effect begins to deform, and the creep model is expressed in Eq. ( 32).
(5) With σ a < σ, the nonlinear accelerated creep element begins to deform, and the creep model is expressed in Eq. (33).

Model validation
The proposed model in this study is verified using the Chengdu clay creep test data.Based on Figs. 6 and 7, the stress calculation results of the isochronous stress ratio-logarithmic strain curve feature points are expressed in Eq. (34), where t L = 8640 min, and t a = 11,520 min.
The creep model proposed in this paper is used to calculate the creep stages of Chengdu clay.The calculation parameters are listed in Table 2, and the comparison between the model calculation results and the creep test results is shown in Fig. 11.
It can be observed that in the elastic deformation stage, the model calculation results are basically consistent with the test results.Between the loading time t interval of 700 min to 1300 min, the calculated strain is greater than the tested strain, and the maximum difference is 0.016%.In each deformation stage after the elastic deformation stage, the model calculation results are highly consistent with the creep test results, which verifies the rationality and applicability of the proposed creep model.( 29)

Discussion
Based on feature points C and F on the isochronous stress ratio-logarithmic strain curve, the effects of curve model parameters k and n on the curve feature points are analyzed, as shown in Fig. 12.It can be observed that under the influence of the curve model parameter k, the curve feature points C and F move parallel to the strain axis as k increases, indicating that the parameter k only affects the strain at the feature point but not the stress at the feature point.Under the influence of the curve model parameter n, the strain at point C increases while the stress decreases as the parameter n increases-meanwhile, the strain and stress at point F nonlinear decreases.Therefore, the isochronous stress ratio-logarithmic strain curve is more sensitive to model parameters n.Based on the isochronous stress-strain curve, an isochronous stress ratio-logarithmic strain curve, which had significant concave-convexity and boundary properties, was proposed.The mathematical model of the isochronous stress ratio-logarithmic strain curve is proposed from the phenomenological perspective, and the calculation methods for the various creep mechanical feature points were provided, including the linear elastic extremum point, initial yield point, long-term strength point, and plastic point.2. Based on the isochronous stress ratio-logarithmic strain curve properties and feature points, the mechanical properties of the clay loading process were considered, and linear elastic, viscoelastic, viscoplastic, and plastic elements were established.In particular, when the loading stress of the clay was between the initial yield stress and the plastic stress, the internal structure of the soil mass sustained hardening and damage at the same time.When the loading stress was below the long-term strength, the soil structure eventually stabilized, at which time the hardening effect dominated.When the loading stress was greater than the longterm strength, the soil structure failed eventually, at which time the damage effect dominated.On this basis, viscoplastic elements considering the load hardening effect and the damage effect were established.

Figure 6 .Figure 7 .
Figure 6.Calculation results of parameters k and n.

Figure 8 .
Figure 8. Elastic modulus fitting curve of the clay in the linear elastic stage.

Figure 9 .
Figure 9. Element model of different creep mechanical property intervals on the isochronous stress ratiologarithmic strain curve.

Figure 10 .
Figure 10.Clay creep model considering the effects of load hardening and damage.

Figure 11 .
Figure 11.Comparison between the test results and the calculated results with the proposed creep model in the various creep stages of Chengdu clay.

3 .
Comprehensively considering the Chengdu clay creep test results and the established element model, the linear elastic element described the transient deformation upon loading, the viscoelastic element and the viscoplastic element considering the load hardening effect describe the decay creep deformation, the viscoplastic element considering the load damage effect describes the steady-state creep deformation, and the plastic element describes the accelerated creep deformation.Based on this, the clay creep model was established, and its rationality was verified by the creep test data.

Figure 12 .
Figure 12.Influence curves of isochronous stress ratio-logarithmic strain curve parameters k and n on creep mechanical feature points C and F.

Table 1 .
Physical and mechanical properties of Chengdu clay.

Table 2 .
Calculation parameters of Chengdu clay in the viscoelastic deformation stage.q is in kPa in the table; E 1 is in MPa; η 1 , η 2 , and η 3 are in MPa•min.