An improved version of barodesy for clay

Barodesy is a constitutive model based on proportional paths and the asymptotic behaviour of soil. It was originally developed for sand in 2009 by Kolymbas, and a version for clay was introduced in 2012. A shortcoming of former barodetic models was that tensile stresses can occur for certain dilative deformations. In this article, an improved version of barodesy for clay and a simpliﬁed calibration procedure are proposed. Basic features are shown, and simulations of element tests are compared with experimental data of several clay types.


Introduction
The constitutive model for soil called barodesy, proposed by Kolymbas [9][10][11][12], is based on the asymptotic behaviour of granulates expressed by the two rules proposed by Goldscheider [6], which have been experimentally confirmed for sand and clay [3,6,26,27]: (1) starting at the stress-free state, T ¼ 0, proportional strain paths 1 lead to proportional stress paths (see footnote 1); and (2) starting at T 6 ¼ 0, proportional strain paths lead asymptotically to the corresponding proportional stress paths starting at T ¼ 0. This means that proportional stress paths function as attractors.
Barodesy exhibits similarities to hypoplasticity and was introduced by Kolymbas [9] for sand in 2009. In 2012, Medicus et al. [20] modified the sand version [10] and introduced barodesy for clay. A major component of barodesy is the so-called R-function, which links proportional strain paths to proportional stress paths and thus acts as a stress-dilatancy relation. Former versions of R in barodesy [4,[9][10][11][12]20] allow proportional stress paths to reach the tensile area. Experiments by Bergholz [1] of saturated reconstituted clay show that, also for highly overconsolidated clay, the stress ratio q/p (in triaxial compression) does not exceed three (i.e., the stress paths stay in the compression regime). The R-function according to Medicus et al. [21] explicitly prohibits tensile stresses and is chosen as one of the equations for the improved version of barodesy for clay. However, an article [21] reviews existing experimental evidence on stress-dilatancy relations and discusses it in the framework of barodesy, but does not provide a constitutive model.
The main differences of barodesy for clay [20] and the version presented here are shown in this article. The tensor R and the scalar quantities f and g have been changed, and the calibration procedure is simplified as compared to [20].
The calibration procedure is explained, and simulations of element tests are compared with experimental data.

Notation
We use the symbolic notation for Cauchy stress T and stretching D. In some cases, the more familiar symbol r i instead of T i is used for the principal stresses. Stress and stretching are defined as negative for compression. Tensors are written in bold capital letters (e.g., X). jXj :¼ ffiffiffiffiffiffiffiffi ffi trX 2 p is the Euclidean norm of X and tr X is the sum of the diagonal components of X. The superscript 0 marks a normalized tensor, i.e., X 0 ¼ X=jXj. 1 denotes the second-order unit tensor. Stresses are considered as effective ones, and the normally used dash is omitted. T is the objective 2 (co-rotational) stress rate resulting from barodesy, and _ T is the time derivative according to Zaremba/Jaumann, which is obtained by _ T ¼ T À WT þ WT, with W being the antimetric part of the velocity gradient. For rectilinear extensions, the objective stress rate T is equal to _ T and will therefore be used in Sect. 4. The stretching tensor D is the symmetric part of the velocity gradient. 3 The mean effective stress is p :¼ À 1 3 tr T, and e vol ¼ tre is the volumetric strain. In this article, d :¼ tr D 0 ¼ _ e vol =j _ ej is used as a dilatancy measure. For axisymmetric conditions, e.g., conventional triaxial or oedometric compression, the axial stress is denoted with r 1 and the radial stress is denoted with r 2 ð¼ r 3 Þ. The associated strains are e 1 and e 2 ¼ e 3 . The void ratio e is the ratio of the volume of the voids V p to the volume of the solids V s . The deviatoric stress is written as q ¼ Àðr 1 À r 3 Þ, and the deviatoric strain reads e q ¼ 2=3 Á ðe 1 À e 3 Þ. The stress ratio K ¼ r 2 =r 1 at critical states equals K c , and for oedometric normal compression K is denoted by K 0 .

Barodesy for clay
Barodesy is expressed by an evolution equation of the rate type T ¼ hðT; D; eÞ. The general form of the constitutive relation is [9]: Note that c 4 equals 1 for clay, and therefore all material constants are dimensionless 4 . R 0 and T 0 are the normalized tensors of proportional stress paths R and actual stress, respectively. R is a tensorial argument 5 of normalized stretching D 0 and is chosen according to Medicus et al. [21]: Fig. 1 In a, b, w _ e and w r according to [8] are defined in the Rendulic plane: i refers to isotropic compression; o to oedometric compression; c to isochoric triaxial compression; and Àc to isochoric triaxial extension. In c, barodesy for clay [20] (2012) is compared with the version presented here (2016), and u c is chosen as 22:6 2 The term objectivity points to the fact that material behaviour is frame-indifferent, i.e., the behaviour is independent of the observers' motion. 3 In general, stretching D is only approximately equivalent to the strain rate _ e. For rectilinear extensions, D equals _ e, with _ e being the logarithmic strain tensor. 4 However, for sand versions of barodesy (e.g., Kolymbas [10]), the exponent is smaller \1. If the exponent does not equal 1, attention should be paid to the dimensions of material constants, cf. [10]. 5 The exponential of the tensor R can be defined by means of its eigenvalues r i : exp R ¼ exp The scalar functions f and g take into account asymptotic states, critical states, the influence of stress level (barotropy) and density (pyknotropy). They are chosen as follows: with the critical void ratio e c e c ¼ exp N À k Ã ln 2p r Ã À 1 ð8Þ and the scalar functions b and K b ¼ À 1 where r Ã is the reference pressure 1 kPa; c 1c 6 are material constants which depend on the soil parameters u c , N, k Ã and j Ã , cf. Table 1; u c is the critical friction angle; N is the ordinate intercept of the isotropic normal compression line (NCL) in the ln p versus lnð1 þ eÞ plot; k Ã is the slope of the NCL; and j Ã is the slope of the unloading line under isotropic compression in the ln p versus lnð1 þ eÞ plot. A detailed description of the soil parameters is given in Sect. 5. Gudehus and Mašín [8] consider the angles w _ e and w r according to Fig. 1a, b, for a graphical representation of proportional strain and stress paths. In Fig. 1c a w _ e -w r plot of barodesy for clay [20] is compared with the improved version presented in this article. The main difference of the models is the choice of the R-function. Note that the Rfunction according to (3)-(5) prohibits tensile stresses. 6 For volume decreasing paths, i.e., À90 \w _ e \90 the models almost coincide, for certain volume increasing paths 7 . Moreover, the proportional stress paths obtained with the 2012 version [20] lie in the tensile stress region (w r \ À 35:3 or w r [ 54:7 ). The consequences are shown in Fig. 2. Deforming soil with proportional strain Table 1 Determination of the constants c 1c 6 on the basis of u c , N, % À3:6213 Starting from a hydrostatic stress state, a highly overconsolidated Weald clay sample is deformed with proportional strain paths in the range of À180 \w _ e \180 . The stress paths approach the corresponding proportional stress paths. In a, certain stress paths, simulated with barodesy for clay [20], approach proportional stress paths in the tensile region. In b, all stress paths, simulated with the improved version of barodesy for clay according to (1)-(10), stay in the compression regime paths, and the stress paths approach the corresponding proportional stress paths. Certain stress paths approach proportional stress paths in the tensile region with the old version [20], cf. Fig. 2a. In Fig. 2b, all stress paths stay in the compressive area. Note that, for simulations with barodesy for sand [4,[9][10][11][12], qualitatively the same is observed for dense samples. However, this shortage for sand is not as drastic as it is for overconsolidated clay in Fig. 2a.

Calibration
In barodesy, material constants are denoted by c i , unless other symbols are established, such as u c , N, k Ã and j Ã in the case of barodesy for clay. All constants c 1c 6 can be determined on the basis of u c , N, k Ã , j Ã , see Table 1. In order to calibrate the four parameters u c , N, k Ã and j Ã a consolidated undrained triaxial test (CU) is sufficient. From consolidation, we get the parameters N, k Ã and j Ã and from undrained compression the critical friction angle u c is obtained. In Sect. 5 the determination of u c , N, k Ã and j Ã is illustrated by element tests. In Table 2, the parameters are shown for several clay types.
Below, the approach for the determination of c 1c 6 is explained.

Constants c 1 and c 2
The constants c 1 and c 2 can be calculated from u c , cf. Table 1. The R-function (Eqs. [3][4][5] includes c 1 and c 2 , captures critical states and Jáky's relation K 0 ¼ 1 À sin u c for oedometric compression and produces similar results to Chu and Lo's relation [3]. Results are presented in Appendix 1. A detailed explanation and further results are given by Medicus et al. [21].

Constants c 4 and function b
The NCL is used for the determination of c 4 and b. 8 r Ã is a reference pressure equal to 1 kPa. The NCL according to Butterfield [2] as well as the critical state line (CSL) are assumed to be linear in the lnð1 þ eÞ -ln p plot, cf. Mašín [15,18]. The constant c 4 and the function b are chosen in order to ensure that a simulation of hydrostatic normal compression with barodesy starting from e ¼ exp N À 1 yields the NCL. A detailed derivation of c 4 and b is shown in Appendix 2.

Constant c 3
Gudehus and Mašín [8] propose the following graphical representation of admissible states with respect to void ratio and proportional stress paths. Figure 3 shows how proportional stress paths (in terms of w r ) are assumed to be connected with p e =p. Hvorslev's equivalent consolidation pressure p e is the value of mean stress on the NCL, which refers to the current specific volume ð1 þ eÞ, cf. Fig. 3a: The distance of a state characterized by e and p from the isotropic normal compression line is therefore indicated by p e =p. For example, for hydrostatic compression it applies p e =p ¼ 1 and w r ¼ 0, and at critical states p e =p c is assumed to be equal to 2. Proportional stretching will eventually lead to constant values of p e =p for compressive stretching, as well as for extensive stretching, so-called asymptotic extension states [7,8,16]. Asymptotic extension states correspond to so-called normal extension lines in the ln p -lnð1 þ eÞ plot [8,16]. Gudehus and Mašín [8] propose the p e =p-w r plot (Fig. 3b) for the directions of proportional stretching in the range of Àd\w _ e \d, according to Fig. 1a. The directions Àd and d are theoretical limits of asymptotic behaviour according to Gudehus and Mašín [8]. Discrete element simulations by Mašín [16] demonstrated that asymptotic states could only be obtained in a narrower range of w _ e . However, in barodesy, asymptotic states are obtained for the whole range of À180 \w _ e \180 . The following procedure for the determination of c 3 is proposed: barodesy predicts for sufficiently long proportional compressive stretching p e =p ¼ const., e.g., for w _ e ¼ 0 the NCL is reached, and p e =p i ¼ 1 (Fig. 3). This also applies for extension paths, which lead to normal extension lines in the ln p -lnð1 þ eÞ plot. In particular for an isotropic extension path, which is denoted with Ài in Fig. 1a, the isotropic normal extension line is reached, see Fig. 3a. It follows that the unloading stiffness in the ln p -lnð1 þ eÞ plot is characterized by the parameter k Ã : With the general form of the barodetic constitutive relation For a proportional isotropic extension paths, Eq. 15 follows from Eqs. 13 and 14 (with d ¼ ffiffi ffi 3 p ): With f, g and b from (6), (7) and (34) and with ð1 þ eÞ=ð1 þ e c Þ ¼ ð2 Á p Ài =p e Þ k Ã , we obtain: 9 Releasing c 3 , yields: Choosing p Ài =p e ¼ 1=1000 in (16) yields: Note that p Ài =p e ¼ 1=1000 is arbitrary and cannot be acquired by experiments. However, the overall performance 10 of barodesy for clay is best by choosing p Ài =p e ¼ 1=1000 and helps to present a calibration procedure which is simple and applicable also for practitioners without performing any least square optimization. Figure 3b shows how w r is related to p e =p in barodesy.

Constant c 5
The constant c 5 has been determined by trial. Setting c 5 ¼ 1=K c gives the best fit concerning overall performance. Setting c 5 ¼ 1 would highly overestimate radial stress under oedometric compression.

Constant c 6
It appears reasonable to require that, under isotropic extension, the stress paths follow the shortest way to the origin regardless of its actual stress state, cf. Fig. 4a. From this requirement, we get for isotropic extension: and with T 1 =T 2 ¼ T 0 1 =T 0 2 and (1) we obtain: (a) (b) Fig. 3 a Schematic plot: sufficiently long proportional stretching sweeps out the memory. The distance from the isotropic normal compression line is defined through p e =p ¼ const, cf. Gudehus and Mašín [8] and Mašín [16]. b p e =p -w r plot according to [8]: simulations with barodesy for clay 9 At isotropic extension, lnð1 þ eÞ equals N À k Ã lnðp e =p Ài Á p=r Ã Þ, cf. Fig. 3a. We therefore get 1 þ e 1 þ e c ¼ expðN À k Ã lnðp e =p Ài Á p=r Ã ÞÞ expðN À k Ã lnð2 Á p=r Ã ÞÞ 10 The parameter c 3 does not only affect extension states, but also shear stiffness. Equation (20) is valid for a proportional isotropic extension or compression paths or if f ¼ 0. Setting f ¼ 0 in (6), we obtain with d ¼ ffiffi ffi 3 p and b from (34): In Fig. 4b, simulations of barodesy show that the stress paths for isotropic extension follow the shortest way to the origin. Response envelopes for London clay are added in Fig. 4b.

Simulations of element tests
In this section, simulations of element tests with and without rotation of principal axes are shown. Element tests in general are an idealization and, as in all experiments inhomogeneities occur. Especially with shearing, localization takes place and the loss of homogeneity is unavoidable. Thus, the comparison of simulations of element tests with experimental data only serves as an approximate reference.

Rectilinear extensions
Isotropic compression: In Fig. 5, an isotropic compression test and its simulation with barodesy is shown. In Fig. 5b, the calibration of the parameters N, k Ã and j Ã is illustrated with experimental data of Dresden clay [1]. As isotropic normal compression is included in the formulation and calibration of barodesy, normally consolidated isotropic compression test results are therefore in agreement with the simulated NCL, see Fig. 5. The unloading stiffness is described through the parameter j Ã . The term ð1 þ eÞ=ð1 þ e c Þ ¼ const. in (22) indicates a straight line (A-B) in the ln p-lnð1 þ eÞ plot in Fig. 5b, cf. [19] 1 þ e 1 þ e c ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi On this line, the tangential unloading stiffness under isotropic extension is _ p= _ e ¼ Àp=ðj Ã ð1 þ eÞÞ with barodesy. Closer to the NCL, the unloading stiffness is slightly higher, and for lower mean stresses p, the stiffness is lower, cf. Fig. 5b and [19].
Triaxial compression: The critical friction angle u c is calibrated with a normally consolidated Weald clay sample and can be obtained from the slope of the critical state line in the p-q plot, cf. Fig. 6a. Test results and numerical simulation with barodesy of a normally consolidated and overconsolidated sample are shown. Note that the simulation of the overconsolidated sample does not allow a higher mobilized friction angle u m than u c . Barodesy therefore underestimates the peak friction (i.e., the maximum mobilized friction) angle in CU tests, cf. Fig. 6b. The simulations of the normally consolidated and overconsolidated samples in the q/p-e 1 plot coincide for CU tests, cf. Fig. 6b. In Fig. 7, limit points of normally consolidated samples obtained by true triaxial tests are shown. The data refer to San Francisco Bay Mud from Lade [13] and are compared with predictions by barodesy. Note that the critical state locus of barodesy practically coincides with the locus according to Matsuoka-Nakai, cf. Fellin and Ostermann [5].
In Fig. 8 drained triaxial compression and extension tests of normally consolidated and overconsolidated Weald clay are shown. The simulations with barodesy are realistic. Contractant behaviour for the normally consolidated samples and dilatant behaviour for the overconsolidated clay is observed. Peak strength is well predicted. ), the stress paths follow the shortest way to the origin regardless of the actual stress state, cf. schematic plot in a. In b, simulations with barodesy are shown: the stress paths for isotropic extension always follow the shortest way to the origin. Response envelopes for London clay are added In Fig. 9, a more general picture of drained triaxial tests simulated with barodesy is shown. Triaxial tests are shown as p-e and p-q plots as well as plots in the normalized stress plane (i.e., p=p e -q=p e ). The paths approach the critical state line in the p-e and p-q plots. Highly overconsolidated samples dilate to approach the CSL in the pe plot, and slightly overconsolidated and normally consolidated samples exhibit contractant behaviour to approach the CSL. Highly overconsolidated samples overshoot the CSL in the p-q plot.
Oedometric compression: In Fig. 10, oedometric compression of London clay is shown. The normal compression behaviour gives reasonable results in the ep plot, as well as in the r 1 -r 2 plot. 11 For unloading, the radial stress is overestimated.

Rotation of principal stress and strain axes
Simple shear test: Figure 11 presents a simulation of a simple shear test with a constant vertical stress of r y ¼ À100 kPa. The evolution of the shear stress s xy is plotted over the shear strain c (in radian). The angle a r denotes the inclination of major principal stress to the horizontal direction x, and a D is the inclination of major principal stretching, respectively. In Fig. 11, a Weald clay sample with K 0 ¼ 1 À sin u c is sheared. The major (a) (b) Fig. 6 Undrained triaxial compression of Weald clay (according to Mašín [18], data by Parry [22]) and numerical simulation with barodesy. The initial states of the normally consolidated and overconsolidated samples are e 0 ¼ 0:622 and e 0 ¼ 0:572. In a, a p-q plot is shown: the slope of the critical state line in the p-q plot is principal stress direction a r is 90 at zero shear strain and decreases to % 45 with ongoing shear strain. The difference between the angles a D and a r , i.e., the angle of noncoaxiality a D À a r becomes very small 12 , i.e., a r % a D % 45 at the critical state. Similar results with hypoplasticity and an elasto-plastic model are shown in Schranz and Fellin [24]. Experiments on sand according to Roscoe et al. [23] and DEM simulations [25,29] yield similar results, cf. Yu [28]. In Fig. 12, the evolution of the angle of non-coaxiality with ongoing shear strain is shown for different initial K 0 values. In Fig. 12a, DEM simulations from Thornton and Zhang [25], Zhang [29] show that the angle of non-coaxiality is small for K 0 ¼ 1. For K 0 ¼ 2, the angle of noncoaxiality decreases with ongoing shear strain to % 0 ; and for K 0 ¼ 0:5 it increases to % 0 . It is stated that noncoaxiality is significant before 10% shear strain [29]. The predictions with barodesy in Fig. 12b are in good agreement with the DEM simulations in Fig. 12a.
The results of the DEM simulations and experiments [23] apply for sand. Therefore, only a qualitative comparison of barodesy for clay (Figs. 11b, 12b) is possible. However, the comparison demonstrates that barodesy is applicable for general deformation, i.e., rotation of principal stress and strain axes.
Appendix 3 summarizes all equations of barodesy for clay.

Summary and conclusions
Barodesy comprises fundamental characteristics of soil behaviour, such as critical states, asympotic states, barotropy, pyknotropy, and a stress-dilatancy relation. Barodesy can be written symbolically as a single equation of the form T ¼ hðT; D; eÞ, i.e., the stress rate is expressed as a function of the stress, stretching and the void ratio. As in basic hypoplastic models, barodesy uses only T and e as memory parameters. This covers many phenomena and is insufficient to capture strong memory effects. Consequently, ratcheting and unrealistic small-strain behaviour are obtained.
In order to calibrate barodesy for clay, four well-known material parameters of soil mechanics, which can be determined from a consolidated undrained compression test, are sufficient. The model provides realistic results, as compared with experimental results of various clay types.  Mašín [18], data by Parry [22]) and numerical simulation with barodesy. The consolidation pressures are r 2 ¼ À69:45 kPa and r 2 ¼ À206:2 kPa with the initial void ratios eini ¼ 0:572 and eini ¼ 0:622. e 1 versus q plot in a, e 1 versus e vol plot in b Fig. 7 Critical stress points of normally consolidated San Francisco Bay Mud from Lade [13] are compared with critical state predictions by barodesy and Mohr-Coulomb. The calculations refer to tr T ¼ À500 kPa and u c ¼ 30:6 . The samples were isotropically consolidated and compressed in conventional and true triaxial tests. The critical stress points are arranged slightly anisotropic, due to the anisotropic orientation of the particles, according to Lade [13]. An isotropic material would leave a rotation in the deviatoric plane by 120 undiscovered 12 At critical states a D À a r % 0:5 . Neglecting the Zaremba/Jaumann expression ÀWT þ WT yields _ T ¼ T. It follows that a D À a r ¼ 0 at failure.  10 Oedometric loading (up to r 1 ¼ À400 kPa), unloading (up to r 1 ¼ À100 kPa) and reloading (r 1 ¼ À266 kPa): experimental results (PhM14) of London Clay (data from Mašín [14]) with eini ¼ 1:476 and numerical simulation with barodesy. lnð1 þ eÞ versus ln p plot in a, r 2 -r 1 plot in b: r 2 is overestimated at unloading (a) (b) Fig. 11 Simple shear test with a constant vertical stress of r y ¼ À100 kPa, the initial radial stress is r x ¼ ð1 À sin u c Þ Á r y ¼ À59:33 kPa. In b directions of principal stress a r and principal stretching a D are shown. Weald Clay with an initial void ratio eini ¼ 0:68 is simulated with barodesy. Fig. 9 Barodesy: simulated paths of drained triaxial tests for London clay. The start points are denoted by a circle, and the end points are denoted by a cross Fig. 12 Evolution of the angle of non-coaxiality in a simple shear test with different initial K 0 values: in a, DEM simulations from Thornton and Zhang [25], Zhang [29] are shown, in b, Weald Clay with eini ¼ 0:68 is simulated with barodesy  Jáky   Fig. 13 Predictions of proportional paths with the R -function [21] in which c 1 and c 2 are calibrated on the basis u c . The results are in good agreement with the relation by Chu and Lo [3]. u c is chosen to 25 in this plot 13 The equation _ e 1þe ¼ tr D holds for incompressible grains. exp c 5 k Ã ln 2