On the friction dependency of the stress Lode angle

The shape of the failure locus of a material is significant for its strength predictions. Even when constitutive models include the same critical stress surface, different critical stress ratios can be predicted for an identical applied isochoric strain path. In this article, we investigate critical stress predictions of different constitutive models, which include the surface according to Matsuoka–Nakai (MN). We perform analytical investigations, true triaxial test simulations with hypoplasticity and barodesy, and discrete element modelling (DEM) simulations to investigate the friction dependency of the stress Lode angle. Our results demonstrate that in hypoplasticity, the direction of the deviatoric stress state at critical state depends solely on the direction of the applied deviatoric strain path. In contrast, in barodesy, the predictions are also dependent on the friction angle of the material. In addition, we compare these results with those obtained with a standard elastoplastic MN model. To validate this friction dependency on the stress Lode angle, we conduct DEM simulations. The DEM results qualitatively support the predictions of barodesy and suggest that a higher friction results in a higher Lode angle at critical stress state.


INTRODUCTION
3][4][5][6] Among others, the critical stress state in these models is described according to MN.With the critical state depending on the Lode angle, even if the models include the same critical stress surfaces, different critical strength predictions can still be obtained for a given mean effective stress state, especially when the strain path is prescribed.In these cases, the stress response has the freedom to approach the critical state at different Lode angles that is governed by the inner working of the models.In the space of their principal values, the deviatoric direction of a strain path is described by the stretching Lode angle  ε and the deviatoric direction of a stress path is described by the stress Lode angle   .In this work, triaxial compression corresponds to the Lode angles   =  ε = 30 • , triaxial extension to   =  ε = −30 major principal stress ( ′ 2 = ( ′ 1 +  ′ 3 )∕2), the stress Lode angle satisfies   = 0 • .Likewise, for the intermediate principal stretching being exactly between minor and major principal stretching,  ε = 0 • .For an isochoric strain path with  ε = 0 • , the resulting deformation is a plane-strain path with ε1 = − ε3 and with the intermediate principal value ε2 = ( ε1 + ε3 )∕2 = 0.The deviatoric direction of a stress state   for a given strain path  ε = const depends on the stress-dilatancy relation of the constitutive model.2][13] In linear elastic-perfectly plastic models, which use a von Mises plastic potential surface, the stress Lode angle   = 0 • is reached at plane-strain failure. 9,12,14However, for other shapes of plastic potential surfaces, different Lode angles at plane-strain failure can be obtained. 9Previous studies have also shown that the plastic potential can significantly affect the results, as demonstrated by Grammatikopoulou et al. 8 who determined the failure heights of an embankment.Huang et al. 15 performed DEM simulations on sand, and found that the stress Lode angle   at the critical state under plane-strain drained conditions was about 10 • .Experimental studies have shown that the deviatoric direction of the strain increment at failure can differ from the deviatoric direction of the stress state in non-axisymmetric conditions.For instance, Goldscheider 16 observed this effect in sand by applying proportional strain paths with fixed deviatoric strain directions.Yamada and Ishihara 17 similarly reported a small deviation between the deviatoric directions of the strain increment and the stress state for non-axisymmetric conditions.Nakai 18 conducted experiments on sand and clay by applying stress paths with various fixed deviatoric directions while keeping the mean effective stress constant.Based on these tests, Nakai concluded that plane-strain failure occurs in a range of 0 • <   < 15 • , which is similar to the range of 0 • <   < 22 • reported by Potts and Zdravković. 7An isochoric plane-strain deformation is characterised by the stretching Lode angle  ε = 0. Therefore, the stress Lode angle of   = 0 • would indicate that the deviatoric direction of the isochoric plane-strain increment coincides with the deviatoric direction of the critical stress state.Wanatowski and Chu 19 investigated sand under drained and undrained plane-strain conditions and found that the stress Lode angle at planestrain failure was always larger than zero.Specifically, for drained conditions, the stress Lode angle was in the narrow range of   = 13 • -15 • , while for undrained conditions, some values were slightly lower.Furthermore, Selvarajah et al. 20 carried out hollow cylinder tests and observed that undrained plane-strain failure is associated with a range of stress Lode angles of 7 This article investigates the critical state behaviour of the constitutive models, barodesy 4,6 and hypoplasticity 2,5 , under non-axisymmetric strain paths using analytical investigations and true triaxial test simulations.Although barodesy and hypoplasticity include the surface according to MN, they differ in their strength predictions.In addition, we compare these results with those of a standard elastoplastic MN model.Our results demonstrate that in hypoplasticity, the direction of the deviatoric stress state at critical state depends solely on the direction of the applied deviatoric strain path, whereas in barodesy, it depends not only on the deviatoric strain path but also on the friction angle.To investigate the friction dependency of the stress Lode angle, we perform DEM simulations and observe how the internal friction of the material affects the Lode angle at critical stress state.DEM simulations results conform with barodesy where the stress Lode angle at critical state is predicted to increase with friction.

Symbols and notation
Second-order tensors are written in bold capital letters, for example, .Fourth-order tensors in bold calligraphic letters, for example, .We use different kinds of tensor operations using the Einstein summation convention.In particular, the indices follow the lexicographic order ( ∶ )  =     .We denote the unit tensor  of second order with   =   using the Kronecker delta   .‖‖ =

√
tr  2 = √  ∶  is the Frobenius norm of a symmetric tensor , tr  is the sum of the diagonal components of .The superscript 0 indicates a normalised tensor, that is,  0 = ∕‖‖.Any symmetric tensor can be represented in eigenspace by its eigenvalues   .We collect them in the vector [ 1 ,  2 ,  3 ] T .The effective Cauchy stress is denoted by  (the normally used prime is omitted) and the stretching by .Compression is defined as negative.In some cases, the more familiar symbol  ′  instead of   is used for the principal stresses.T is the objective stress rate.Note that with  ′ = − 1 3 tr , the mean effective stress is positive for compression.The stretching tensor  is the symmetric part of the velocity gradient.In general, the stretching  is only approximately equivalent to the strain rate ε.For rectilinear extensions however,  equals ε, with the logarithmic strain .The void ratio  is the ratio of the volume of the voids   to the volume of the solids   .The deviatoric stress  for rectilinear extensions is and the deviatoric strain reads The directions of stress or stretching as vectors in the principal stress or stretching space measured in the deviatoric plane are given by the stress Lode angle   and stretching Lode angle  ε, respectively: with the deviatoric stress with the deviatoric stretching In this work,   is used for MN equivalent friction angle.If a given critical stress state is described by   , the equivalent MN surface includes this particular stress state.

COAXIALITY AND ALIGNMENT AT CRITICAL STATE
In elastoplastic models, the Lode angle at critical state depends on the shape of the plastic potential and on the kinematic constraint of a strain path. 8,9,12,13Hypoplastic and barodetic models do not use the standard notions of elastoplastic models (such as elastic region, yield surface, plastic potential surface and so forth).Therefore, the Lode angle at the critical state must be derived in a different way for these models.For this purpose, we begin with a brief outline of the definitions of coaxiality, alignment and critical state and then carry out related analytical investigations for hypoplasticity and barodesy.

Coaxiality and alignment
Two tensors are called coaxial if they share the same eigenvectors.These eigenvectors can be used as columns of a tensor  with which both tensors can be diagonalised simultaneously.For the tensors  and  used here, this further means that by fixing the first eigenvalue  1 of , the first column of  is defined and with that the first eigenvalue  1 of  is specified.
The axis  1 of the principal stress space and the axis  1 of the principal strain space thus point in the same direction.Alignment 21,22 means that the direction of the deviatoric stress coincides with the direction of the stretching, for which in critical state tr  = 0 and with that  * =  holds.This means that the vector representing the normalised deviatoric stress  0 * and the vector representing the deviatoric stretching  0 * are equal: T , where the minus sign is mechanically meaningless.This also means that the Lode angles of  and  are equal:   =  ε.
Simpler for derivations is the condition which means with a scalar factor  > 0. In mathematics, this is also called colinear.Using Equation (6) for the Lode angle for the stress (3) yields

Critical state
Within the realm of constitutive modelling, the critical state is often envisaged as an asymptotic state towards which the stress  and the void ratio  are approaching for monotonic shearing, either isochoric or with freely evolving volume.The void ratio approaches the critical void ratio   , which is a function of the critical stress state   towards which the stress is approaching.This behaviour can be modelled with asymptotic models where, strictly speaking, the critical state will never be reached in a monotonic shear deformation starting with arbitrary initial stresses and initial void ratio.However, we can start our investigation at critical state, that is, set the stress  =   and the void ratio  =   and perform an isochoric deformation.That means to apply a stretching  =   for which the volume remains constant, that is, tr   = 0.For such a deformation, the stress remains at the critical stress, that is, The evolution of the void ratio  is given by if grain crushing is excluded.The condition tr   = 0 yields ė = 0 and  remains at the critical void ratio.
As tr   = 0, we have

Hypoplasticity
Hypoplasticity is a constitutive framework introduced by Kolymbas. 23The hypoplastic models according to von Wolffersdorff 2 and Mašín 5 apply for sand and clays, respectively.Concepts from critical state soil mechanics are included in the mathematical formulation.Sand hypoplasticity 2 and clay hypoplasticity 5 include the critical stress surface according to MN.The general form of hypoplasticity reads The tensor function  is homogeneous in  of a certain degree and positively homogeneous of first degree with respect to .The fourth-order tensor  and the second-order tensor  depend on the stress and the void ratio.At critical state, holds, which indicates that the critical stress   is an implicit function of the isochoric stretching   =  *  ≠ .The void ratio at critical state is  =   > 0. The critical void ratio is a function of the mean effective stress  ′ , which remains constant tr   at critical state.With that we can write also Since we are studying critical states where the void ratio is a unique function of the stress, we can limit ourselves to functions of  and  without losing generality.Any objective three-dimensional symmetric function (, ) can be represented according to the general representation theorem 24 quoted in Truesdell and Noll [25, eq.(13.7)] The scalar functions   are invariants and joint invariants of  and .Since hypoplasticity can be written as only terms linear in  can occur in Equation ( 14), that is,  4 =  6 =  8 = 0.In the versions of von Wolffersdorff 2 and Mašín, 5 also terms with  2 , , ,  2  and  2 do not appear, that is,  3 =  5 =  7 = 0 and hypoplasticity at critical state (13) reads

Coaxiality at critical state
The stress at critical state can be written in diagonal form with the eigenvalues  , and the associated eigenvectors as columns of the tensor .Multiplication of Equation ( 16) with  T from the left and with  from the right and setting  =  T    yields which can only hold if the off-diagonal elements of  are zero, that is, the diagonal elements are the eigenvalues of   .In other words,   and   can be diagonalised with the same tensor , that is, they are coaxial.

Alignment at critical state
From Equation ( 16), it follows that with the scalar factors  and .The deviatoric part of   is which means alignment, see Equation (6).

Barodesy
8][29] In contrast to hypoplastic models, the basic structure of barodesy has not been derived on the basis of the representation theorem of a tensor function, but is based on the two rules introduced by Goldscheider, 16 which are related to the so-called proportional paths and the asymptotic soil behaviour, proven experimentally for sand 16,30 and clay. 31Proportional stress paths are stress paths for which the ratios of the principal stresses remain constant:  ′ 1 ∕ ′ 2 = const,  ′ 2 ∕ ′ 3 = const and thus ∕ ′ = const.In the same way, paths with constant ratios of the principal strains are called proportional strain paths.
As in hypoplasticity, the objective effective stress rate is formulated as a function of effective stress, stretching and void ratio.Barodesy also comprises concepts from critical state soil mechanics, like a stress-dependent critical void ratio   .Unlike hypoplasticity, analytic splitting of barodesy into parts that are linear and nonlinear in  is not possible. 32he core of barodesy is the so-called  function, 4 which maps the directions of the stretching to directions of the stress state.There a several versions for that function, however, they have a common formulation at critical state The constant  can be computed from the critical state friction angle   : with   = 1 − sin   1 + sin   .
In critical state, holds, which means with  < 0. Equation ( 24) defines a cone in the principal stress space, the so-called critical stress-states, which practically coincides with the MN cone, as shown by Fellin and Ostermann. 33

Coaxiality at critical state
The stretching at critical state can be written as with the eigenvalues  , and  consisting of the eigenvectors.The stress at critical state follows from Equation ( 24) with the definition of the matrix exponential which implies that   and   share the same eigenvectors and are, therefore, coaxial.

Alignment at critical state
Since   and   are coaxial, we can consider eigenvalues of these tensors in the following derivation without loss of generality, and Equation ( 24) reads: The deviatoric stress at critical state is Alignment (6) requires Here we used that tr   = 0 and therefore  0 * , =  0 , .With  = Without loss of generality, we consider the following ordering of the eigenvalues  ,1 ≤  ,2 ≤  ,3 < 0 and compute from Equation ( 30) and which can be interpreted as two negative slopes of secants of the logarithm function through (− ,1 , ln(− ,1 )).Since they are equal, this implies  ,2 =  ,3 , that is, axisymmetric triaxial conditions.With tr   =  ,1 + 2 ,2 it follows that This yields For alignment, the same must hold for the stretching, that is, and This is true for axisymmetric triaxial compression and for triaxial extension, with another ordering of the eigenvalues.
Barodesy shows alignment at critical state only for axisymmetric triaxial compression and extension.Alignment is not fulfilled for all other deformations.

CRITICAL STATES IN THE DEVIATORIC PLANE -TRUE TRIAXIAL TESTS WITH HYPOPLASTICITY AND BARODESY
Numerical simulations of true triaxial tests enable us to investigate different deviatoric directions of stress and strain paths.In this section, we investigate the behaviour of hypoplasticity and barodesy under different, fixed deviatoric directions  ε.Specifically, we analyse the models' predictions for the stress Lode angle at critical state.Additionally, we compare our findings to those obtained with an elastoplastic MN model 9 , as the failure criterion is then consistent with hypoplasticity and barodesy.Other studies on elastoplastic models investigate in detail the influence of the plane-strain constraint and different shapes of the plastic potential and yield surfaces. 8,9,12,13

Simulations with hypoplasticity
We carry out true triaxial tests with constant mean effective stress as well as isochoric paths with clay hypoplasticity 5 and sand hypoplasticity. 2The used material parameters are given in Tables 1 and 2.

Paths with constant mean effective stress and fixed deviatoric directions of strain paths
The simulations presented in Figure 1 employ clay hypoplasticity 5 and investigate the deviatoric directions of stress paths for a given deviatoric direction strain path.The mean effective stress is constant at  ′ = 150 kPa during all tests, and the initial state is highly overconsolidated, with an initial overconsolidation ratio (OCR) of 4. The deviatoric directions  ε are are set to  ε ∈ {−30, −20, … , 30} • in the individual tests.Figure 1A displays the stress ratio ∕ ′ versus deviatoric strain   response for the different deviatoric directions, which reveals variations in the critical stress ratio due to the non-circular shape of MN.The MN surface in terms of principal stresses is described by see also Refs. 1, 36.At critical state, the MN parameter  MN equals  MN,c , and thus the critical stress surface is defined: Figure 1B illustrates the evolution of the MN parameter (36) with increasing deviatoric strain.It is noteworthy that the peak states nearly lie on a MN equivalent surface.At critical state,  MN approaches by definition  MN,c for any deviatoric direction, as shown in Figure 1D, which displays the intersection of the deviatoric plane with the critical stress surface according to MN for  ′ = 150 kPa. Figure 1C shows the evolution of the Lode angle   with ongoing deviatoric strain.Throughout the simulation with hypoplasticity,   remains equal to  ε.Although an anisotropic stress state would initially affect the evolution of the Lode angle   for a given  ε, at critical state, it holds   =  ε in hypoplasticity, irrespective of the initial stress state and initial density.

Isochoric paths with fixed deviatoric directions of strain paths
In Figure 2, simulations of true triaxial tests with tr  = 0 using sand hypoplasticity 2 are presented.The initial conditions are described by a hydrostatic stress state  ′ ini = 150 kPa and by an initial void ratio 0.85, indicating an initially loose state.The deviatoric directions  ε are again kept constant at  ε ∈ {−30, −20, … , 30} • .Figure 2A displays the stress paths.Figure 2B shows the evolution of the MN parameter  MN .As for the clay model,   =  ε at critical state, see Figure 2C,D.

Simulations with barodesy
In this section, we revisit the simulations conducted in Section 3.1 using barodesy for clay 6 and barodesy for sand 4 .The material parameters are presented in Tables 1 and 3.The resulting predictions are similar to those obtained with hypoplasticity, as depicted in Figures 1 and 2. However, one notable difference between the two models is observed in the Lode angle at critical stress state, which exhibits some deviation in the case of barodesy.

Paths with constant mean effective stress and fixed deviatoric directions of strain paths
Figure 3 shows simulations with barodesy for clay, where different deviatoric directions  ε are applied, ranging from −30 • to 30 • .Again the mean effective stress is kept constant throughout the simulations at  ′ = 150 kPa. Figure 3A shows the stress ratio ∕ ′ versus deviatoric strain   response.Similar to the hypoplastic predictions, the difference in the critical stress ratio for different deviatoric directions is due to the non-circular shape of the critical stress surface following MN.In Figure 3B, the evolution of the MN parameter with deviatoric strain is shown.The  MN versus   responses practically overlap, indicating that peak states lie on a MN equivalent surface if the initial stress state is hydrostatic.At critical state,  MN coincides by definition with  MN,c for any deviatoric direction, as seen in Figure 3D. Figure 3C displays the evolution of the Lode angle   with ongoing deviatoric strain.In contrast to hypoplasticity,   differs from  ε at critical state in the simulations with barodesy.In all simulations, the Lode angle   exceeds  ε at critical state.

Isochoric paths with fixed deviatoric directions of strain paths
In Figure 4, simulations of true triaxial tests on initially loose sand samples with tr  = 0 are presented using barodesy for sand. 4As with the previous simulations, the deviatoric directions  ε are kept constant at  ε ∈ {−30, −20, … , 30} • .The stress paths are shown in Figure 4A, while the MN parameter  MN evolution with deviatoric strain is shown in Figure 4B.The stress-strain response of all tests virtually overlap, indicating that for an initially isotropically consolidated state, the predictions are on a MN equivalent surface for a given   = const.As with barodesy for clay, in all simulations with barodesy for sand, the stress Lode angle   exceeds the stretching Lode angle  ε at critical sate, as shown in Figure 4C.Once again,   at critical state is independent of the initial conditions, that is, initial void ratio and initial stress state.However, different values of   are approached for barodesy for sand and clay.This is because, in barodesy, the Lode angle   for a given  ε is friction-dependent and, therefore, depends on the parameter   .

Friction dependency of the stress Lode angle in barodesy
In Figure 5, the critical stress states in the deviatoric plane are depicted: the full line represents the cross section of the critical stress surface according to MN with a deviatoric plane.The arrows indicate the directions of isochoric strain In barodesy (solid lines), the results are friction-dependent:   increases with increasing   for a constant value of  ε .The results marked with  ε = 0 • ( * ) refer to plane-strain isochoric predictions.The dashed line in (A) refers to hypoplastic predictions.
increments and are characterised through  ε.The directions of proportional stress paths in the deviatoric plane are described by the Lode angle   .In hypoplasticity,   =  ε for isochoric conditions, see Figure 5A.However, in barodesy,   differs from  ε, as illustrated in Figure 5B,C.Figure 5C shows the critical stress states according MN for different values of   .For a given  ε, (e.g. ε = 0 marked with * ), the Lode angles   differ for different values of   .Figure 6 presents an analysis of the friction dependence of stress Lode angles in barodesy.Figure 6A shows a  ε versus   plot, which demonstrates that for a given value of  ε,   increases with increasing friction angle   .This trend is further illustrated in Figure 6A,B.Figure 6C compares the Mohr-Coulomb equivalent friction angle  MC with the MN equivalent friction angle   .If a given critical stress state is described by  MC , the equivalent Mohr-Coulomb surface includes this particular stress state.If a given critical stress state is described by   , the equivalent MN surface includes this particular stress state.It is worth noting that the difference between MN and Mohr-Coulomb is most significant (as indicated by the maximum values in Figure 6C, for plane-strain critical states predicted with barodesy ( ε = 0, marked with * ).

Elastoplastic Matsuoka-Nakai model
In addition to the predictions made using hypoplasticity and barodesy, we also compare our results to those obtained using a standard elastoplastic model.Lagioia and Panteghini 9 investigated an elastoplastic model that employs MN for both the yield as well as the plastic potential surfaces.By using a non-associative flow rule for this model, the stress Lode angle at plane-strain failure can be determined from the angle of dilatancy using Equation (38), as shown in Figure 7: We simplified Equation ( 38) from an equation originally proposed in Lagioia and Panteghini (in their work labelled as Equation ( 37)). 9 For critical state, the value  = 0 • in Equation (38), leads to a Lode angle   of zero degrees.For  = 0 • , the MN plastic potential surface becomes cylindrical (i.e. the von Mises criterion), and therefore the isochoric-plane strain path is perpendicular to the von Mises plastic potential surface only at   = 0 • , see also Refs.9, 12.This result is consistent with the predictions obtained from hypoplasticity under plane-strain, isochoric conditions.A simple linearisation of Equation ( 38) is possible, since for  = 90 • , we get   = 30 • : This is a very good approximation, with accuracy better than 0.25 • , as shown in Figure 7.

(A) (B)
F I G U R E 7 (A) Stress Lode angle for an elastoplastic Matsuoka-Nakai model at plane-strain failure in dependence of the angle of dilatancy  according to Lagioia and Panteghini. 9(B) The difference between Equation ( 38) and its approximation (39) is smaller 0.25 • for 0 ≤  ≤ 90 • .The bullet in (A) marks the here considered critical state situation with  = 0.

DEM SIMULATIONS VERSUS CONSTITUTIVE MODELS' PREDICTIONS
To the authors' knowledge, previous studies have not yet investigated the friction dependency on the stress Lode angle.
To validate the predictions of the constitutive models and to enhance the understanding of geomechanics in this respect, discrete element modelling (DEM) simulations are suitable for qualitative analysis.In this section, we simulate isochoric, plane-strain compression tests using DEM and compare the results with the predictions made by barodesy and hypoplasticity.The plane-strain path is chosen as representative strain path for non-axisymmetric conditions because it is between axisymmetric compression and extension.Moreover, boundary value problems are often idealised under plane-strain conditions, making this strain path particularly relevant in geotechnical engineering.During isochoric (tr  = 0) plane-strain compression, the value of  ε is zero.In the simulations, friction is varied to investigate its effect on the Lode angle   .

DEM simulations
A series of 3D DEM simulations have been carried out on assemblies of 32,000 polydisperesed, cohesionless spherical particles and considering a range inter-particle friction values.The radii of the particles are uniformly distributed between [0.67, 1.23] mm.Periodic boundaries were adopted to prevent heterogeneities associated with shear bands formation.
During sample preparation and hydrostatic compression stage, the inter-particle friction coefficient,   , was set to an artificially low value of 0.08 to achieve a relatively dense sample.The sample was first consolidated to a hydrostatic stress of 300 kPa.After reaching equilibrium at the end of the consolidation stage, the inter-particle friction was increased to its final value and the deviatoric loading was conducted under strain-controlled plane-strain isochoric conditions.The contact stiffness along the normal and tangential directions,   and   , and the rolling stiffness,   , are calculated as follows: where   is a stiffness parameter, and  1 and  2 are the radii of the particles in contact.In this study, we adopt the calibrated values in Sibille et al. 37 who represented Hostun sand response with the following values:   = 500 MPa,   = 0.3 and   = 5.0.The rolling moment and the tangential forces at contact are limited by rolling resistance and inter-particle friction, respectively.The limit for plastic rolling to occur is given by     min{ 1 ,  2 } where   is the magnitude of the normal force at the contact and   is the rolling resistance which is kept at 0.55.The inter-particle friction coefficient,   , limiting the tangential force is varied systematically in the range   ∈ [0.13, 0.25] to observe its effect on stresses at the critical state.The straining rate was kept low enough to ensure a quasi-static condition with the ratio of average unbalanced force to the average contact force kept below 10 −3 .All deviatoric tests started from an identical consolidated sample.The results of the DEM simulations are presented in Figure 8.The samples are subjected to isochoric, plane-strain paths, while the inter-particle friction is varied for different simulations.A change in friction affects the shape of the stress paths as shown in Figure 8A.As the friction is increased, the initial contractive response (manifested as a decrease in  ′ ) and the ensuing phase transformation disappears to be replaced by a mostly dilative response where  ′ increases throughout the deviatoric loading.
Increasing the friction also leads to higher stress ratios ∕ ′ and MN parameter values  MN , as illustrated in Figure 8B

Barodesy and hypoplasticity
In both models, barodesy and hypoplasticity, the critical stress predictions are controlled by the critical friction angle   .The influence of varying   on the critical strength predictions can be observed in Figures 9 and 10.An increase in the friction angle leads to an increase in the stress ratio ∕ ′ as well as the MN parameter  MN , as shown in Figures 9A,B and 10A,B.Figures 9C and 10C display the stress Lode angle under critical state conditions for plane-strain deformation.In the case of barodesy (Figure 9C), an increase in   leads to an increase in the stress Lode angle, which is corroborated by the results of the DEM simulations in Figure 8.The DEM results are of qualitative value and are appropriate for tendencies rather than precise quantitative predictions.Nevertheless, the simulations clearly show that the Lode angle at the critical stress state increases with friction, which is consistent with the predictions of barodesy.Conversely, hypoplasticity shows no dependence of the Lode angle on   , as  ε = 0 implies that   = 0, as shown in Figure 10C.To further illustrate the relationship between the stress Lode angle and friction angle, Figure 11

Discussion
The investigated constitutive models, barodesy and hypoplasticity, share similarities in their mathematical formulation and concepts, often leading to comparable predictions.The question arises whether the differences in the predicted stress Lode angle between the investigated constitutive models, barodesy and hypoplasticity, are significant.Thus, it is of interest to quantify this difference, particularly for plane-strain predictions.
In Figure 12, the failure surface according to MN is depicted for different values of   .Figure 12B displays the stress ratio ∕ ′ for Lode angles in the range of −30 • ≤   ≤ 30 • .For triaxial extension (  = −30 • ), ∕ ′ is 6 sin   ∕(3 + sin   ), for triaxial compression (  = 30 • ), ∕ ′ is 6 sin   ∕(3 − sin   ), for   = 0 • , the stress ratio is 38 : At high values of   , a change in the Lode angle   has a more pronounced effect on variations in the stress ratio and critical strength predictions according to MN.For a critical friction angle of   = 20 • , for instance, the difference in the predicted stress Lode angle between barodesy and hypoplasticity is Δ  = 6.7 • , with barodesy approaching   = 6.7 • and hypoplasticity predicting   = 0 • under plane-strain conditions.However, the resulting variation in stress ratio Δ(∕ ′ ) due to the change in Lode angle is relatively small at Δ(∕ ′ ) = 0.027, representing only a 4% increase, and can be considered negligible, as shown in Figure 12.However, for larger values of   , the impact of variation in stress Lode angle becomes more significant.For example, for   = 30 • , the difference in Lode angle Δ  = 10 • under plane-strain critical states with barodesy leads to a 9% increase of Δ(∕ ′ ) = 0.086 in stress ratio compared to hypoplasticity.Similarly, for   = 40 • , the difference in Lode angle Δ  = 13.4 • leads to a 16% increase (Δ(∕ ′ ) = 0.19) in predictions with barodesy compared to hypoplasticity.For plane-strain critical stress states, the stress ratio predictions of hypoplasticity are always conservative compared to barodesy.Our investigations have revealed that a difference in the predicted Lode angle   under plane-strain conditions using MN can have a significant impact on strength predictions, particularly for larger values of   .

SUMMARY AND CONCLUSIONS
This study examined the critical state behaviour of the constitutive models hypoplasticity and barodesy, as well as DEM simulations, under isochoric, proportional strain paths.In both constitutive models, the stress Lode angle at critical state is independent of the initial stress state and density, as critical states are asymptotic states.Our study demonstrated, both analytically and numerically, that the Lode angles of stress and stretching coincide for hypoplasticity at critical state, which means that the normalised deviatoric stress and stretching are aligned.In contrast, the Lode angle   in barodesy depends on friction and always exceeds  ε at critical state.This is consistent with the DEM simulations presented in the article, which show that the stress Lode angle at critical state for an applied isochoric plain-strain path increases with friction.The insights gained from this study will help deepen our understanding of the critical state behaviour of soils and will have an influence on the development of constitutive models.

A U T H O R C O N T R I B U T I O N S
G.M. led the conceptualisation of the study, with input from all authors.G.M. performed the simulations with the constitutive models, and wrote Sections 1, 3-5 (except Section 4.1).M.P. performed the DEM simulations and wrote Section 4.1.A.O. formulated equations (38) and (39).A.O. and W.F. performed the analytical calculations and wrote Sections 2.1 and 2.2.All authors discussed the results and edited the whole manuscript.

1
True triaxial test simulations with clay hypoplasticity: the initial OCR ini equals 4, and the mean effective stress  ′ = 150 kPa is kept constant in the simulations.Various paths with  ε = const are applied.

4 5
True triaxial test simulations with barodesy for sand: the initial void ration  ini equals 0.85, the initial mean effective stress  ′ ini = 150 kPa.Isochoric strain paths with various deviatoric directions  ε = const are applied.Critical states in the deviatoric plane: the full line is the cross section of the critical stress surface according to Matsuoka-Nakai with a deviatoric plane, the arrows denote the deviatoric directions of proportional strain paths and are characterised through  ε .For isochoric (tr  = 0) plane-strain compression, it follows  ε = 0 • .(A) In hypoplasticity,   =  ε for isochoric conditions.(B), (C) In barodesy,   differs from  ε .(A) and (B) are modified from Medicus et al.14

F I G U R E 8 F I G U R E 9
Results of DEM simulations of isochoric, plane-strain paths (i.e.tr  = 0,  ε = 0 • ) with varying friction.The Matsuoka-Nakai parameter  MN is in the range of 10.2-11.3 at critical state.The related Matsuoka-Nakai equivalent critical friction angles are 21.0 • -28.4 • .The results show that the stress Lode angle depends on friction.Simulations with barodesy for sand of isochoric, plane-strain paths with varying critical friction angle (i.e.tr  = 0,  ε = 0 • ) with varying friction   ∈ {15, 20, … , 45} • .The results show a friction dependency of the stress Lode angle.
,C.The critical state is characterised by  MN values in the range of 10.2-11.3,corresponding to MN equivalent critical friction angles between 21.0 • and 28.4 • .In Figure 8D, the evolution of the stress Lode angle with ongoing deviatoric strain is presented.It is observed that an increase in friction leads to a corresponding increase in the Lode angle at critical stress state.Specifically, when deformation is restricted in the plane-strain direction ( 2 = 0), materials with higher friction develop a lower value of the normalised intermediate principal stress | ′ 2 |∕ ′ , which corresponds to higher values of Lode angle,   .The commonly used relation suggested by Jáky suggests a relation between the coefficient of lateral earth pressure at rest,  0 , and the friction angle at the critical state: An increase in   leads to a decrease in the radial stresses | ′ 2 | = | ′ 3 | (with  2 =  3 = 0) during oedometric compression.Similarly, under isochoric plane-strain deformation, an increase in   results in a reduction of the intermediate principal stress | ′ 2 |, which corresponds to the direction with  2 = 0.

F I G U R E 1 0F I G U R E 1 1
Simulations with sand hypoplasticity of isochoric, plane-strain paths (i.e.tr  = 0,  ε = 0 • ) with varying critical friction angle   ∈ {15, 20, … , 45} • .The stress Lode angle is independent on friction.Comparison of the stress Lode angle predictions for hypoplasticity and barodesy under plane-strain conditions with varying friction.The results show that an increase in friction leads to an increase in the Lode angle for barodesy and the DEM simulations, while hypoplasticity and the elastoplastic MN model for  = 0 predict a zero Lode angle for any plane-strain critical stress state.
compares the Lode angle with increasing MN equivalent friction angle   for barodesy, hypoplasticity, the elastoplastic MN model for  = 0 • and the DEM results.The DEM results support the friction dependency of the stress Lode angle.The predictions of hypoplasticity coincide with the elastoplastic model presented in Section 3.4.

2
Effect of the stress Lode angle under plane-strain conditions on strength predictions.The plane-strain critical stress states are marked with * for barodesy and with • for hypoplasticity.The larger   is, the higher is the difference between the two models.The filled, red squares refer to the DEM results.
Parameter set for the clay models.
F I G U R E 3True triaxial test simulations with barodesy for clay: the initial overconsolidation ratio OCR ini equals 4, and the mean effective stress  ′ = 150 kPa is kept constant in the simulations, various paths with  ε = const are applied.