An Elastic-Viscoplastic Model for Time-Dependent Behavior of Soft Clay and Its Application on Rheological Consolidation

To describe the time-dependent behavior of soft clay, this paper extended one-dimensional Nishihara model to three-dimensional stress state based on the framework of Perzyna’s overstress theory and modified cam-clay model. The yield criterion of modified cam-clay model was used to describe the plastic properties of soft clay, and the overstress theory was used to describe the strain rate effect. Triaxial rheological tests were carried out on Ningbo soft clay and the rheological characteristics were studied. Based on laboratory results, the parameters of proposedmodel were determined by curve fitting, which show that thismodel is suitable for the rheological characteristics of Ningbo soft clay. The analysis of parameters shows that, the value of parameters changes slightly with different deviatoric stress when the confining pressure was constant, but changes notably with the increase of confining pressure. A user material subroutine of the proposed constitutive mode was coded on the platform of the FEM software ABAQUS and verified by triaxial compression of soil column. A plain strain problem was computed to analyze the rheological consolidation properties of soft clay, in which the rheological effect and the finite strain effect were considered.


Introduction
According to the classic consolidation theory presented by Terzaghi, when soft clay is consolidated under external loads, the excessive pore pressure dissipates and the effective stress increases, while the deformation grows until the excessive pore pressure dissipates completely.This theory has laid the foundation of soft soil mechanics and even became the main theoretical basis of the computation of soft soil foundation.However, in the former works, Buisman first discovered the long term deformation of soils after the excessive pore pressure was dissipated in experiments, which is also reported in the works of Zeevaert [1], Leonards and Ramiah [2], and also in lots of engineering practices.This long term deformation which cannot be explained by classic theories is usually called rheology, which contains creep, stress relaxation, long-term strength, elastic aftereffect, hysteresis effect, and so on.Because of strict requirement of post-construction settlement, researchers are paying increasing attention on the rheology of soft clay.To propose a suitable constitutive model to describe the rheology of soft clay and to determine the model parameters accurately by laboratory tests, these are the key points of the topic.Since 1960s, a number of constitutive models were proposed for time-dependent behavior of soft clay.Commonly used rheological constitutive models can be divided into three categories: empirical models, component models, and microscopic models, in which the first two categories were most widely used.Empirical models such as Singh-Mitchell model [3], Mesri et al. model [4,5] are mostly one dimensional models, which were often used to calculate the rheological deformation in one dimensional problem but were limited in complex stress state.Component models are intuitive in physical meaning, which can describe numerous rheological models by different combinations of components and can be easily extended to multidimensional problems to apply in numeral computation for complex stress state, such as the work of Adachi and Oka [6], Fodil et al. [7], Yin et al. [8,9], Hinchberger and 2 Mathematical Problems in Engineering Rowe [10], Hinchberger and Qu [11].Based on the framework of Perzyna's overstress theory [12] and modified cam-clay model, this paper proposed a three-dimensional rheological model.The model was verified by laboratory triaxial rheological tests of Ningbo soft clay, and the model parameters were determined by curve fitting.On the platform of ABAQUS, a user material subroutine of the model was coded and verified.

Model for Soft Clay
As one of commonly used component models, Nishihara model was connected by Hooke body, Kelvin body, and Bingham body, which can describe elastic deformation, creep, elastic aftereffect, and viscous flow comprehensively, such as Figure 1.When a constant external pressure  was given, the stress and strain of the model can be written as follows: In the above equation,   ,  V ,  V denote the stress of Hooke body, Kelvin body, and Bingham body, respectively; ,   ,  V ,  V denote the total strain, the strain of Hooke body, Kelvin body, and Bingham body, respectively.
Substituting the stress-strain relationship of Kelvin body and Bingham body into (1a) and (1b), the constitutive equation of one-dimensional Nishihara model can be obtained as follows: where  0 denotes the elastic modulus of Hooke body,  1 and  1 denote the elastic modulus and viscosity coefficient of Kelvin body,  2 and   denote the viscosity coefficient and yield stress of Bingham body, ⟨ −   ⟩ denotes a switch function, which can be defined as follows: Equation ( 2) contains the first two stages of creep, which are often called primary creep stage and steady creep stage.In fact, for complex stress state, (2) cannot be used to calculate the deformation of soil, it is necessary to extend (2) to three dimensional, so some hypotheses can be presented as follows: (1) the soil material is isotropic; (2) the volume deformation is just caused by spherical stress, and it is unrelated to shearing stress; (3) the Poisson's ratio of soil is constant and does not change with the stress or time.
Based on the above hypotheses, the stress-strain relationship of Hooke body is given as where  1 is the first stress invariant and  and  0 denote the bulk modulus and shear modulus of Hooke body.For Kelvin body, the stress-strain relationship is given as where  1 and  1 denote the shear modulus and 3D viscosity coefficient of Kelvin body, in which  1 = 2(1 + )/ 1 ,  is the Poisson's ratio of soil.For Bingham body, the stress-strain relationship is given as In (6),  2 denotes 3D viscosity coefficient of Bingham body, and ⟨()⟩ is a switch function to judge whether plastic yield is occurring and whether the amplitude of the plastic yield occurred,  is the plastic potential function, when the associated flow rule is adopted,  = , where  is the yield function.() can be defined as follows: where  is the current yield function,  0 is the initial yield function, and  is a constant derived from experiments; according to the work of Wei [13], the value of  can be assigned to 1.0 approximately for soft clay.
According to the work of Adachi and Oka [6] and Yin et al. [8,9], we define a static yield criterion  0 , which represents a reference yield surface for the material, such as Figure 2. Its initial shape depends on the consolidation pressure    .The expansion of the static yield surface, which describes the hardening of the material, is expressed by the variation of the consolidation pressure due to the inelastic volumetric strain  V V as follows: In (8),  0 is the intercept of initial static yield surface in the   axis, for soft clay, it can be assigned to the biggest consolidation stress in history;  is the slope of normal consolidation curve in  ∼ ln   plane;  is the slope of recompression curve in  ∼ ln   plane.
A dynamic yield criterion   is defined to describe the current state of stress as follows: From ( 8) and ( 9), the  and  0 in (7) can be written as So, we can write (7) as follows: Combining ( 6)∼( 11), ( 6) can be written as follows: From ( 4), (5), and (12), when the stress is constant, the total strain of the presented model can be written as follows:  The soil investigated in the experiments is Ningbo soft clay, which is a kind of problematic soil for low strength, high compressibility, and time-dependent behavior and deposits in Hangzhou Bay, China.The basic properties of the soil were summarized in Table 1.Both moisture content and moist unit weight of this material are significantly less than those of typical natural sedimentary deposits.

Experiments and Computation
The experimental equipment was refitted on the platform of strain controlled triaxial apparatus, the original strain loading method was changed to stress loading method, while the confining pressure system, the back pressure loading system, and the measurement system were reserved.The soil was cut into replicate specimens with a diameter of 39.8 mm and a height of 80.0 mm and placed in the triaxial test apparatus.Both soil specimens were consolidated 24 hours under a confining pressure of 100 kPa and 200 kPa, respectively; then, keep the confining pressure as constant and apply the deviatoric stress increment until the specimens were damaged or the total strain exceeds 15%.In the whole process, the free drainage conditions were kept.Table 2 showed the different loading rates applied in the experiments.In order to make sure that the creep deformation can be fully developed, each load increment lasted no less than 7 days until the deformation of the specimen is lower than 0.01 mm within 24 h.

Experimental Results.
The experimental results were presented to observe the rheological behavior of the soft clay.It is considered that deformation occurs during the whole process, but only the deformation that occurs after the dissipation of the excessive pore water pressure can be regarded as creep only.The full shearing process stress-straintime curve of each specimen is illustrated in Figure 3, and the mutistage strain-time curves of each specimen are illustrated in Figure 4.The figures show that the creep process under triaxial loading exhibits attenuation characteristics when the deviatoric stress is small, which contains two stages of creep: primary stage and steady stage; while when the deviatoric stress is large, it may exhibit accelerated characteristics, which contains three stages of creep: primary stage, steady stage and accelerated stage, such as the last load increment of specimen 2.
3.4.Discussion.Table 3 shows that the model parameters of specimen 1 ( 3 = 100 kPa) change slightly with the change of deviatoric stress except the first stage; Table 4 shows that the model parameters of specimen 2 ( 3 = 200 kPa) change slightly with the change of deviatoric stress except  0 .Compared with the model parameters of specimen 1 and specimen 2, it can be found that both  0 and  1 increase notably with the increase of  3 , which means that the elastic strain reduces when  3 increases, and  1 increases notably while  2 is reduced notably with the increase of  3 , which means that the viscoelastic strain rate and viscoplastic strain rate change significantly when  3 increases.For the lack of more experiment data, the further variation of model parameters was unable to find out, which needed further study.
Based on the parameters in Tables 3 and 4, the components of elastic strain, viscoelastic strain, and viscoplastic strain at any time can be calculated with (13).Tables 5 and 6 show    11   Viscous strain (%) Time (h) VES, q = 60 kPa VPS, q = 60 kPa VES, q = 120 kPa VPS, q = 120 kPa VES, q = 180 kPa VPS, q = 180 kPa VES, q = 80 kPa VPS, q = 80 kPa VES, q = 160 kPa VPS, q = 160 kPa VES, q = 240 kPa VPS, q = 240 kPa 1 10 be found that the elastic strain accounted for the major part of the total strain, but the viscoelastic strain and viscoplastic strain increase with the increase of  3 , which takes about 15%∼30% for each loading increment of the total and cannot be neglected.Figure 7 shows several typical curves of the development of computed viscoelastic strain and viscoplastic strain of specimen 1 and specimen 2. It can be found that the viscoelastic strain increases quickly in the early stage and then reaches steady in a short time (about 30 hours), while the viscoplastic strain increases the whole process, which indicates that the long-term deformation after consolidation of soft clay foundation is mainly caused by the viscoplasticity rather than viscoelasticity

Derivation and Verification of UMAT
4.1.Derivation of UMAT.ABAQUS is a FEM software which has powerful computing capabilities of strongly nonlinear problems.A user material subroutine (UMAT) is coded for the proposed model on the platform of ABAQUS to study the elastic-viscoplastic consolidation properties of soft clay.According to ABAQUS, the stiffness coefficient matrix (Jacobian matrix) and the equation to update the stress and strain should be defined in the UAMT.The Jacobian matrix D is defined as follows: where Δ and Δ denote the stress increment matrix and strain increment matrix, respectively.For the proposed model, ( 4)∼( 6) can be written as follows in 2D or 3D problems: In (16a), (16b), and (16c)   ,   ,  V  , and  V  denote the Voigt matrix form of stress and strain component tensors; A denotes the flexibility matrix, for plane strain problem and 3D problems, and it can be written as follows: Suppose that the time step of the th increment is Δ  , so the viscoelastic strain increment in the increment step can be written as follows: where 18) can be written as follows: The viscoplastic strain increment in the increment step can be written as follows: In order to calculate the viscoplastic strain rate at the end of the th increment, apply the Taylor series expansion and ignore the second order trace; then we can obtain the following equation: Suppose that H  =  ε V  /; (21) can be written as Substituting ( 22) into (21), we can obtain the following equation: From (16c), H  can be written as where a = /.According to (9), a can be written as follows: where From (25), the following equation can be derived: where , and M 3 are matrices written as follows: Combining the above equations, the value of H  can be derived.Substituting H  into (22), the viscoplastic strain increment Δ V  is obtained.For current increment, Δ  is calculated by (28) as follows: Substituting ( 19) and ( 23) into (28), we can obtain the following equation: Equation ( 29) is the relationship between stress increment and strain increment, in which matrix D is defined as follows: Based on the above equations, a user material subroutine was coded for the proposed constitutive model on the platform of ABAQUS.

Verification of UMAT.
In order to verify the validity of the UMAT, we study a soil column of triaxial compression such as Figure 8, the size of the column is 0.1 m × 0.1 m × 0.2 m.Supposing that  1 = 100 kPa,  2 =  3 = 50kPa, and the model parameters are given in Table 7.The column was cut into one element, the element type is C3D8, and the consolidation effect was not considered here.Embedding the UMAT coded in this paper on the platform of ABAQUS, we calculated the triaxial compression property of the column; the results were summarized in Figure 9∼Figure 10.
From Figure 9(a), we can see the deformation of the column exhibits apparent rheology, the rheological deformation takes about 37% of the whole deformation.Figure 9(b) shows that the static yield surface gradually expanded with the time until it reached steady state, when it coincided with the dynamic yield surface.From Figures 9(c) and 9(d), we can see that the viscoelastic strain and the viscoplastic strain both gradually increased with time until they reached their steady state; accordingly, the viscoelastic strain rate and the viscoplastic strain rate were gradually deduced to nearly zero.The above results fit the property of Nishihara model and the hardening law of modified Cambridge model very well, which proved the correctness of the UMAT.
In the proposed model,  1 ,  2 , and  0 are the most important parameters; Figure 10 showed the sensitivity of these parameters.From Figures 10(a) and 10(b), we can see the final viscous deformation is a constant regardless of the different values of  1 ,  2 , but the time when it reaches steady varies greatly, so, the accuracy of  1 and  2 is directly related to the accuracy of computed displacement in the rheological process.Figure 10(c) shows that the value of  0 has a great influence on the viscoplastic strain: the smaller  0 is, the more rapidly the viscoplastic strain develops and the larger the final viscoplastic strain occurs, and so, the accuracy of  0 is not only directly related to the accuracy of computed displacement in the rheological process, but also affected the final displacement, which needs to be paid more attention.

Analysis on Rheological Consolidation of Soft Clay
5.1.Computed Model.Consolidation problems can be classified as small strain theory and finite strain theory according to the difference of geometric assumption.Small strain theory assumed that the strain is little, but the consolidation deformation is usually very large for deep soft clay foundation.Take dredger fill for example, the deformation may exceed 80%, for which the small strain theory is no longer suitable [14,15].As the pioneer researchers studied on finite strain consolidation theory, Mikasa [16] and Gibson et al. [17] derived their one-dimensional finite strain consolidation equations, respectively, in 1960s, and many scholars developed further research on the topic in the following decades.Because of strong nonlinearity, most of existing researches on finite strain consolidation adopted some simplified assumptions without considering rheological characteristics.By means of ABAQUS, this paper studied the rheological consolidation of soft clay considering finite strain effect.The computed model is plane strain problem as shown in Figure 11.The soil foundation is homogeneous, anisotropic, saturated, and normally consolidated soft clay with a thickness of 10.0 m, and only the top boundary is permeable.In addition, the semilogarithmic empirical formula proposed by Tavenas et al. [18] was adopted to consider the variation of permeability in the process of consolidation as follows: where  V0 is the initial permeability coefficients, and  V is the permeability coefficient in the consolidation process;  0 is the initial void ratio, and  is the void ratio in the consolidation process;   is the permeability index.According to the above work, the model parameters are summarized in Table 8 approximately.

Consolidation Property of Soft Clay.
From Figure 12 we can see the differences of the consolidation process between finite strain and small strain.Figure 12(a) shows that for any location the excess pore pressure of finite strain is less than that of small strain at the same time, and the difference is increasing with the increase of time in beginning period.Figure 12(b) shows that the displacement of finite strain is less than that of small strain and the difference increased significantly if the load is increasing.For example, when  = 200 kPa, the difference may take about 20% while it is nearly 0 when  = 50 kPa.So, it is essential to consider the error caused by small strain if the load is large and the compressibility of the soil is large.
From Figure 13 we can see the influence of rheology on the consolidation process.Figure 13(a) shows that the dissipation of excess pore pressure is significantly slower if the rheological effect is considered; for example, when  = 1000 d, the excess pore pressure is nearly 0 if the rheological effect is not considered, but it is still larger than 20 kPa if the rheological effect is considered.Figure 13(a) shows that the displacement is influenced by rheological effect too; for example, in the early stage ( ≤ 100 d), the displacement is nearly the same for two cases, but the final displacement may have a difference of more than 35%, which indicates that the displacement in the early stage is mainly caused by consolidation, while the displacement in the late stage is mainly caused by rheology.SS, q = 100 kPa FS, q = 50 kPa SS, q = 50 kPa SS, q = 200 kPa FS, q = 200 kPa FS, q = 100 kPa (4) In the proposed model, the viscoelastic viscosity coefficient  1 , the viscoplastic viscosity coefficient  2 , and the initial yield surface  0 are the most important parameters. 1 and  2 have important influence on the time of viscous strain developed, and  0 have important influence on both the viscoplastic strain rate and the final viscoplastic strain.

Conclusions
(5) The existence of rheology will significantly slow down the dissipation of the excess pore pressure.Considering the rheological effect, the total displacement is larger, and the displacement in the early stage is mainly caused by consolidation, while in the late stage it is mainly caused by rheology.
(6) For consolidation problems of high compressibility soft clay foundation, great error will be caused by small strain assumption, so finite strain effect should be considered.When finite strain effect is considered, the consolidation process is developed more rapidly, and the computed consolidation displacement will be less than that of small strain.

Figure 9 :
Figure 9: Triaxial compression property of soil column.

( 1 )
An elastic-viscoplastic model developed to describe the time-dependent behavior of normally consolidated soft clay was presented in the paper on the

Table 1 :
Property of testing material.

Table 2 :
Loading scheme of triaxial rheological tests.
, we can see that the model fits testing data of primary creep stage and steady creep stage very well, while for the accelerated creep stage, it needs further study.According to the value of

Table 7 :
Material parameters of verification model.