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

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. dimensional oedometer tests. , E , K G , n , M 2 , and p r ′ can be determined from the triaxial shear tests. d v , 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.


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. erefore, 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]. e 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 di erential approach [11] and engineering theories of creep [12]; (3) general stress-strain-time models, such as overstress models [13,14] and nonstationary ow surface (NSFS) models [15]; (4) other contributions, such as Borja's model [16] and subloading surface model [17]. e rst two categories can help further the phenomenological understanding of the creep properties of soft clay. erefore, the rst 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 speci c boundaries and loading conditions [1]. On the contrary, the general stress-strain-time models can consider time e ects and in-viscid deformation under any possible loading conditions, which prompts them to be widely used, especially in finite element analysis. However, most stressstrain-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][19][20].
In this paper, a practical model for the time-dependent behaviors of soft clay is developed. e model's main characteristics are as follows: ( 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. e presented model is then put into the finite element program Abaqus and validated with a variety of soft clays.

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.
where _ ε is the elastic strain rates, _ ε vp ij is the viscoplastic strain rates, _ ε vp vij is the volumetric creep strain rates, and _ ε vp sij is the shear creep strain rates. e 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. e elastic strain rates _ ε e ij in (1) are supposed to obey generalized Hooke's law and can be calculated: where C ijkj 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 E e and the Poisson ratios v e . e viscoplastic strain rates _ ε vp ij in (1) are defined according to the nonstationary flow surface theory presented by Olszak and Perzyna [13]. Flow surface can be represented as follows: where F is the loading function of σ ij ′ , K s , and β; σ ij ′ is the current effective stress; K S is a strain-hardening parameter related to the viscoplastic strains ε vp ij ; and β is a time variable associated with clay proper. _ ε vp ij can be determined based on the assumption of associative normality flow rule, where the loading function F is equal to the potential function g.

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 matrix C ijkj in (3) can be written as follows: 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 2 Advances in Civil Engineering clay, Poisson's ratio usually equals 0.3, due primarily to the smaller proportion of elastic strain in total strain [26]).
where K G and n are material parameters, and p a is the atmospheric pressure and is introduced to make parameters dimensionless. e relationship of both G/p a and p/p a is a straight line in dual logarithmic coordinate system, and the intercept and slope of the line are K G and n, respectively.

Compression
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 p c ′ is the hardening function taking the plastic volume strain ε p v as hardening parameter. e yield curve of the modified Cam-clay model is an ellipse in the p ′ − q ′ plane, as indicated in Figure 1. e 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.
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 ′ + p r ′ extended modified Cam-clay model to the condition that cohesion force exists. Equation (11) can be rewritten as follows: where is the hardening function; like p c ′ , p r ′ 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]. e yield curve of the modified model is depicted in Figure 2.
e 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.
For the modified Cam-clay model, the e− lnp relationship is available in the construction of p c ′ − ε p v relation. However, the relationship is not linear when p ′ value is small, especially for compacted soils, as shown in Figure 3. e linear relation between e and lnp will produce certain difficulties in mathematical treatment. erefore, it is a necessity to modify the approach for the modified Cam-clay model [29] formulated the relationship between p c ′ and ε v by a hyperbolic function on the basis of experiment analyses.

Advances in Civil Engineering
isobaric conditions, where the paper takes the atmospheric pressure p a for model parameters dimensionless. e relation between p c ′ /p a and ε v shown in Figure 4(a) is nearly linear in p c ′ /p a ε v − p c ′ /p a coordinates, as depicted in Figure 4(b). e p c ′ /p a − ε v curve can be expressed as follows: where k is the slope of the straight line, l � 1/ε a v is the intercept of the straight line on the horizontal axis, ε a v is the asymptotic value of ε v , andE i is the initial tangent slope of the hyperbola.
e relationship between ε p v and p c ′ 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. en, (13) can be also written as follows: e strain-hardening function p c ′ can be expressed as follows: By substituting (15) into (13) and replacing M in (13) with M 1 , we can obtain where M 1 � 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 ε t v over time, as indicated in Figure 5. It is important to note that the timedependent volumetric strain ε t v does not include the instantaneous volume strain ε p v . Many scholars describe the ε t v − t relationship with a hyperbola. e hyperbolic function can be given as follows:   where c v � 1/h, h is the initial tangent slope of the hyperbola, Supposing that the ratio of ε u ] and ε p v is constant under different stress conditions, i.e., ε u v /ε p v � ξ 1 , ε t v can be further expressed by h and ξ 1 as follows: After rearranging (17), the expression for ε p v can be obtained: By substituting (19) into (16), the following equation can be obtained: en, the compression creep surface function can be given as follows: where H V is the strain-hardening function related to the viscoplastic volumetric strains ε t v , 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.

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 ε p s is selected as the hardening parameter of the shear creep yield criterion.
Hsieh et al. (1990) recommended a hardening law for shear creep based on triaxial tests for clay as follows: where G is the elastic shear modulus, M 2 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 doubleyield surface model [26], as indicated in Figure 6. Triaxial shear creep experiments are done to obtain the variation of the shear strain ε t S over time, as indicated in Figure 7. A hyperbolic function for ε t S and t can be proposed as follows: where c s � 1/ς is a material parameter, ς is the initial tangent slope of the hyperbola, d s � 1/ε u s is a material parameter, and ε u s is the asymptotic value of ε t S . It is reasonable to suppose that the ratio of ε u s and ε p s is constant under different stress conditions; i.e., ε u s /ε p s � ξ 2 . ε t s can be further expressed by ς and ξ 2 as follows: After rearranging (26), the expression for ε p S can be obtained: By substituting (27) into (24), the shear creep surface function can be given as follows: where H S is the strain-hardening function related to the viscoplastic shear strain ε t s , and β S is the time variable associated with soil property. e constitutive model is also formulated by applying the associated flow theory. e compression creep function FV and the shear creep function FS are illustrated in Figure 8. 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.
us, the EVP model, which is a comprehensive model used Advances in Civil Engineering 5 in soft clay, can simulate the dilatancy and time-dependent behaviors of soft clay.

Viscoplasticity Flexibility Matrix.
It can be concluded from (7) that the viscoplastic flexibility matrix C tp is key to determining the viscoplastic strain rate _ ε vp ij ; how to acquire C tp will be detailed in this section.
Based on the above analysis, C tp can be further written as follows: where A � zF V /zε t v , B � zF S /zε t s . e different terms in (30) can be expressed as follows:

Parameter Determination.
ere are fifteen parameters in the proposed model, including ], Poisson's ratio; E, elastic modulus; K G and n, elastic material parameters; M 1 and M 2 , modified strength parameters; p r ′ , 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; d v , h, and ξ 1 , the parameters related to compression creep; and ξ 2 and ς, the parameters associated with shear creep.
Among the fifteen parameters, a, M 1 , l, k, and p r ′ in (16) can be obtained from the triaxial compression tests or one-Advances in Civil Engineering dimensional oedometer tests. ], E, K G , n, M 2 , and p r ′ in (24) can be determined from the triaxial shear tests. d v , 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, K G , n, M 2 , and p r ′ in (24), many triaxial shear tests are required to gain a collection of q ′ − ε s curves, as shown in Figure 9. G s � 3G i is the initial tangent of q ′ − ε s curves, and G i is the initial tangent modulus. G can be calculated according to (32). K G can be obtained by drawing the relationship between G/p a and p ′ /p a in dual logarithmic coordinate system. v � 0.3, and E can be determined using (10). M and p r ′ can be gained by drawing the CSL in p ′ − q ′ plane. M 2 can be evaluated utilizing (33).
where R f � q f ′ /q ult ′ is the failure ratio, q f ′ is the failure stress and is gained from triaxial shear tests, and q ult ′ is the asymptotic value of q. e parameters l, k, and p r ′ 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.15 d, d is the slope of the last part of ε v − ε a curve, as illustrated in Figure 10. e parameters d v , 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. e parameters ξ 2 and ς can be determined utilizing triaxial shear creep tests, and Section 3.3 gives their detailed calculation process.

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

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. is was a mucky clay whose natural moisture content, plasticity index, density, and permeability coefficient were 40%-67%, 20-23, 16-17 kN/m 3 , 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. e 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. e parameters are listed in Table 1. e experimental results show that the Tianjin Coastal clay exhibits obvious strain-hardening behavior. Figure 11 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. e 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.
In Figure 11(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. e 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.

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. e 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. e testing confined pressure changed from 100 kPa to 200 kPa to 300 kPa. e 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. e test temperature was strictly controlled at 24 ± 1°C during the whole testing period. e parameters are listed in Table 3. Figure 12 shows good agreement between the simulations using the double-yield nonstationary flow surface EVP model and experiment for undrained triaxial compression creep tests. e developed model generates a larger error when the stress level is small. e 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.

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. e basic properties of silted soil are listed in Table 4. 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. e load applied to soil samples was applied in four    stages. e following cell pressures were applied: 200 kPa, 300 kPa. e parameters are listed in Table 5.
In Figure 13, the practical measurement data for silted soil in front of dam are compared with the simulation values for the EVP model. e 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.

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. e difference between the two flow surfaces is shown in Figure 14. e compression yield function f V � 0 and the shear function f S � 0 define a nonstationary compression surface and a nonstationary shear      Advances in Civil Engineering surface, respectively. All possible stress states in the p ′ − q ′ plane lie in or on the two surfaces. When f V < 0 and f S < 0, the clay lies within the two flow surfaces, and there are only elastic strains. When f V � 0, f S � 0, and a load is applied, both elastic and viscoplastic strains present.

Conclusions
In this paper, a new EVP model is formulated to describe the stress-strain-time properties of soft clay. e 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. To better describe soft clay's time-dependent behaviors, subsequent improvements to the EVP model with the nonassociated flow rule should be promoted.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.