# An EVP Double-Yield Nonstationary Flow Surface Model for Soft Clay

Authors: Liucheng Chang; Hongyu Wang; Ya Wang
Journal: Advances in Civil Engineering
 (2022)

---

## Abstract

According to experimental findings, soft clay has compression and shear creep. The goal of the study is to develop an elastic-viscoplastic (EVP) model for the viscous behavior of soft clay. The new EVP model is deduced from the associated flow rule, Yin’s double-yield surface theory, and the critical state soil mechanics. The nonstationary flow surface model with double-yield surface, compared to other existing models, considers the compression and shear creep properties of soft clay. The model parameters can be identified from two sorts of tests (i.e., the general triaxial tests and triaxial rheological tests). Different types of soft clay, such as Tianjin Coastal clay, Shantou dredger fill, and silted soil in front of dam, were tested to verify the accuracy of the constructed model. All comparisons of computed and experimental findings revealed that the new EVP model accurately reflected the time-dependent behaviors of soft clay in both shear and compression creep.


---

## Body

## 1. Introduction

As many engineering cases showed, soft clay obviously shows time-dependent behaviors [1]. For lack of concern for the rheological behaviors of soft clay during the design and construction period [2], many safety problems have occurred in large infrastructure projects constructed in soft ground conditions in the past decades. For example, more than 3000 people were killed in a massive landslide in the south slope of Vaiont reservoir in 1963, which was due primarily to the creep deformation of the weak layer [3]. Moreover, Kanto Airport’s postconstruction settlement, which opened in 1994, grew at a rate of 7 cm per year, having attained a 50-year design value.Soft clay’s engineering properties (e.g., low-intensity strength and large postconstruction settlement) still make projects built in soft clay distribution areas increasingly challenging. Therefore, it is extremely urgent to develop an available constitutive model for the viscous materials so as to model the time-dependent behaviors of soft clay.Many models have been proposed in geotechnical engineering to capture the time-dependent properties of soft clay since the creep characteristics of soft clay became a hot topic. Chaboche [4] and Karim and Gnanendran [5] put forward a concise overview of constitutive theories and introduced simply the features of EVP models and their advantages and shortcomings. Besides, the relevant reviews can also be found in the literature of Owen and Hinton [6] and Desai and Siriwardane [7]. The EVP models can be primarily divided into four categories [1]: (1) empirical models, such as primary empirical relations [2, 8] and semiempirical secondary relations [9, 10]; (2) rheological models, such as differential approach [11] and engineering theories of creep [12]; (3) general stress-strain-time models, such as overstress models [13, 14] and nonstationary flow surface (NSFS) models [15]; (4) other contributions, such as Borja’s model [16] and subloading surface model [17]. The first two categories can help further the phenomenological understanding of the creep properties of soft clay. Therefore, the first two models are commonly regarded as the foundation of three-dimensional constitutive relations. However, empirical and rheological models are limited to only simulating the strain-stress-time relation of soft clay with specific boundaries and loading conditions [1]. On the contrary, the general stress-strain-time models can consider time effects and in-viscid deformation under any possible loading conditions, which prompts them to be widely used, especially in finite element analysis. However, most stress-strain-time models have only one yield surface, resulting in the insufficient accuracy of existing available creep models in describing the creep behaviors of soft clay. Furthermore, hydrostatic pressure and shear stress actually coexist between soil particles and bound water; the time-dependent response of soft clay can be demonstrated by two different characteristics: shear creep and compression creep, which can be better described using double-yield surfaces [18–20].In this paper, a practical model for the time-dependent behaviors of soft clay is developed. The model’s main characteristics are as follows:(1)
The developed model is a double-yield surface hardening model based on the nonstationary flow surface theory and Yin’s theoretical framework. General triaxial and triaxial rheological tests can be used to directly determine model parameters.(2)
The developed model, compared to the prior EVP model, takes into account both compression creep and shear creep appropriateness for time-dependent properties in various confining pressures.In the following, the double-yield nonstationary flow surface EVP model is presented using the theory of the nonstationary flow surface and Yin’s double-yield surface theory. The presented model is then put into the finite element program Abaqus and validated with a variety of soft clays.


## 2. The Fundamental Assumption in the EVP Model

Reviewing the literature, we find that Skrzypek [21] divided total strains into instantaneous elastic segments and time-dependent components, which have been widely accepted by scholars in geotechnical engineering for decades. For time-dependent deformation, there are currently two main hypotheses: Borja and Kavazanjian [16] and [22] proposed that total strain rates are subdivided into elastic strain rates, plastic strain rates, and viscoplastic strain rates. In comparison, Desai and Zhang [23], Yin and Graham [24], and Kabwe et al. [25] argue that total strain rates can be divided into elastic strain rates and viscoplastic strain rates. Afterward, [10] and Kong et al. [18] believed that viscoplastic strain rates could be further subdivided into volumetric strain rates and deviatoric strain rates. In this paper, the latter hypothesis is accepted; i.e., total strain rates are the sum of elastic components, volumetric components, and deviatoric components.(1)ε˙=ε˙ije+ε˙ijvp,(2)ε˙ijvp=ε˙vijvp+ε˙sijvp,where ε˙ is the elastic strain rates, ε˙ijvp is the viscoplastic strain rates, ε˙vijvp is the volumetric creep strain rates, and ε˙sijvp is the shear creep strain rates. The superscripts e and vp and the subscripts v and s denote elastic strain, viscoplastic strain, volumetric creep, and shear creep, respectively; the overdot denotes the rate of a variable. The elastic strain rates ε˙ije in (1) are supposed to obey generalized Hooke’s law and can be calculated:(3)ε˙ije=Cijklσ˙kl′,where Cijkj is the fourth-order flexibility tensor and σ˙kl′ is an effective stress rate. Supposing soft clay is an isotropic material, there are only two elastic constants in the flexibility matrix: the elastic modulus Ee and the Poisson ratios ve.The viscoplastic strain ratesε˙ijvp in (1) are defined according to the nonstationary flow surface theory presented by Olszak and Perzyna [13]. Flow surface can be represented as follows:(4)F=Fσij′,Ks,β,where F is the loading function of σij′, Ks, and β; σij′ is the current effective stress; KS is a strain-hardening parameter related to the viscoplastic strains εijvp; and β is a time variable associated with clay proper. ε˙ijvp can be determined based on the assumption of associative normality flow rule, where the loading function F is equal to the potential function g.(5)ε˙ijvp=∂g∂σij′=∂F∂σij′,where g is the potential function.The viscoplastic deformation can be obtained by differentiating (4):(6)F˙=∂F∂σij′σ˙ij′+∂F∂ε˙klvp∂ε˙klvp+∂F∂ββ˙.The constitutive equations are based on the nonstationary flow surface theory by the use of (1), (3), and (5), and the generalized Hooke’s law can be written as follows:(7)ε˙=Cijklσ˙kj′+Ctpdσij′dt,where Ctp=−∂F/∂σij′/∂F/∂εklvp−∂F/∂β∂β/∂t/∂F/∂εklvp is the viscoplastic flexibility matrix.


## 3. Model Description

### 3.1. Elastic Constitutive Relation

Supposing soft clay is an isotropic material, the stress-strain relationship submits general Hooke’s law, and the strain rates are dependent on effective stress rates; then, the viscoplastic flexibility matrixCijkj in (3) can be written as follows:(8)Cijkl=1E,i=j,i,j=1,2,3,vE,i≠j,i,j=1,2,3,G,i=j,i,j=4,5,6,where G is the shear modulus and can be calculated according to (9), E is the elastic modulus and can be determined utilizing (10), and v is the Poisson’s ratio. For all clay, Poisson’s ratio usually equals 0.3, due primarily to the smaller proportion of elastic strain in total strain [26]).(9)G=KGpappan,(10)E=2G1+v,where KG and n are material parameters, and pa is the atmospheric pressure and is introduced to make parameters dimensionless. The relationship of both G/pa and p/pa is a straight line in dual logarithmic coordinate system, and the intercept and slope of the line are KG and n, respectively.


### 3.2. Compression Creep Constitutive Relation

By combining associated flow rule and critical state soil mechanics theory, Roscoe et al. [27] proposed the Cam-clay model with a bullet yield surface by triaxial drained or undrained tests on normally consolidated and feeble-overconsolidated clay. The Cam-clay model can exceedingly calculate elastic-plastic deformation, particularly plastic volume deformation. However, it cannot reflect shear deformation. As a result, Roscoe and Burland [28] mended the Cam-clay model and proposed the modified Cam-clay model. The modified Cam-clay model is expressed as follows:(11)pc′εvp=p′+q′2M2p′,where p′=σ8′ is the average effective stress and q′=3/2τ8′ is the generalized shear stress. More detailed expressions for p′ and q′ can be found in the elastoplastic mechanics. In the principal stress space, σ8′ and τ8′ are referred to as octahedral normal and shear stresses, respectively; M is the slope of the critical state line (CSL); and pc′ is the hardening function taking the plastic volume strain εvp as hardening parameter. The yield curve of the modified Cam-clay model is an ellipse in the p′−q′ plane, as indicated in Figure 1. The shape of the yield curve (i.e., the ratio of the long axis and short axis of the ellipse) depends only on M but does not correlate with the deformation behaviors of clay, which is not consistent with actual situations for much soft clay.Figure 1
The yield curves of modified Cam-clay model.In general, the effective cohesion of normal consolidated clay under drained shear conditions is equal to zero. In practice, the clay layer is pre-consolidated to a certain extent, so the cohesion force is not equal to zero. For this purpose [29], who replaced p′ with p′+pr′ extended modified Cam-clay model to the condition that cohesion force exists. Equation (11) can be rewritten as follows:(12)p0′εvp=p′+q′2M2p′+pM′,where p0′εvp=pc′εvp−pr′ is the hardening function; like pc′, pr′ is the intercept of the CSL on the transverse axis in the p′−q′ plane. Equation (12) is the modified Cam-clay model proposed by [29]. The yield curve of the modified model is depicted in Figure 2. The modified Cam-clay model can simulate the stress-strain properties of clay exceptionally under different conditions. Consequently, the modified Cam-clay model promoted by [29] is essentially used in the paper.Figure 2
The yield curve of the modified Cam-clay model developed by [29].For the modified Cam-clay model, thee−lnp relationship is available in the construction of pc′−εvp relation. However, the relationship is not linear when p′ value is small, especially for compacted soils, as shown in Figure 3. The linear relation between e and lnp will produce certain difficulties in mathematical treatment. Therefore, it is a necessity to modify the approach for the modified Cam-clay model [29] formulated the relationship between pc′ and εv by a hyperbolic function on the basis of experiment analyses. Figure 4(a) shows the pc′/pa−εv relation curve in triaxial isobaric conditions, where the paper takes the atmospheric pressure pa for model parameters dimensionless.Figure 3
Changing curves ofe and lnp′.Figure 4
The experimental data for thepc′−εv relationship [30]: (a) pc′/pa vs. εv; (b) pc′/paεv vs. pc′/pa.
(a)(b)The relation betweenpc′/pa and εv shown in Figure 4(a) is nearly linear in pc′/paεv−pc′/pa coordinates, as depicted in Figure 4(b). The pc′/pa−εv curve can be expressed as follows:(13)pc′paεv=l+kpc′pa=Ei+pc′paεv,where k is the slope of the straight line, l=1/εva is the intercept of the straight line on the horizontal axis, εva is the asymptotic value of εv, andEi is the initial tangent slope of the hyperbola.The relationship betweenεvp and pc′ can be described with a hyperbolic function. Compared with plastic volumetric strain, elastic volume strain is very small; hence, it is possible to assume that plastic volumetric strain is approximately equal to the total volumetric strain. Then, (13) can be also written as follows:(14)εv≈εvp=pc′kpc′+lpa.The strain-hardening functionpc′ can be expressed as follows:(15)pc′=lεvp1−kεvppa.By substituting (15) into (13) and replacing M in (13) with M1, we can obtain(16)p′+q′2M12p′+pr′−lεvp1−kεvppa=0,where M1=1.2M is a parameter that reflects the variety of stress-strain curves. Equation (16) is Yin’s double-yield surface model [26].Triaxial compression creep tests are done to obtain the changeable rule of the volumetric strainsεvt over time, as indicated in Figure 5. It is important to note that the time-dependent volumetric strain εvt does not include the instantaneous volume strain εvp. Many scholars describe the εvt−t relationship with a hyperbola. The hyperbolic function can be given as follows:(17)εvt=tcv+dvt,where cv=1/h, h is the initial tangent slope of the hyperbola, dv=1/εvu, and εvu is the asymptotic value of εvt. Supposing that the ratio of ενu and εvp is constant under different stress conditions, i.e., εvu/εvp=ξ1, εvt can be further expressed by h and ξ1 as follows:(18)εvt=t1/h+t/ξ1εvp.Figure 5
Relation curveεvt and t under triaxial compression creep conditions.After rearranging (17), the expression for εvp can be obtained:(19)εvp=ξ1th−εvthtεvt.By substituting (19) into (16), the following equation can be obtained:(20)p′+q′2M12p′+pr′+lξ1th−εvthktξ1−htεvt−kξ1εvtpa=0.Then, the compression creep surface function can be given as follows:(21)FVp′,q′,εvt,t=p′+q′2M12p′,+pr′+HV+βV,(22)HV=−lξ1εvthktξ1−htξ1−kξ1εvtpa,βV=lξ1thdv+lξ1dvhktξ1−kξ1−htpa,where HV is the strain-hardening function related to the viscoplastic volumetric strains εvt, and βV is a time variable associated with soil property. In order to simplify the model structure and reduce the number of model parameters, the constitutive model is formulated by applying the associated flow theory, where the loading function is identical to the plastic potential function.


### 3.3. Shear Creep Constitutive Relation

Volume expansion causes two dilemmas: (1) the decrease of mean normal stresses; (2) the production of shear stresses. Plastic expansion is only dependent on shear strains without considering tensile stresses. It is reasonable to hypothesize that the plastic deviatoric strainεsp is selected as the hardening parameter of the shear creep yield criterion.(23)εsp=∫dεsp=∫23A,A=dε1p−dε2p2+dε2p−dε3p2+dε3p−dε1p2.Hsieh et al. (1990) recommended a hardening law for shear creep based on triaxial tests for clay as follows:(24)Fs=aq′Gq′M2p′p′+pr′−q′−εsp,where G is the elastic shear modulus, M2 is the parameter that is lower than M, anda is a material parameter as indicators of the magnitude for shear dilation or shear contraction. Equation (24) is the parabolic yield function in Yin’s double-yield surface model [26], as indicated in Figure 6.Figure 6
Conceptual shear yield trajectory in thep′−q′ plane.Triaxial shear creep experiments are done to obtain the variation of the shear strainεSt over time, as indicated in Figure 7. A hyperbolic function for εSt and t can be proposed as follows:(25)εSt=tcs+dst,where cs=1/ς is a material parameter, ς is the initial tangent slope of the hyperbola, ds=1/εsu is a material parameter, and εsu is the asymptotic value of εSt. It is reasonable to suppose that the ratio of εsu and εsp is constant under different stress conditions; i.e., εsu/εsp=ξ2. εst can be further expressed by ς and ξ2 as follows:(26)εst=t1/ς+t/ξ2εsp.Figure 7
Relation curveεst and t under triaxial conditions.After rearranging (26), the expression for εSp can be obtained:(27)εSp=ξ2tς−εstςεstt.By substituting (27) into (24), the shear creep surface function can be given as follows:(28)FSp′p′,q′,εst,t=aq′Gq′M2p′+q′−q′+HS+βS,(29)HS=ξ2εst,βS=−ξ2ςt,where HS is the strain-hardening function related to the viscoplastic shear strain εst, and βS is the time variable associated with soil property. The constitutive model is also formulated by applying the associated flow theory.The compression creep function FV and the shear creep function FS are illustrated in Figure8. It can be seen that the two yield curves separate the p′−q′ plane into four regions: A1 indicates elastic zone; A4 indicates the viscoplastic strains associated with the compression yield; A2 indicates the viscoplastic deformation related to the shear yield; A3 indicates the viscoplastic flow linked to the two yield criteria. Thus, the EVP model, which is a comprehensive model used in soft clay, can simulate the dilatancy and time-dependent behaviors of soft clay.Figure 8
Volume yield surface and shear creep shear yield surface in thep′−q′ plane.


### 3.4. Viscoplasticity Flexibility Matrix

It can be concluded from (7) that the viscoplastic flexibility matrix Ctp is key to determining the viscoplastic strain rate ε˙ijvp; how to acquire Ctp will be detailed in this section.Based on the above analysis,Ctp can be further written as follows:(30)Ctp=−∂FV/∂p′∂p′/∂σij′A−∂FV/∂p′∂p′/∂σij′A−∂FV/∂βV∂βV/∂tA−∂FS/∂p′∂p′/∂σij′B−∂FS/∂p′∂p′/∂σij′B−∂FS/∂βS∂βS/∂tB,where A=∂FV/∂εvt, B=∂FS/∂εst. The different terms in (30) can be expressed as follows:(31)∂FV∂q′=2q′2M12p′+pr′2,∂FV∂p′=1−q′2M12p′+pr′2,∂FS∂q′=−aq′2M22GM2p′+M2pr′−q′2q′/M2p′+pr′−q′,A=∂FV∂εvt=−ξ1A2A3A1,A1=hpaltkξ12Q+βVthk−th−kξ1,A2=ξ1htQ+βV+kQ+βV+lpa,A3=t2Q+βV2k−1h2+ktξ1Q+βV2−2h−Q+βV2k2+l2pa2ξ12,Q=p′+q′2M12p′+pr′,B=∂FS∂εSt=−S+βS2ξ2,S=aq′Gq′M2p′+pr′−q′,∂FV∂βV∂βV∂t=−lξ1dvhkξ1+dvk2ξ1−htdvhkξ1−ht−kξ12,∂FS∂βS∂βS∂t=−ξ2ςt2,∂p′∂σij′=δij,whereδij=13,ifi=j,0,ifi≠j,∂q′∂σij′=32q′σij′−p′,ifi=j,3q′σij′,ifi≠j,


### 3.5. Parameter Determination

There are fifteen parameters in the proposed model, includingν, Poisson’s ratio; E, elastic modulus; KG and n, elastic material parameters; M1 and M2, modified strength parameters; pr′, the intercept of the CSL on the transverse axis in the p′−qσ plane; a, the parameter as indicator of magnitude for shear dilation or shear contraction; l and k, the parameters related to triaxial compression tests; dv, h, and ξ1, the parameters related to compression creep; and ξ2 and ς, the parameters associated with shear creep.Among the fifteen parameters,a, M1, l, k, and pr′ in (16) can be obtained from the triaxial compression tests or one-dimensional oedometer tests. ν, E, KG, n, M2, and pr′ in (24) can be determined from the triaxial shear tests. dv, h, and ξ1 can be determined using oedometer creep tests or triaxial compression creep tests. ξ2 and ς can be determined utilizing triaxial shear creep tests.To determine the parametersν, E, KG, n, M2, and pr′ in (24), many triaxial shear tests are required to gain a collection of q′−εs curves, as shown in Figure 9. Gs=3Gi is the initial tangent of q′−εs curves, and Gi is the initial tangent modulus. G can be calculated according to (32). KG can be obtained by drawing the relationship between G/pa and p′/pa in dual logarithmic coordinate system. v=0.3, and E can be determined using (10). M and pr′ can be gained by drawing the CSL in p′−q′ plane. M2 can be evaluated utilizing (33).(32)G=23Gs,(33)M2=MRf0.25,where Rf=qf′/qult′ is the failure ratio, qf′ is the failure stress and is gained from triaxial shear tests, and qult′ is the asymptotic value of q.Figure 9
Schematic stress-strain curves for the triaxial shear test.The parametersl, k, and pr′ in (16) should be obtained from the triaxial compression or one-dimensional oedometer tests. Section 3.2 explains how they arrived at their values. a=0.25−0.15d, d is the slope of the last part of εv−εa curve, as illustrated in Figure 10.Figure 10
Theεa−εv curve.The parametersdv, h, and ξ1 can be determined using oedometer creep tests or triaxial compression creep tests, and their determination has been presented in Section 3.2. The parameters ξ2 and ς can be determined utilizing triaxial shear creep tests, and Section 3.3 gives their detailed calculation process.


## 3.1. Elastic Constitutive Relation

Supposing soft clay is an isotropic material, the stress-strain relationship submits general Hooke’s law, and the strain rates are dependent on effective stress rates; then, the viscoplastic flexibility matrixCijkj in (3) can be written as follows:(8)Cijkl=1E,i=j,i,j=1,2,3,vE,i≠j,i,j=1,2,3,G,i=j,i,j=4,5,6,where G is the shear modulus and can be calculated according to (9), E is the elastic modulus and can be determined utilizing (10), and v is the Poisson’s ratio. For all clay, Poisson’s ratio usually equals 0.3, due primarily to the smaller proportion of elastic strain in total strain [26]).(9)G=KGpappan,(10)E=2G1+v,where KG and n are material parameters, and pa is the atmospheric pressure and is introduced to make parameters dimensionless. The relationship of both G/pa and p/pa is a straight line in dual logarithmic coordinate system, and the intercept and slope of the line are KG and n, respectively.


## 3.2. Compression Creep Constitutive Relation

By combining associated flow rule and critical state soil mechanics theory, Roscoe et al. [27] proposed the Cam-clay model with a bullet yield surface by triaxial drained or undrained tests on normally consolidated and feeble-overconsolidated clay. The Cam-clay model can exceedingly calculate elastic-plastic deformation, particularly plastic volume deformation. However, it cannot reflect shear deformation. As a result, Roscoe and Burland [28] mended the Cam-clay model and proposed the modified Cam-clay model. The modified Cam-clay model is expressed as follows:(11)pc′εvp=p′+q′2M2p′,where p′=σ8′ is the average effective stress and q′=3/2τ8′ is the generalized shear stress. More detailed expressions for p′ and q′ can be found in the elastoplastic mechanics. In the principal stress space, σ8′ and τ8′ are referred to as octahedral normal and shear stresses, respectively; M is the slope of the critical state line (CSL); and pc′ is the hardening function taking the plastic volume strain εvp as hardening parameter. The yield curve of the modified Cam-clay model is an ellipse in the p′−q′ plane, as indicated in Figure 1. The shape of the yield curve (i.e., the ratio of the long axis and short axis of the ellipse) depends only on M but does not correlate with the deformation behaviors of clay, which is not consistent with actual situations for much soft clay.Figure 1
The yield curves of modified Cam-clay model.In general, the effective cohesion of normal consolidated clay under drained shear conditions is equal to zero. In practice, the clay layer is pre-consolidated to a certain extent, so the cohesion force is not equal to zero. For this purpose [29], who replaced p′ with p′+pr′ extended modified Cam-clay model to the condition that cohesion force exists. Equation (11) can be rewritten as follows:(12)p0′εvp=p′+q′2M2p′+pM′,where p0′εvp=pc′εvp−pr′ is the hardening function; like pc′, pr′ is the intercept of the CSL on the transverse axis in the p′−q′ plane. Equation (12) is the modified Cam-clay model proposed by [29]. The yield curve of the modified model is depicted in Figure 2. The modified Cam-clay model can simulate the stress-strain properties of clay exceptionally under different conditions. Consequently, the modified Cam-clay model promoted by [29] is essentially used in the paper.Figure 2
The yield curve of the modified Cam-clay model developed by [29].For the modified Cam-clay model, thee−lnp relationship is available in the construction of pc′−εvp relation. However, the relationship is not linear when p′ value is small, especially for compacted soils, as shown in Figure 3. The linear relation between e and lnp will produce certain difficulties in mathematical treatment. Therefore, it is a necessity to modify the approach for the modified Cam-clay model [29] formulated the relationship between pc′ and εv by a hyperbolic function on the basis of experiment analyses. Figure 4(a) shows the pc′/pa−εv relation curve in triaxial isobaric conditions, where the paper takes the atmospheric pressure pa for model parameters dimensionless.Figure 3
Changing curves ofe and lnp′.Figure 4
The experimental data for thepc′−εv relationship [30]: (a) pc′/pa vs. εv; (b) pc′/paεv vs. pc′/pa.
(a)(b)The relation betweenpc′/pa and εv shown in Figure 4(a) is nearly linear in pc′/paεv−pc′/pa coordinates, as depicted in Figure 4(b). The pc′/pa−εv curve can be expressed as follows:(13)pc′paεv=l+kpc′pa=Ei+pc′paεv,where k is the slope of the straight line, l=1/εva is the intercept of the straight line on the horizontal axis, εva is the asymptotic value of εv, andEi is the initial tangent slope of the hyperbola.The relationship betweenεvp and pc′ can be described with a hyperbolic function. Compared with plastic volumetric strain, elastic volume strain is very small; hence, it is possible to assume that plastic volumetric strain is approximately equal to the total volumetric strain. Then, (13) can be also written as follows:(14)εv≈εvp=pc′kpc′+lpa.The strain-hardening functionpc′ can be expressed as follows:(15)pc′=lεvp1−kεvppa.By substituting (15) into (13) and replacing M in (13) with M1, we can obtain(16)p′+q′2M12p′+pr′−lεvp1−kεvppa=0,where M1=1.2M is a parameter that reflects the variety of stress-strain curves. Equation (16) is Yin’s double-yield surface model [26].Triaxial compression creep tests are done to obtain the changeable rule of the volumetric strainsεvt over time, as indicated in Figure 5. It is important to note that the time-dependent volumetric strain εvt does not include the instantaneous volume strain εvp. Many scholars describe the εvt−t relationship with a hyperbola. The hyperbolic function can be given as follows:(17)εvt=tcv+dvt,where cv=1/h, h is the initial tangent slope of the hyperbola, dv=1/εvu, and εvu is the asymptotic value of εvt. Supposing that the ratio of ενu and εvp is constant under different stress conditions, i.e., εvu/εvp=ξ1, εvt can be further expressed by h and ξ1 as follows:(18)εvt=t1/h+t/ξ1εvp.Figure 5
Relation curveεvt and t under triaxial compression creep conditions.After rearranging (17), the expression for εvp can be obtained:(19)εvp=ξ1th−εvthtεvt.By substituting (19) into (16), the following equation can be obtained:(20)p′+q′2M12p′+pr′+lξ1th−εvthktξ1−htεvt−kξ1εvtpa=0.Then, the compression creep surface function can be given as follows:(21)FVp′,q′,εvt,t=p′+q′2M12p′,+pr′+HV+βV,(22)HV=−lξ1εvthktξ1−htξ1−kξ1εvtpa,βV=lξ1thdv+lξ1dvhktξ1−kξ1−htpa,where HV is the strain-hardening function related to the viscoplastic volumetric strains εvt, and βV is a time variable associated with soil property. In order to simplify the model structure and reduce the number of model parameters, the constitutive model is formulated by applying the associated flow theory, where the loading function is identical to the plastic potential function.


## 3.3. Shear Creep Constitutive Relation

Volume expansion causes two dilemmas: (1) the decrease of mean normal stresses; (2) the production of shear stresses. Plastic expansion is only dependent on shear strains without considering tensile stresses. It is reasonable to hypothesize that the plastic deviatoric strainεsp is selected as the hardening parameter of the shear creep yield criterion.(23)εsp=∫dεsp=∫23A,A=dε1p−dε2p2+dε2p−dε3p2+dε3p−dε1p2.Hsieh et al. (1990) recommended a hardening law for shear creep based on triaxial tests for clay as follows:(24)Fs=aq′Gq′M2p′p′+pr′−q′−εsp,where G is the elastic shear modulus, M2 is the parameter that is lower than M, anda is a material parameter as indicators of the magnitude for shear dilation or shear contraction. Equation (24) is the parabolic yield function in Yin’s double-yield surface model [26], as indicated in Figure 6.Figure 6
Conceptual shear yield trajectory in thep′−q′ plane.Triaxial shear creep experiments are done to obtain the variation of the shear strainεSt over time, as indicated in Figure 7. A hyperbolic function for εSt and t can be proposed as follows:(25)εSt=tcs+dst,where cs=1/ς is a material parameter, ς is the initial tangent slope of the hyperbola, ds=1/εsu is a material parameter, and εsu is the asymptotic value of εSt. It is reasonable to suppose that the ratio of εsu and εsp is constant under different stress conditions; i.e., εsu/εsp=ξ2. εst can be further expressed by ς and ξ2 as follows:(26)εst=t1/ς+t/ξ2εsp.Figure 7
Relation curveεst and t under triaxial conditions.After rearranging (26), the expression for εSp can be obtained:(27)εSp=ξ2tς−εstςεstt.By substituting (27) into (24), the shear creep surface function can be given as follows:(28)FSp′p′,q′,εst,t=aq′Gq′M2p′+q′−q′+HS+βS,(29)HS=ξ2εst,βS=−ξ2ςt,where HS is the strain-hardening function related to the viscoplastic shear strain εst, and βS is the time variable associated with soil property. The constitutive model is also formulated by applying the associated flow theory.The compression creep function FV and the shear creep function FS are illustrated in Figure8. It can be seen that the two yield curves separate the p′−q′ plane into four regions: A1 indicates elastic zone; A4 indicates the viscoplastic strains associated with the compression yield; A2 indicates the viscoplastic deformation related to the shear yield; A3 indicates the viscoplastic flow linked to the two yield criteria. Thus, the EVP model, which is a comprehensive model used in soft clay, can simulate the dilatancy and time-dependent behaviors of soft clay.Figure 8
Volume yield surface and shear creep shear yield surface in thep′−q′ plane.


## 3.4. Viscoplasticity Flexibility Matrix

It can be concluded from (7) that the viscoplastic flexibility matrix Ctp is key to determining the viscoplastic strain rate ε˙ijvp; how to acquire Ctp will be detailed in this section.Based on the above analysis,Ctp can be further written as follows:(30)Ctp=−∂FV/∂p′∂p′/∂σij′A−∂FV/∂p′∂p′/∂σij′A−∂FV/∂βV∂βV/∂tA−∂FS/∂p′∂p′/∂σij′B−∂FS/∂p′∂p′/∂σij′B−∂FS/∂βS∂βS/∂tB,where A=∂FV/∂εvt, B=∂FS/∂εst. The different terms in (30) can be expressed as follows:(31)∂FV∂q′=2q′2M12p′+pr′2,∂FV∂p′=1−q′2M12p′+pr′2,∂FS∂q′=−aq′2M22GM2p′+M2pr′−q′2q′/M2p′+pr′−q′,A=∂FV∂εvt=−ξ1A2A3A1,A1=hpaltkξ12Q+βVthk−th−kξ1,A2=ξ1htQ+βV+kQ+βV+lpa,A3=t2Q+βV2k−1h2+ktξ1Q+βV2−2h−Q+βV2k2+l2pa2ξ12,Q=p′+q′2M12p′+pr′,B=∂FS∂εSt=−S+βS2ξ2,S=aq′Gq′M2p′+pr′−q′,∂FV∂βV∂βV∂t=−lξ1dvhkξ1+dvk2ξ1−htdvhkξ1−ht−kξ12,∂FS∂βS∂βS∂t=−ξ2ςt2,∂p′∂σij′=δij,whereδij=13,ifi=j,0,ifi≠j,∂q′∂σij′=32q′σij′−p′,ifi=j,3q′σij′,ifi≠j,


## 3.5. Parameter Determination

There are fifteen parameters in the proposed model, includingν, Poisson’s ratio; E, elastic modulus; KG and n, elastic material parameters; M1 and M2, modified strength parameters; pr′, the intercept of the CSL on the transverse axis in the p′−qσ plane; a, the parameter as indicator of magnitude for shear dilation or shear contraction; l and k, the parameters related to triaxial compression tests; dv, h, and ξ1, the parameters related to compression creep; and ξ2 and ς, the parameters associated with shear creep.Among the fifteen parameters,a, M1, l, k, and pr′ in (16) can be obtained from the triaxial compression tests or one-dimensional oedometer tests. ν, E, KG, n, M2, and pr′ in (24) can be determined from the triaxial shear tests. dv, h, and ξ1 can be determined using oedometer creep tests or triaxial compression creep tests. ξ2 and ς can be determined utilizing triaxial shear creep tests.To determine the parametersν, E, KG, n, M2, and pr′ in (24), many triaxial shear tests are required to gain a collection of q′−εs curves, as shown in Figure 9. Gs=3Gi is the initial tangent of q′−εs curves, and Gi is the initial tangent modulus. G can be calculated according to (32). KG can be obtained by drawing the relationship between G/pa and p′/pa in dual logarithmic coordinate system. v=0.3, and E can be determined using (10). M and pr′ can be gained by drawing the CSL in p′−q′ plane. M2 can be evaluated utilizing (33).(32)G=23Gs,(33)M2=MRf0.25,where Rf=qf′/qult′ is the failure ratio, qf′ is the failure stress and is gained from triaxial shear tests, and qult′ is the asymptotic value of q.Figure 9
Schematic stress-strain curves for the triaxial shear test.The parametersl, k, and pr′ in (16) should be obtained from the triaxial compression or one-dimensional oedometer tests. Section 3.2 explains how they arrived at their values. a=0.25−0.15d, d is the slope of the last part of εv−εa curve, as illustrated in Figure 10.Figure 10
Theεa−εv curve.The parametersdv, h, and ξ1 can be determined using oedometer creep tests or triaxial compression creep tests, and their determination has been presented in Section 3.2. The parameters ξ2 and ς can be determined utilizing triaxial shear creep tests, and Section 3.3 gives their detailed calculation process.


## 4. Model Validation and Discussion

The presented model has been implemented in the finite element program Abaqus, and verifications have been carried out by comparing computed results to experimental data.


### 4.1. Drained Triaxial Shear Creep Tests for Tianjin Coastal Clay

Wang et al. [31] carried out drained triaxial shear creep tests for the reconstituted Tianjin Coastal clay. This was a mucky clay whose natural moisture content, plasticity index, density, and permeability coefficient were 40%–67%, 20–23, 16–17 kN/m3, and 5.0 × 10−8 cm/s, respectively. Its specific gravity was 2.64–2.70. In these tests, the confining pressures were 100 kPa, 150 kPa, and 200 kPa. The deviatoric stresses were 60 kPa, 120 kPa, 180 kPa, and 240 kPa. In multistage loading tests, the duration of each load level was 3 d. The parameters are listed in Table 1. The experimental results show that the Tianjin Coastal clay exhibits obvious strain-hardening behavior.Table 1
The parameters for Tianjin Coastal clay.
M1M2alknhdvE (MPa)KGpr′ (kPa)vξ1ξ2ς1.21.220.21543.50.671.850.0743.7255.3860.10.32.723.721.82Figure11 shows the experimental results and the predicted value for the Tianjin Coastal clay at the confining pressure of 150 kPa. As shown in Figure 11(a), the new EVP model can fit the experimental results reasonably well except for the deviatoric stress of 240 kPa. The calculated value is higher than the experimental results in the early stages of creep (0–20 min). Jiang et al. [32] suggested that such problems could be solved by using the hysteretic response equation developed by [33]. To simplify the mathematical structure of the present model, the hysteretic response equation concept is not introduced here.Figure 11
Comparison of the calculated values and the measured values of drained compression creep tests for Tianjin Coastal clay under different stress levels (σ3′=150kPa): (a) axial strains vs. creep time; (b) volumetric strains vs. creep time.
(a)(b)In Figure11(b), the predicted value of volumetric creep strains obtained from the EVP model is compared with the test data for the Tianjin Coastal clay. It is pretty evident that the gas between the predicted value and the measured data grows smaller as deviatoric stress increases. The difference may be due to applying the associated flow rule to the EVP model under low deviatoric stress conditions. Islam et al. [3] demonstrated that the EVP model with a nonassociated flow rule could capture the time-dependent properties of clay better compared to the associated flow type EVP models.


### 4.2. Undrained Triaxial Compression Creep Tests for Shantou Dredger Fill

Luo et al. [34] carried out undrained triaxial compression creep tests for the reconstituted Shantou dredger fill. The Shantou dredger fill property compositions are listed in Table 2. Finally, the diameter and height of specimens were 39.1 mm and 80 mm, respectively. The testing confined pressure changed from 100 kPa to 200 kPa to 300 kPa. The vertical stress σ1 had to change accordingly so as to achieve the deviatoric stress levels D = 0.2, 0.4, 0.6, and 0.8. The test temperature was strictly controlled at 24 ± 1°C during the whole testing period. The parameters are listed in Table 3.Table 2
Fundamental characteristic indexes of the Shantou dredger fill.
Maximum dry densityOptimum water contentDry densityEffective cohesive strengthEffective internal friction angleSpecific gravityLiquid limitPlastic limitPlasticity index(g/cm3)(%)(g/cm3)(kPa)(°)(%)(%)(%)1.6720.11.672328.882.7129.6∼36.113.2∼17.014.9∼19.1Table 3
The parameters for Shantou dredger fill.
Confining pressureM1E (MPa)KGvlkhndvpr′ (kPa)ξ11001.5716.645.30.3112.02.31.360.720.3336.32.452001.7915.252.60.3123.21.51.220.530.1742.43.723002.3312.958.10.3133.71.21.760.690.0755.85.39Figure12 shows good agreement between the simulations using the double-yield nonstationary flow surface EVP model and experiment for undrained triaxial compression creep tests. The developed model generates a larger error when the stress level is small. The deviation between predicted and measured results decreases with the increase in stress level. Overall, the comparisons reveal that the creep behaviors of Shantou dredger fill are sufficient to be reliably predicted.Figure 12
Volumetric strains versus creep time in the experiment and simulation for Shantou dredger fill: (a)σ3′=100kPa; (b) σ3′=200kPa; (c) σ3′=300kPa.
(a)(b)(c)


### 4.3. Drained Triaxial Shear Creep Tests for Silted Soil in Front of Dam

Drained triaxial shear creep tests for silted soil in front of dam were carried out in this paper. The basic properties of silted soil are listed in Table4. Soil samples were collected at a depth of 4-5 m in front of Xiaowan Reservoir, located in the southern mountainous area of Ningxia, and had a diameter of 50 mm and a height of 100 mm. The load applied to soil samples was applied in four stages. The following cell pressures were applied: 200 kPa, 300 kPa. The parameters are listed in Table 5.Table 4
Fundamental characteristic indexes of the silted soil in front of dam.
Moisture content (%)Void ratio (%)Initial density (g/cm3)Saturation (%)Specific gravityLiquid limit (%)Plastic limit (%)Plasticity indexLiquid index25.82.7281.37886.82.7553.126.426.72.3Table 5
The parameters for silted soil in front dam.
Confining pressureM1M2alknhdvE (MPa)KGpr′ (kPa)vξ1ξ2ς2001.281.210.31131.450.491.570.5837.960420.33.62.681.573001.331.330.31241.230.531.480.3130.7253.245.30.32.72.531.32In Figure13, the practical measurement data for silted soil in front of dam are compared with the simulation values for the EVP model. The comparison analysis shows that the EVP model can reasonably characterize damp and stable creep. During accelerated creep, there is a noticeable discrepancy between the estimates and the test data.Figure 13
Comparison of the calculated values and the measured values of drained triaxial creep tests on silted soil in front of dam: (a) axial strain vs. creep time (σ3′=200kPa); (b) volumetric creep strain vs. creep time (σ3′=200kPa); (c) axial strain vs. creep time (σ3′=300kPa); (d) volumetric creep strain vs. creep time (σ3′=300kPa).
(a)(b)(c)(d)


### 4.4. Discussions

A phenomenon was found in (21) and (28); that is, the yield surface changed over time even when viscoplastic strains remained the same. In this sense, the yield surface can be considered an unsteady curve compared to the classical flow surface. The difference between the two flow surfaces is shown in Figure 14. The compression yield function fV = 0 and the shear function fS = 0 define a nonstationary compression surface and a nonstationary shear surface, respectively. All possible stress states in the p′−q′ plane lie in or on the two surfaces. When fV < 0 and fS < 0, the clay lies within the two flow surfaces, and there are only elastic strains. When fV = 0, fS = 0, and a load is applied, both elastic and viscoplastic strains present.Figure 14
The schematic diagram of both nonstationary compression and shear surfaces. For a given viscoplastic strain, the stress point which lies in a flow surface will be reached at different points B1, B2, or B3 over time (t). When soil is seen as an elastoplastic material and viscoplastic strains remain constant, the stress point will be reached at the same point (such as B1) along a given load path, which is independent of time.


## 4.1. Drained Triaxial Shear Creep Tests for Tianjin Coastal Clay

Wang et al. [31] carried out drained triaxial shear creep tests for the reconstituted Tianjin Coastal clay. This was a mucky clay whose natural moisture content, plasticity index, density, and permeability coefficient were 40%–67%, 20–23, 16–17 kN/m3, and 5.0 × 10−8 cm/s, respectively. Its specific gravity was 2.64–2.70. In these tests, the confining pressures were 100 kPa, 150 kPa, and 200 kPa. The deviatoric stresses were 60 kPa, 120 kPa, 180 kPa, and 240 kPa. In multistage loading tests, the duration of each load level was 3 d. The parameters are listed in Table 1. The experimental results show that the Tianjin Coastal clay exhibits obvious strain-hardening behavior.Table 1
The parameters for Tianjin Coastal clay.
M1M2alknhdvE (MPa)KGpr′ (kPa)vξ1ξ2ς1.21.220.21543.50.671.850.0743.7255.3860.10.32.723.721.82Figure11 shows the experimental results and the predicted value for the Tianjin Coastal clay at the confining pressure of 150 kPa. As shown in Figure 11(a), the new EVP model can fit the experimental results reasonably well except for the deviatoric stress of 240 kPa. The calculated value is higher than the experimental results in the early stages of creep (0–20 min). Jiang et al. [32] suggested that such problems could be solved by using the hysteretic response equation developed by [33]. To simplify the mathematical structure of the present model, the hysteretic response equation concept is not introduced here.Figure 11
Comparison of the calculated values and the measured values of drained compression creep tests for Tianjin Coastal clay under different stress levels (σ3′=150kPa): (a) axial strains vs. creep time; (b) volumetric strains vs. creep time.
(a)(b)In Figure11(b), the predicted value of volumetric creep strains obtained from the EVP model is compared with the test data for the Tianjin Coastal clay. It is pretty evident that the gas between the predicted value and the measured data grows smaller as deviatoric stress increases. The difference may be due to applying the associated flow rule to the EVP model under low deviatoric stress conditions. Islam et al. [3] demonstrated that the EVP model with a nonassociated flow rule could capture the time-dependent properties of clay better compared to the associated flow type EVP models.


## 4.2. Undrained Triaxial Compression Creep Tests for Shantou Dredger Fill

Luo et al. [34] carried out undrained triaxial compression creep tests for the reconstituted Shantou dredger fill. The Shantou dredger fill property compositions are listed in Table 2. Finally, the diameter and height of specimens were 39.1 mm and 80 mm, respectively. The testing confined pressure changed from 100 kPa to 200 kPa to 300 kPa. The vertical stress σ1 had to change accordingly so as to achieve the deviatoric stress levels D = 0.2, 0.4, 0.6, and 0.8. The test temperature was strictly controlled at 24 ± 1°C during the whole testing period. The parameters are listed in Table 3.Table 2
Fundamental characteristic indexes of the Shantou dredger fill.
Maximum dry densityOptimum water contentDry densityEffective cohesive strengthEffective internal friction angleSpecific gravityLiquid limitPlastic limitPlasticity index(g/cm3)(%)(g/cm3)(kPa)(°)(%)(%)(%)1.6720.11.672328.882.7129.6∼36.113.2∼17.014.9∼19.1Table 3
The parameters for Shantou dredger fill.
Confining pressureM1E (MPa)KGvlkhndvpr′ (kPa)ξ11001.5716.645.30.3112.02.31.360.720.3336.32.452001.7915.252.60.3123.21.51.220.530.1742.43.723002.3312.958.10.3133.71.21.760.690.0755.85.39Figure12 shows good agreement between the simulations using the double-yield nonstationary flow surface EVP model and experiment for undrained triaxial compression creep tests. The developed model generates a larger error when the stress level is small. The deviation between predicted and measured results decreases with the increase in stress level. Overall, the comparisons reveal that the creep behaviors of Shantou dredger fill are sufficient to be reliably predicted.Figure 12
Volumetric strains versus creep time in the experiment and simulation for Shantou dredger fill: (a)σ3′=100kPa; (b) σ3′=200kPa; (c) σ3′=300kPa.
(a)(b)(c)


## 4.3. Drained Triaxial Shear Creep Tests for Silted Soil in Front of Dam

Drained triaxial shear creep tests for silted soil in front of dam were carried out in this paper. The basic properties of silted soil are listed in Table4. Soil samples were collected at a depth of 4-5 m in front of Xiaowan Reservoir, located in the southern mountainous area of Ningxia, and had a diameter of 50 mm and a height of 100 mm. The load applied to soil samples was applied in four stages. The following cell pressures were applied: 200 kPa, 300 kPa. The parameters are listed in Table 5.Table 4
Fundamental characteristic indexes of the silted soil in front of dam.
Moisture content (%)Void ratio (%)Initial density (g/cm3)Saturation (%)Specific gravityLiquid limit (%)Plastic limit (%)Plasticity indexLiquid index25.82.7281.37886.82.7553.126.426.72.3Table 5
The parameters for silted soil in front dam.
Confining pressureM1M2alknhdvE (MPa)KGpr′ (kPa)vξ1ξ2ς2001.281.210.31131.450.491.570.5837.960420.33.62.681.573001.331.330.31241.230.531.480.3130.7253.245.30.32.72.531.32In Figure13, the practical measurement data for silted soil in front of dam are compared with the simulation values for the EVP model. The comparison analysis shows that the EVP model can reasonably characterize damp and stable creep. During accelerated creep, there is a noticeable discrepancy between the estimates and the test data.Figure 13
Comparison of the calculated values and the measured values of drained triaxial creep tests on silted soil in front of dam: (a) axial strain vs. creep time (σ3′=200kPa); (b) volumetric creep strain vs. creep time (σ3′=200kPa); (c) axial strain vs. creep time (σ3′=300kPa); (d) volumetric creep strain vs. creep time (σ3′=300kPa).
(a)(b)(c)(d)


## 4.4. Discussions

A phenomenon was found in (21) and (28); that is, the yield surface changed over time even when viscoplastic strains remained the same. In this sense, the yield surface can be considered an unsteady curve compared to the classical flow surface. The difference between the two flow surfaces is shown in Figure 14. The compression yield function fV = 0 and the shear function fS = 0 define a nonstationary compression surface and a nonstationary shear surface, respectively. All possible stress states in the p′−q′ plane lie in or on the two surfaces. When fV < 0 and fS < 0, the clay lies within the two flow surfaces, and there are only elastic strains. When fV = 0, fS = 0, and a load is applied, both elastic and viscoplastic strains present.Figure 14
The schematic diagram of both nonstationary compression and shear surfaces. For a given viscoplastic strain, the stress point which lies in a flow surface will be reached at different points B1, B2, or B3 over time (t). When soil is seen as an elastoplastic material and viscoplastic strains remain constant, the stress point will be reached at the same point (such as B1) along a given load path, which is independent of time.


## 5. Conclusions

In this paper, a new EVP model is formulated to describe the stress-strain-time properties of soft clay. The following are the main characteristics of the EVP model with double-yield surfaces:(1)
Plastic strains are used as the hardening parameters for both yield surfaces in the double-yield nonstationary flow surface EVP model.(2)
The model can model the compression and shear creep of soft clay, including decay creep and steady creep.(3)
Compared with Yin’s double-yield surface mode, seven more parameters need to be defined for the accurate description of both shear creep and compression creep using triaxial rheological tests. The rest of the parameters can be determined by general triaxial tests.To better describe soft clay’s time-dependent behaviors, subsequent improvements to the EVP model with the nonassociated flow rule should be promoted.



---

*Source: 5892890-2022-05-23.xml*