Extension of a basic hypoplastic model for overconsolidated clays

This paper presents a hypoplastic constitutive model for overconsolidated clays. The model is developed based on a basic model proposed recently for sand. New density and sti(cid:11)ness factors are introduced to account for history dependence. The model is capable of predicting a desirable asymptotic state boundary surface for overconsolidated clays. Validation with experimental results shows this model can properly describe the main features of overconsolidated clays.


Introduction
Constitutive modeling of overconsolidated (OC) clays has a long history, pioneered by the Modified Cam-clay (MCC) model based on the critical state concept [32,33]. The MCC model is capable of describing the behavior of normally consolidated and lightly OC clays. The capacity of the model for OC clays has been largely enhanced by introducing various concepts, such as bounding surface [6] and subloading surface [12]. These theories can be integrated with various dilatancy relations to form constitutive models for overconsolidated soils [7,47,48,49]. An inspection of the relevant literature reveals, however, that most of the existing models for OC clays are based on elastoplasticity.
Hypoplastic models are developed as an alternative approach to the prevailing elastoplastic models in describing the overconsolidation behavior. Following the pioneer work by Kolymbas and Wu [17,18,40,41], different versions of hypoplastic models were developed for sand [4,9,37] . These models served as the basis for the model development of clay. The first attempt to extend hypoplastic model for clay was made by Niemunis [26,27], who combined the critical state concept with hypoplasticity to describe the behavior of OC clays. Meanwhile, Herle and Kolymbas [14] modified the model by von Wolfferdorff [37] (VW model) for soils with low friction angle. Based on these works, Mašín [21] proposed the first hypoplastic model for clay (DM model), which employs a density factor to describe the consolidation history.
However, later Mašín himself found that the DM model gives rise to an unrealistic asymptotic state boundary surface (ASBS) when a certain combination of parameters is used [23]. The reason for this deficiency is that the adopted density factor is not only dependent on the consolidation history but also related to other material parameters. To address this issue, Mašín [22,23] proposed a new form of hypoplastic models with explicit ASBS. Nevertheless, the original form is still valuable for further developments and enhancements. For instance, Wu et al. [45] reported a basic hypoplastic model for sand. This model employs the basic form of previous hypoplastic models but additionally includes a new tensorial term, which provides an opportunity to refine the modelling of the dilatant behaviour and asymptotic state.
In this work, we will extend the basic model [45] for OC clays. In the following, we first inspect the limitation of the basic model. Next, a detail process for developing the model for OC clays is presented. A new density factor is introduced to account for consolidation history. The employment of the density factor allows the model to predict reasonable ASBSs for OC clays. The characteristics of the model are then investigated by a parametric study. At the end, model performance is demonstrated by simulating several element tests, and comparing the prediction with experiments on different OC clays.

General framework
To gain perspective, we start with the rate-form constitutive equation by Wu and Kolymbas [41], which has the following general form: where L and N are fourth-and second-order tensors, respectively. D is the strain rate (stretching) tensor. D = tr(D 2 ) stands for the norm of the strain rate tensor.T is the Jaumann stress rate defined by: whereṪ is the material time-derivative of the Cauchy stress tensor and W the spin tensor. Bear in mind that the failure criterion and flow rule of the model (1) is the interactional outcome of L and N. This implies the intensity and direction of flow are mutually dependent. To remedy this shortcoming, Niemunis [26] proposed a simple rearrangement of the constitutive framework (1) by introducing a tensor function B: where the tensor B includes both the flow direction and intensity of the flow by the function Y and a tensor m, respectively.
With the above rearrangement, a generalized form of hypoplastic models is obtained as follows: where f s and f d are the stiffness and density factors, accounting for the effects of stress level and density, respectively [21,26].

A basic hypoplastic model
To develop a new constitutive model for OC clays, we start by considering the following basic model proposed by Wu et al. [45]: in which C i (i = 1, 2, 3, 4) are material constants. T * is the deviatoric stress tensor defined as T * = T − I/3trT with I being the unit tensor. The material constants can be related to some well established parameters in soil mechanics, such as the initial tangential modulus E i , the initial Poisson's ratio ν i , the critical friction angle φ c , and dilatancy angle ψ. Moreover, for critical state with ψ = 0, the model (5) can be further rearranged as the following form: whereT is the normalized stress defined asT = T/tr(T); the multiplier c 2 permits a refined modelling of the dilatancy behavior. The constant a accounts for the critical state through: Note that the equation (6) is a special case of equation (5) with ψ = 0. This basic model has several merits such as simplicity and robustness in numerical simulation. These advantages make it an interesting stand-alone model for many applications [30,31,36,38,46] and further developments [19,20,39].

Limitation of the basic model for OC clays
The performance of the basic model is well documented in the literature [45]. The model is capable of modeling the basic behavior of normally consolidated clay, such as failure and contractive deformation. However, there exist some shortcomings for modeling OC clays. The most important ones are summarized as follows: • The consolidation history is not incorporated in this model, thus it cannot describe strain softening and volume dilatancy of OC clays.
• This model does not include the swelling index κ * . Therefore, the directional stiffness for loading and unloading is not adjustable.
• This model has an constant stiffness factor f s . It is adequate to describe the initial shear modulus of normally consolidated soil but it will underestimate significantly the initial shear moduli of OC clays.
Note that the VW model and DM model have similar forms with model (6) except that the second term tr(D)T is not included. These two models address some of the above mentioned shortcomings. However, the VW model is seriously flawed because the underlying basic equation does not perform properly [45]. The DM model also suffers from some shortcomings [23].

Proposed model for OC clays
To address the aforementioned shortcomings of the basic model (6), we proceed to extend it by considering some salient features of OC clays.

Model framework
We first consider the linear behavior of the model (6) by manipulating the linear terms. To doing so, the model (6) can be expressed as follows: where L i (i = 1, 2, 3) are the linear tensorial function in equation (6) with c i (i = 1, 2, 3) being the material constants. The effects of the linear terms can be best appreciated by considering the response envelope. For a rosette of strain rates with a constant magnitude, the response envelope at a given initial stress state T 0 is the surface spanned by all stresses T = T 0 + ∆T [10]. Figure 1 shows some response envelopes depicted on the Rendulic plane with T 22 = T 33 . the response envelope of the model (6) represents an ellipse. Obviously, the multiplier c 1 scales the magnitude of D, and thus can expand or shrink the response envelopes. The multipliers c 2 and c 3 have similar effects, that is increase c 2 and c 3 gives rise to longer major axis while does not change the minor axis of the ellipse.
Note that the distance from the initial stress T 0 to any point on the response envelope represents the stiffness along certain strain rate. As shown in Figure 1(a), the vectors | oa| and | ob| represent, respectively, the directional stiffness of loading and unloading at a isotropic stress state. | oc| and | od| denote the directional stiffness with isochoric strain. For clays, the ratio | oa|/| ob| is specified by the loading and swelling indexes, i.e., κ * and λ * . It is known that a desirable ratio of | oa| to | oc| can be achieved by implementing two scale multipliers with L 1 and L 3 [14]. Analogously, with the help of the term c 2 L 2 , an adjustable ratio of loading to unloading stiffness is available by introducing the swelling index κ * .
Another remark concerns with the enhancement of the initial shear stiffness, G i , for OC clays. This subject has been investigated by several experiments [1,2]. The experimental results indicate that the dependence of G i on stress level and OCR can be described by exponential functions [11]. However, hypoplastic models usually underestimate G i of OC clays, which is an intrinsic deficiency of many hypoplastic models. As  Figure 1, the initial shear stiffness can be enhanced by multiplying the whole model. Following this way, we introduce an additional multiplier into the model where m R represents the stiffness enhancement of OC clays. In the following, the elements for the model (9) are derived. The complete formulations are given in the appendix.

Tensor function L
The linear tensorial function of model (6) is considered for the new model. It can be written as: where I and 1 are the fourth-order and the second-order unit tensor, respectively. Next, we proceed to determine the constants c 1 , c 2 , c 3 by assuming that the proposed model has the same form with the model (6) at an isotropic stress state. Therefore, the model (9) can be expressed with the following decomposed forms: whereṗ andq are the rates of the mean stress and the deviatoric stress, respectively. q = q/tr(T) is the normalized deviatoric stress. D v and D q are the volumetric strain and deviatoric strain, respectively. For an isotropic stress state, we have initiallyq = 0, D q = 0, and D v < 0, and thus equation (11a) becomes: where

3/3)/3 represents the initial bulk modulus of the soil for loading.
Since the modification of the model does not influence the formulation of the model (6) at an isotropic stress state, and therefore we have: For unloading in an isotropic compression test, we have D v > 0. The rate of mean stress can be gained as follows:

3/3)/3 represents the initial bulk modulus of the soil for unloading.
On the other hand, in a drained test with constant volume change, i.e., D v = 0 (isochoric shear), the equation (11b) is recast to: wherein G i = 0.5f s c 1 is the initial tangential shear modulus for isochoric shear. Comparison between the initial tangential bulk modulus in loading and unloading conditions yields In addition, the ratio of K + i to G i can be considered as a material constant v i . Combining the equations (12 to 15) and making use of the above conditions, we obtain the formulations for the multipliers c 1 , c 2 : with γ being the ratio (λ * + κ * )/(λ * − κ * ). The third constant c 3 is then obtained from equation (13).
3.3. Flow rule m and limit stress condition Y The hypoplastic flow rule m indicates the flow direction, which should obey the conditions m ∼T * at a critical state, and m = I/ √ 3 at an isotropic stress state [26]. The tensorial function B of the model (6) fulfils these requirements. In this paper, a similar form including Matsuoka-Nakai failure criterion is adopted: where F is an interpolation function to account for the Matsuoka-Nakai failure criterion [21,37]. The failure criterion is obtained by incorporating the Matsuoka-Nakai failure criterion into the non-linearity Y . In this work, a similar Y as that used in [21] is adopted. It reads: 3.4. The density factor f d The factor f d is used to consider the consolidation history. Bearing in mind that f d = 1 shall be guaranteed for normal consolidation (OCR = 1), we proceed to introduce the following simple expression for the density factor: where α is a model parameter; R is the stress ratio representing the overconsolidation [47]. A smaller value of R corresponds to a larger OCR to give: in which the p + e and p * e are the pre-consolidation pressure in MCC model and the Hvorslev equivalent pressure respectively. Their definitions are shown in Figure 2 with the formulations given as follows: in which the inclination M corresponds to the stress ratio at the critical state. It is interpolated with the function F to define the critical state line at different Lode angles θ, i.e., M = 6F sinφ c /(3−sinφ c ); N stands for the logarithmic value of the specific volume at the reference stress p r = 1 kPa; the index λ * is the slope of the normal compression line in the double logarithmic ln(1 + e) − lnp plane. The stiffness factor f s The factor f s can be obtained from the pre-defined isotropic normal compression line. After Butterfield [5], the normal compression line is expressed as: Time differentiation of equation (22) gives rise to the rate form: Comparing equations (23) and (12), and making use of D v =ė/(1+e), the expression for f s is obtained: It can be seen that the stiffness factor f s includes both the strength parameter a and the deformation parameters (λ * , κ * ) through a simple formulation.
3.6. The multiplier m R In this paper, we propose a simple formulation for m R : where β is a constant controlling the influence of consolidation history on the initial shear modulus. Obviously, we have m R = 1 for normal consolidation with R = 1; for extreme case with R → 0, the multiplier m R = exp(β). If needed, β can be considered as a material parameter. In the present simulation, however, a fixed value β = 1 is used.

Model performance
In this section, we discuss the basic properties of the proposed model. First, a parametric study is performed to facilitate the selection of proper parameters in numerical simulation. Next, the performance of the proposed model is presented. Additionally, the influence of parameters on the shape and size of the ASBS is investigated.

Model parameters
The proposed hypoplastic model contains six constitutive parameters, i.e., φ c , N , λ * , κ * , v i , and α. The first five have the same physical interpretation as those used in the MCC model: φ c is the critical state friction angle; N is the value of ln(1 + e) at the isotropic normal compression line for p = 1 kPa; λ * is the slope of the isotropic normal compression line on the ln(1 + e) − lnp plane; κ * denotes the slope of unloading line in the same plane, and v i is ratio between the initial bulk and shear moduli. Since the first five parameters of our model have the same physical meaning as those in the DM model. In the following, therefore, a series of triaxial tests are simulated to show the effects of α. In the simulations, two different conditions with OCR = 5 and 1.2 are considered. The confining pressure is assumed to be 100 kPa. The parameters used in the simulations are listed in Table 1.  Figure 3 shows the stress -strain relationship and effective stress paths of the model. A close inspection of the results reveals that α mainly influences the peak friction angle  Figure 3(a, c), The parameter α governs the rate to achieve the critical state; a smaller α leads to a faster convergence to the critical state. In addition, α influences the effective stress paths in triaxial tests. As shown in Figure 3(b, d), the normalized stress paths expand out with the increase of α in triaxial compression tests. However, this influence is less significant in triaxial extension tests. Note that in the DM model α is a variable related to λ * , κ * , and φ c . In a recent model by Mašín [23], a fixed value α = 2 is adopted for numerical simulations.

Influence of parameters on the ASBS
Asymptotic behaviors are among the most significant mechanical properties of clays. For instance, one of the most important asymptotic states of clay is the well -known critical state behavior, which describes the asymptotic behavior of clay with constant volume during shearing [22,23,24]. The asymptotic states of a hypoplastic model can be achieved after a sufficiently long proportional stretching, i.e. stretching with a constant direction of the strain rate. This leads to the ASBS that bounds all admissible asymptotic states in the stress -void ratio space.
The size of ASBS should be controlled by the critical state friction angle φ c . As shown in Figure 4, the ASBS of the proposed model expands with the increase of φ c , This expansion does not leads the ASBS to the tensile stress region. It is found that the ASBS of our model is close to the elliptic shape of the yield surface of the MCC model [32] in the range p > p * e /2. The elliptic shape, however, becomes narrow in the range p < p * e /2, coinciding with the nonlinear Hvorslev surface on the dry side. This property ensures the model predicts properly the peak strength of heavily OC clays [3,34].  [24] reported that the ratio λ * /κ * has a significant influence on the shape of the ASBS predicted by the DM model. As shown in Figure 5(a), for high ratios of λ * /κ * , the DM model predicts unrealistic shape of the ASBS. However, this shortcoming does not exist in the proposed model. As shown in Figure 5(b), the proposed model predicts very similar ASBS with varying ratios of λ * /κ * , even with high values of λ * /κ * the shape of ASBS is reasonable.

Mašín and Herle
For the DM model, the dependence of ASBS on the ratio of λ * /κ * is due to the connection between λ * /κ * and the parameter α. As α is not an independent parameter in the DM model, high ratio of λ * /κ * will give rise to small value of α. In the proposed model, however, the parameter α is given as an independent material constant, and thus the ratio of λ * /κ * does not influence α. As shown in Figure 6, increase of α gives rise to expansion in the wet side (p/p e > 0.5), and shrinkage in the dry side (p/p e < 0.5). This is consistent with the results presented in Figure 3, that the increase of α leads to increased peak friction angle of OC clays, and this influence becomes more obvious with low value α. However, the shape of ASBS remains almost unchanged with α increasing from 1.0 to 3.0.

Validation with experimental data
In this section, model predictions are compared with experimental results on both remoulded and undisturbed OC clays. Different types of triaxial tests, including isotropically consolidated undrained triaxial tests (CIU), isotropically consolidated drained triaxial tests with constant mean pressure (CIDP) and constant confining pressure (CID), are simulated to validate this model. Parameters used in the simulations are summarised  Table 2. The first five parameters are collected from the literature, while the parameters v i and α are obtained through optimization methods [16,50,51] to gain the best match with experimental results.

CIU test, Remoulded Wenzhou clay
The Wenzhou clays are mainly marine deposits from Wenzhou, eastern coastal China. The deposit usually reaches a depth of 150 m and exhibits evident heavily OC behavior. A series of triaxial compression tests on remoulded clay was carried out by Gu et al. [8] using a GDS device. The samples contain clay fraction of 63.8% and a fine fraction of 93.4%. The natural water content can reach to w L = 67.7% with plasticity index I p = 37.2 and initial void ratio, e 0 , ranging from 1.55 -1.59 [8].  Figure 7: Comparison between experimental and numerical undrained triaxial tests on OC Wenzhou clay. Experimental data from Gu et al. [8] Four samples were initially subjected to pre-consolidation pressure (p c )of 100, 200, 400, 800 kPa. After the completion of pre-consolidation, the effective confining pressure was reduced to p 0 = 100 kPa in order to produce OC samples with specific OCR values. The test results are presented in Figure 7 together with model prediction. As shown in Figure 7(a), the samples with different OCRs exhibit various shear strengths, this history dependence can be predicted by the proposed model. The peak strength can also be captured quite accurately. However, the stress -strain curves agree not very well at the bends of the curves, as it underestimates the initial stiffness and also the degradation rate of the stiffness. Figure 7(b) shows that the predicted effective stress paths agree well with experimental results for heavy consolidation, while less well prediction is observed for normal consolidation.

CIDP test, Remoulded Fujinomori clay
Experimental data on reconstituted Fujinomori clay by Nakai and Hinoko [25] is a classical set of experiments used in the literature to evaluate different constitutive models within the framework of elastoplasticity [47]. The test was carried out with constant mean effective stress. Both drained triaxial compression and extension tests were performed on OC Fujinomori clay with OCR from 1 to 8; The sample with OCR = 8 was tested at mean effective stress of p 0 = 98 kPa, while other samples were tested at mean effective stress of 196 kPa.
The predicted stress -strain and volumetric change curves using the proposed hypoplastic model are shown in Figure 8. Good agreement in the stress ratio curves is achieved in the compression tests, while the model underestimates the peak strength in the extension tests. A more significant difference is observed in the results of volumetric strain, especially in the extension tests. Nevertheless, the main features of the variation of stress ratio and volumetric strain can be reasonably captured by the proposed model.

CIU and CID test, Remoulded Weald clay
A series of tests have been carried out by Parry [29] on reconstituted Weald clay with both normal consolidated (p c = 206.8 kPa, OCR = 1) and heavy consolidated (p c = 827.4 kPa, OCR = 12) states. The experiments are simulated and the model prediction are compared with experimental data. Parameters for this simulation are collected from the literature [23] by Mašín, who compared his two hypoplastic models for clay using this set of data.
Predictions of the drained and undrained triaxial compression and extension tests, together with experimental results, are shown in Figure 9. It can be seen that the model predictions agree very well with the experimental results. Especially, excellent agreements are achieved in the stress -strain curves of the drained tests. On the other hand, this model underestimates the shear strength in undrained triaxial extension test for OCR = 12.

CID and CIU test, undisturbed Bangkok stiff clay
In the previous evaluation, the model is validated with different types of reconstituted clays. This model is able to capture most properties of reconstituted OC clay. To  Figure 8: Comparison between experimental and numerical drained constant mean pressure triaxial compression and extension tests on OC Fujinomori clay, data from Nakai and Hinoko [25] broaden the validation of the proposed model, the final set of modeling will focus on an undisturbed clay, Bangkok stiff clay, from Bong Ngoo Hao site and Rangsit site. According to the original test results [13], although taken from different sites, these two soils exhibit similar overconsolidation behavior. In this study, selected triaxial compression test are simulated using the proposed model. The reported peak friction angle and cohesion for samples from Bong Ngoo Hao site are φ p = 26 • and c = 30 kPa, respectively. This friction angle is apparently higher than the critical friction angle. Considering the high cohesion value and taking into account that the failure envelope for stiff, OC clays is curvilinear, a critical friction angle φ c = 31.5 • is adopted for the Bong Ngoo Hao stiff clay. Results from drained and undrained triaxial tests on these samples show that a slightly lower critical friction angle was obtained for samples from the Rangsit site. These results are adopted in our simulation. The parameters λ * and κ * are calibrated from oedometer tests presented in [35], and the other parameters are listed in Table 1.
Comparison between numerical and experimental results of the drained compression tests on undisturbed Bangkok clay from Bong Ngoo Hao site is presented in Figure  − Figure 9: Comparison between experimental and numerical drained (a,b) and undrained (c,d) triaxial tests on remoulded Weald clay. Experimental data from Parry and Aust [29] 10. The adopted pre-consolidation pressure for the simulations is 700 kPa with various confining pressures of 35 Figure 11. Note that the results for confining pressures of 140, 276, and 345 kPa are not shown in the p/q − ε 1 curves because these testing results showed significant deviation from the general experimental trend. The corresponding effective stress paths are repented in the p − q curves with solid symbols. An excellent agreement with experimental results is achieved. This again suggests that the proposed model can  Figure 10: Comparison between experimental and numerical drained triaxial compression tests on undisturbed Bangkok stiff clay at Bong Ngoo Hao site. Experimental data from Hassan [13] correctly predict the behaviors of OC clays.

Conclusions
In this paper, we present extend a basic hypoplastic constitutive model for OC clays. New density and stiffness factors are introduced to account for consolidation history. The model is capable of predicting a desirable ASBS. The shape of the ASBS is close to the elliptic shape of the yield surface of the MCC model on the wet side, while similar to the nonlinear Hvorslev surface on the dry side. This feature ensures that the model can properly predict the peak and critical strengths of OC clays. This model is validated by experiments on both reconstituted and undisturbed OC clays. The validation shows that the proposed model can capture the most common features of OC clays.  Figure 11: Comparison between experimental and numerical undrained triaxial compression tests on undisturbed Bangkok stiff clay at Rangsit site. Experimental data from Hassan [13] The first author wishes to thank the Otto Pregl Foundation for financial support. All numerical simulations were performed using the Triax element test driver. The provider Dr. David Mašín is acknowledged.
The hypoplastic flow rule m is with the interpolation factor F given as where the Lode angle θ is defined as: The limit stress condition Y is given as: where I 1 , I 2 , and I 3 are stress invariants expressed as: The multipliers f s , f d , and m R are: with the stress ratio R being defined as: where p + e is the pre-consoldiation pressure in the MCC model, and p * e is the Hvorslev equivalent pressure.
In the above equations, γ, a, and M are material-related constants: Finally, there are six parameters for the proposed model: φ c , λ * , κ * , N , v i , and α.