On stress ratio equations in three-dimensional stress space for modelling soil behaviour

a

centre of circle in proposed stress ratio function

b

alternative third stress invariant

CAX4

Four-node axisymmetric element

f

yield surface

G

shear modulus

g

tangent line to yield surface in π-plane

H₀

NorSand constant plastic hardening modulus

H_ψ

NorSand state-based plastic hardening modulus

k

ratio of triaxial extension to compression stress ratios

M

stress ratio

M_tc

stress ratio at the triaxial compression meridian

M_te

stress ratio at the triaxial extension meridian

N

volumetric coupling coefficient

p

mean stress invariant

q

deviatoric stress invariant

R

radius of circular arc in proposed stress ratio function

STOL

Stress tolerance

TXC

Triaxial compression

TXE

Triaxial extension

UMAT

User material

x, y

Cartesian coordinates

α

parameter for the Van Eekelen (1980) stress ratio function

Γ

intercept of semi-log critical state line at 1 kPa

θ

Lode angle

κ₁

parameter for the Lade-Duncan (1975) stress ratio function

λ

slope of semi-log critical state line

ν

Poisson’s ratio

σ

second-order Cauchy stress tensor

σ

normal stress acting on a plane

σ₁

major principal stress

σ₂

intermediate principal stress

σ₃

minor principal stress

τ

shear stress acting on a plane

χ_tc

state-dilatancy parameter defined at triaxial compression

ϕ

friction angle

ψ

soil state parameter

Constitutive models are a set of equations that connect a material’s elastic and plastic strains to its state of stress and the changes of its state. These models are commonly written in terms of stress and strain invariants, which are scalar proxies for the stress and strain tensors, or more familiarly their principal values. There are various sets of invariants that can be used, but the most common ones represent a hydrostatic component, a shear or deviatoric component, and a third invariant that represents the relative magnitude of principal strains or stresses. It is customary to develop the constitutive equations with the first two invariants under triaxial conditions where the third invariant is constant. Once the equations are developed, various model parameters are tied to the third invariant, making the model capable of operating in the three-dimensional stress space.

In geomechanics for example, constitutive equations are written on p, the mean stress, and q, the deviatoric stress, under triaxial compression conditions given its ease of understanding, abundance of data, and fully known states of stress and strain. A third invariant, the Lode angle θ, is then introduced to extend the equations to the three-dimensional stress space by incorporating the effect of intermediate stress on behaviour. When the mobilized, peak, or critical state stress ratios become functions of the Lode angle, the yield surface and its plastic strain twin, the plastic potential surface, assume a spatial shape such as a cylinder, a cone, a bullet shape, or an ellipsoid.

Tests on soil in three-dimensional apparatuses have repeatedly shown that soil behaviour changes as the intermediate principal stress varies from being equal to the minor principal stress (triaxial compression) to the major principal stress (triaxial extension). This observation has been repeated using cubical (Sutherland and Mesdary, 1969; Lade and Duncan, 1973), plane strain (Cornforth, 1961; Cornforth, 1964; Wanatowski and Chu, 2007), and hollow cylinder specimens (Yang, 2013).

Isolating the influence of the intermediate stress from that of dilation and other elements of soil behaviour limits the amount of usable data to critical state conditions, where other effects have been exhausted. Published test data that reliably reached the critical state are largely concentrated around strength conditions near triaxial compression on a small selection of sands (Jefferies and Been, 2015).

The most common failure criterion used in geomechanics is probably Mohr-Coulomb (1776). It is also popular as an elastic-perfectly plastic model when the failure surface takes the form of a yield surface. In the shear stress-normal stress space the yield/failure envelope is linear with a slope represented by the friction angle, ϕ. This friction angle is a strength parameter which can also be represented by a stress ratio at failure conditions. The model sometimes includes apparent cohesion to represent the increase in dilatancy due to over-consolidation at lower confining stresses, which is often taken to be zero in the case of granular soils. The shear stress-normal stress space is three-dimensional in nature and the Mohr-Coulomb (1776) failure criterion produces one of the many existing functions for extending all constitutive models, simple and sophisticated, to three-dimensions.

This paper examines of some of the popular equations that are used in geomechanics to generalize constitutive models to any stress state and discusses what properties such equations should ideally exhibit. A new numerically stable and experimentally compatible function is proposed.

The stresses at a point in a material can be represented via the second order Cauchy stress tensor σ, with positive compressive stresses being the sign convention used in this work, as is typical in geomechanics. This symmetric tensor has six independent components which represent the full state of stress on any three orthogonal planes with both normal and shear components. The eigenvalues and eigenvectors of this tensor represent the principal stresses and their respective directions. The principal stresses are invariants of the stress tensor in the sense that they retain the same magnitude and direction regardless of the orientation of any infinitesimal element at the material point. Sets of scalar invariants defined using the magnitudes of principal stresses are often the building blocks of constitutive models. p, q, and θ as defined by Equations 1 through Equations 3 are popular in geomechanics, having the benefit of a comprehensible physical meaning and the fact that they can be visualized in the principal stress space easily. The mean stress, p, is the average of the three principal stresses and represents the hydrostatic component of the stress tensor. The deviatoric stress, q, is a function of the differences in principal stresses, and represents the magnitude of the stress causing distortion. The Lode angle, θ, is the third stress invariant and represents the relative magnitudes of principal stresses.

If the three principal stresses are plotted orthogonally to form the stress space shown in Figure 1(a), the stress invariants all have unique visualizations. A bisector of the space where the three coordinates are equal is the hydrostatic axis. This axis is also normal to the π-plane (octahedral or deviatoric plane) and intersects the π-plane at the origin of its polar coordinates. p is proportional to the distance of a given π-plane to the origin of the stress space, q is proportional to the polar coordinate radius of the π-plane to the yield surface, and θ is the angular coordinate measured from the angle bisector as illustrated in Figure 1(b).

The Lode angle, θ is by definition restricted to the interval of [−π/6, π/6] radians because it is one third of an arcsine as shown in Equation 3. A value of θ = π/6 represents triaxial compression (TXC) and triaxial extension (TXE) is represented by θ = −π/6. The reference sector with the θ axis shown in Figure 1(b), is where σ₁ ≥ σ₂ ≥ σ₃ and σ₁, σ₂, and σ₃ are the major, intermediate, and minor principal stresses, respectively. The other five sectors represent other orders of σ₁, σ₂, and σ₃.

The formulation of elasto-plastic models has a few important components that are almost ubiquitously shared among all models, albeit with different levels of complexity. The most common component is the yield surface, which is the locus of stress states satisfying the equation f = 0 representing the threshold between elastic and elasto-plastic behaviour. This surface is defined in the reference sector, and then reflected around the Lode angle limits for the other five sectors in stress space.

Other important components include the plastic potential surface which dictates the direction of plastic strains (flow rule) relative to the stress state, a hardening rule that controls how the size and shape of both yield and plastic potential surfaces evolve, and the consistency condition. The consistency condition closes the set of constitutive equations so that the model is not under-defined by taking the total differential of the yield surface and setting it to zero.

Advanced constitutive models are often defined in the p–q space with a stress ratio M linked to various components of the model. This stress ratio is a constant for triaxial stress states and ratio of deviatoric to mean stress q/p becomes equal to M when the sample has failed. The stress ratio itself can take various forms, including the critical state friction ratio, mobilized friction ratio, dilatancy line, bounding line, and so on in various models (Ishihara, Tatsuoka and Yasuda, 1975; Dafalias and Manzari, 2004; Jefferies and Been, 2015). To make models three-dimensional and represent the impact of intermediate principal stress on soil strength, the stress ratio is made into a functional form of the stress state, turning the various (ratio) lines into surfaces. This functional form of the stress ratio is often user-defined in constitutive modelling, as the function can be multivariable and made up of any combination of the stresses acting on an element as well as material parameters (Brinkgreve, 1994; Halabian et al., 2022). Typically for soils, stress ratio functions assume the form of a single-variable function of the Lode angle, M(θ). In the rest of this work, all discussions and formulation are equally applicable to all plasticity surfaces that are tied to the stress ratio. Experimental data at critical state are used to assess various equations, given their availability and ubiquity among models. The term “yield surface” will be used in a general sense to represent all such surfaces (e.g. yield, plastic potential, bounding, dilatancy, etc.). The functional form of the stress ratio M(θ) will be used to represent the shape of these surfaces (including the yield surface) on the π-plane. The values of M(θ) at triaxial extension and triaxial compression meridians are specifically denoted as M_te and M_tc, and their ratio (M_te/M_tc) is represented by the variable k.

A well-behaved constitutive model needs a yield surface that is differentiable (smooth) and convex. A continuously differentiable yield surface is necessary for the flow rule to be defined everywhere and for the consistency condition to be computable. Assuming that the yield surface is differentiable means that the partial differential of the generalized yield surface on the π-plane with respect to the Lode angle, ∂M(θ)/∂θ, needs to be continuous. The continuity is not satisfied where the yield surface reflects in polar space unless ∂M(θ)/∂θ = 0 at those reflection meridians.

Drucker’s (1959) postulate, which establishes that loading a material followed by unloading to the original stress level dissipates energy, is a manifestation of the second law of thermodynamics in order to avoid the creation of energy. It leads to the postulate of maximum plastic work, which requires the yield surface to be convex. The mathematical condition for a yield surface function M(θ) to maintain convexity at a given point on the π-plane is derived in  Appendix A.

The smoothness and convexity of the yield surface, which are controlled by the stress ratio function, are most easily viewed on the π-plane. As shown in Figure 2, the cross-section of yield surfaces with the π-plane takes on various shapes including the familiar irregular hexagon of Mohr-Coulomb (1776) or the circle of Drucker-Prager (1952).

A 3-D generalization for the stress ratio would ideally satisfy the following conditions:

1. Match available experimental data, especially under triaxial conditions M_te and M_tc;

2. Be simple to numerically implement;

3. Have a continuous derivative; and

4. Remain convex throughout the Lode angle domain.

Some of the existing generalizations of M(θ) are summarized in Table 1. The equations are presented in a consistent normalized form. The values of ∂M(θ)/∂θ at the limits are summarized, with zero representing a continuous derivative. Their convexity conditions were also computed using the equations described in  Appendix A. Figure 3 compares these equations, normalized to M_tc to make all the entries comparable, with experimental data and their strengths and shortcomings are discussed below. Where equations allow it, k = M_te/M_tc is set to 0.75 as an intermediate value based on the experimental data. For functions that only have the parameter M_tc in their formulation, a value of 1.3 is used based on the experimental data. Laboratory testing for truly three-dimensional stress conditions in soils where the complete stress state is known is rare, and rarer still is this kind of test data run monotonically to high strain where the soil reaches a clearly identifiable critical state. Moreover, nearly all of the available test data that is suitable is for positive values of Lode angle, with much of it being between the values of 15° and 30° (Jefferies and Been, 2015). Plane strain tests by Cornforth (1961, 1964) on Brasted sand (M_tc = 1.31) and by Wanatowski and Chu (2007) on Changi sand (M_tc = 1.35), as well as hollow cylinder tests by Yang (2013) on Leighton Buzzard sand (M_tc = 1.27) are plotted in Figure 3 along with the stress ratio plots. Importantly, the data are those interpreted to be at the critical state to make sure dilation effects are not being mixed into the influence of θ on M. The limits of the domain at ±π/6 are the triaxial compression and extension meridians, M_tc and M_te, respectively, where the function reflects on the π-plane to form the adjacent sectors. Sharp corners on the π-plane, such as those seen in the Mohr-Coulomb (1776) envelope in Figure 2, will appear as non-flat ends in Figure 3. For an equation of M to be smooth at the transition between sectors on π-plane, its derivatives with regards to θ, ∂M(θ)/∂θ, must be zero at θ = ±π/6.

The basis for many of the stress ratio formulations proposed for the cross-section of the yield surface (or stress ratio in general) on the π-plane is the Mohr-Coulomb (1776) failure criterion. The criterion itself is a simple manifestation of the fundamental mechanics of friction that, for a given plane, linearly relate the normal stress, σ, and shear stress, τ, through the coefficient of friction, tan(ϕ), where ϕ is the friction angle. Mohr-Coulomb’s assumption of a constant friction angle produces a value of k for a given value of ϕ, making the shape of the hexagon in Figure 2(a) function of the friction angle.

Friction angle ϕ and stress ratio M are directly related through the parameter b, which is an alternative third stress invariant. Equation 4(a) shows this relation, and Equation 4(b) shows the definition of b as a ratio of principal stress differences. For example, for ϕ = 30°, at triaxial compression (b = 0), M_tc = 1.2; and at triaxial extension (b = 1), M_te = 0.86, resulting in k = 0.71.

The irregular hexagon of Mohr-Coulomb (1776) has sharp corners (Figure 2), violating the third condition for an ideal stress ratio function. These sharp corners are evident when viewing how the model approaches the triaxial extension and compression limits at π/6 and −π/6 in Figure 3 as well.

One simple example of a stress ratio function entirely free of sharp corners is the Drucker-Prager (1952) criterion which appears as a circle on π-plane (Figure 2), indicating a constant stress ratio in Figure 3. This criterion is a cone in the principal stress space centred around the hydrostatic axis, which means that material yielding is dependent on p and q, but not θ. The Drucker-Prager criterion forces M_te to be equal to M_tc, which does not follow the trend established by data as shown in Figure 3.

The issue of finding a suitable cross-section free of sharp corners has been addressed by multiple researchers (Gudehus, 1973; Argyris et al., 1974; Willam and Warnke, 1974; Jiang and Pietruszczak, 1988) who created smooth functions in the π-plane as shown in Figure 4. These functions vary in their mathematical complexity but only use one or two parameters including k which allows calibration to M_te as well within the confines of convexity; a capability the Mohr-Coulomb equation does not have.

Gudehus (1973) and Argyris et al. (1974) proposed a simple equation with smooth transitions across meridians but as demonstrated in Figure 3 it overshoots nearly all the available data. It also requires k ≥ 0.778 to remain convex (Table 1), which is near the upper limit of the range of k for typical soils and is too high a practical limit (Cornforth, 1961; Cornforth, 1964). The equation proposed by Willam and Warnke (1974) is fairly complicated, albeit differentiable, but it acts as an upper bound to the available data. Jiang and Pietruszczak’s (1988) equation is fairly simple and differentiable, but also crosses the upper range of the data in Figure 3.

Approximating the Mohr-Coulomb (1976) criterion has not been the only path for generalizing constitutive models to three-dimensional stress conditions. Lade and Duncan (1975) described yielding based on simple combinations of principal stresses and a single material parameter κ₁, which is used to anchor the yield surface to the reference triaxial compression condition (Lade and Musante, 1977). The parameter κ₁ affects both the size and shape of the yield surface, and the curve cannot be adjusted to fit triaxial compression and extension data simultaneously. The Lade-Duncan yield surface was described initially using only combinations of principal stresses, but Yang et al. (2006) were able to present it using the stress invariants. The function does not fit the data presented in Figure 3, even though it was created as an equation that reproduced Lade and Duncan’s (1975) “failure surface” data. This discrepancy highlights the importance of separating other influences, such as hardening and dilatancy, from yield or plastic potential functions by, for example, using the critical state as a reference, as was done here.

Van Eekelen (1980) developed a three-parameter stress ratio function and examined its applicability to Lade and Duncan’s (1975) data and convexity using various combinations of parameters. By examining the results, Van Eekelen then simplified the equation down to one parameter: α. However, this parameter must be numerically computed after anchoring the function to triaxial compression, much like with the Lade-Duncan surface. Similar to other one-parameter models (e.g. Drucker-Prager (1952) and Lade and Duncan (1975)) when a value of M_tc is selected as an anchor, the rest of the yield surface is automatically determined with no control over the range of the function including the extension condition. When α is selected, k is automatically determined by whatever value it takes. In this case, the model does reasonably well over positive Lode angles (Figure 3), but underestimates behaviour in extension.

The Lade-Duncan and Van Eekelen stress ratio functions are plotted on the π-plane in Figure 5, with their parameters selected for the case of k = 0.75. The two functions nearly coincide at this selected value of k, but take on other characteristic shapes as their calibration parameter is adjusted. As k decreases, the Lade-Duncan surface takes on more of a circular shape, while the Van Eekelen surface becomes similar to a smoothed triangle. Both surfaces are differentiable and show smooth transitions over the compression and extension meridians.

Matsuoka and Nakai (1974) developed a stress ratio function formulation by expanding the mobilized friction angle framework into principal stress space and identifying two mobilized friction angle parameters that represent the behaviour in the direction of the intermediate principal stress. Their equation does not have a closed-form solution, which makes it difficult to implement and its shape a function of M_tc, in a similar fashion to Mohr-Coulomb, as discussed earlier.

Based on observations of Cornforth’s (1961, 1964) plane strain tests on Brasted sand, Jefferies and Shuttle (2011) noted that, for soil, the Mohr-Coulomb (1776) equation underestimates the true relation of stress ratio and Lode angle while the Matsuoka and Nakai (1974) formula overestimates this relation (see Figure 3). While it was pointed out that an averaging scheme could be implemented, they proposed a simple closed-form equation which is closer to the data near triaxial compression and reproduces the Mohr-Coulomb (1776) M_te, but it is not smooth at the compression meridian and thusly not differentiable throughout the entire Lode angle domain. The Matsuoka-Nakai and Jefferies-Shuttle stress ratio functions are plotted on the π-plane in Figure 6.

It is readily apparent from Figure 3 and Table 1 that many of the presented curves for the stress ratio do not follow experimental data from laboratory tests. Four of the functions (Mohr-Coulomb (1776), Drucker-Prager (1952), Matsuoka and Nakai (1974), and Jefferies and Shuttle (2011)) have M_tc as their only calibration parameter and thusly do not allow M_te to be calibrated independently of M_tc. Convexity ranges between unconditional to conditions that restrict calibration parameters (e.g. k). In the case of Gudehus (1973) and Argyris et al. (1974), the limit on k may be too restrictive, limiting the ability of the equation to capture some observed material behaviour. Some of the presented functions have additional difficulties such as complicated or implicit equations, or having sharp corners (not differentiable) on π-plane. A more appropriate stress ratio function should be free of these complications while passing through the available data and allowing for the definition of M(θ) at both compression and extension meridians.

Equation 5, which is a sinusoidal function sin(3θ) scaled and shifted using the ratio, k, fulfils most of the listed criteria albeit with a high convexity limit of k ≥ 0.819 and is a poor match for the available data, as shown in Figure 7.

Mathematically, sin(3θ) is a simple equation that is easy to implement and has the familiar and easily interpretable shape of a sinusoid, with flat ends at ±π/6. To correct for the poor matching of experimental data while ensuring the other conditions remain satisfied, the sine curve can be stretched downward by adjusting its argument. One suitable function that pulls the sine wave downward to fit the data better is an arc from a circle as shown in Figure 8. The arc is tied at the endpoints of the Lode angle domain to ensure that the function will remain flat at θ = ±π/6.

Implementing the circular arc into the sine function argument produces Equation 6:

a is the positive value of the coordinate of the centre of the circle, which is a simple function of R due to the anchoring of the arc to the domain endpoints:

For the proposed function, the derivative at θ = ±π/6 is zero (Equation 7). The reflections on π-plane are smooth as seen in Figure 9.

The convexity condition limits R to [1.51, 6.26] at k = 0.75. Figure 10 illustrates the range of permissible R values as a function of k. The convexity is visible in Figure 9 and reaches its critical values near the extension meridian.

With k = 0.75 and R = 1.51, the stress ratio function proposed is a better fit to the hollow cylinder data of Yang (2013) at zero Lode angle and passes through more of the clustered plane strain data collected by Cornforth (1961, 1964) and Wanatowski and Chu (2007), rather than being an upper bound as many of the existing functions act as for three-dimensional stresses with positive Lode angle. The function is smooth and convex over a wide range of its parameters (k and R), and is simple to implement and has the benefit of being based on understandable geometry. The stress ratio function proposed also has the benefit of flexibility: while R = 1.51 passes through the current data that is available, the dataset is fairly sparse due to a lack of material testing and there is no critical state failure data for negative Lode angles other than at the extension meridian at −30°. This function has the advantage of being calibratable to produce a wide range of other shapes for when more M(θ) vs θ data becomes available in the literature. The convexity limitations on k and R illustrated in Figure 10 may become problematic in some materials, but for the limited data presented here, the function appears to capture experimental data well.

Lack of high-quality experimental data that cover a wide range of Lode angles, especially negative values, and dilatancy conditions is what limits further validation and improvements to this aspect of constitutive modelling. As three-dimensional models become more common in practise, this often taken-for-granted step in development of constitutive models deserves more attention.

To validate the applicability and efficiency of the proposed equation, drained triaxial compression tests were modelled using the NorSand critical state-based constitutive model (Jefferies, 1993; Jefferies and Shuttle, 2005). Triaxial compression conditions are ideal for this purpose because the triaxial compression meridian is where certain stress ratio equations, such as NorSand’s recommended Jefferies and Shuttle (2011) equation, lack continuous derivatives.

The NorSand model was implemented (Mozaffari and Ghafghazi, 2017) in a UMAT user material subroutine in the Finite Element package Abaqus. For three-dimensional stress–strain analysis, the recommended stress ratio function for the NorSand constitutive model is the Jefferies and Shuttle (2011) equation, though NorSand is compatible and has been used with other stress ratios such as the Matsuoka and Nakai (1974) equation and the Gudehus (1973) and Argyris et al. (1974) equation (Cheong, 2006; Cheong et al., 2011; Roy, 2012). For this work, some modifications were added (Liu et al., 2023), and the modified Euler method with automatic error control (Sloan, Abbo and Sheng, 2001) was used for stress integration. The recommended NorSand stress ratio function (Jefferies and Shuttle (2011)) and Equation 6 were implemented as options and compared in otherwise identical multi-element triaxial compression simulations. The axisymmetric computational domain was divided into 800 four-node axisymmetric quadrilateral elements (CAX4). The vertical displacement of the bottom boundary was restricted, while constant pressure was applied on the lateral boundary. The model was sheared by applying a constant vertical displacement on the top boundary, until an axial strain of 10% was reached. The initial conditions and input NorSand parameters were identical to those for dense sand given in Mozaffari and Ghafghazi (2017), as summarized in Table 2. For the stress ratio function in Equation 6, R = 1.51 and k = 0.75 were used in the analysis. The simulated results are compared with the existing VBA solution (Jefferies and Been, 2015) as shown in Figure 11. Ideally, both models should perfectly agree with the single element (Jefferies and Been, 2015) model, given that under triaxial compression M(θ) = M_tc regardless of the equation used. However, the multi-element model using the Jefferies and Shuttle (2011) equation exhibits numerical instability after about 4% axial strain (Figure 11). The instability is also manifested in the deformed shape of the model in Figure 12 by the time the models reach 10% axial strain (Figure 12(c)), while the model using Equation 6 (Figure 12(b)) maintained the expected uniformly deforming shape.

The model using Equation 6 is not only more stable, but also significantly more computationally efficient. When using a stress tolerance (STOL) of 10⁻⁴, Equation 6 only required one iteration to apply an axial strain increment of 0.02% at 4% axial strain while more than 80 iterations were needed with the Jefferies and Shuttle (2011) equation. With a stricter STOL of 10⁻⁶, the required numbers of iterations became 6 and over 9000, respectively. This led to more than 30 times longer computational time with a 4.7 GHz 16-core CPU to run the same triaxial simulation.

Geomechanical constitutive models are almost ubiquitously developed under triaxial compression conditions. They are then generalized into the three-dimensional stress space by making their components, typically various forms of the stress ratio, functions of the Lode angle.

A three-dimensional generalization of a constitutive model is created by using a stress ratio that is a function of Lode angle. A high-quality generalization should ideally match available experimental data, be simple to numerically implement, retain convexity and have a continuous derivative across the entire domain of the Lode angle. Nine popular stress ratio functions were discussed, and it was shown that some do not fit the existing data well because they lack the necessary number of parameters, or their shapes do not allow good fits. Some functions, like the classic Mohr-Coulomb (1776), have sharp corners on π-plane around triaxial conditions, resulting in their derivative and hence gradients not being continuous and causing singularities in their flow rules and consistency conditions. Others have complicated formulations or implicit solutions. Of the existing equations, the equation proposed by Jiang and Pietruszczak’s (1988) satisfies all the conditions laid out as it is fairly simple, differentiable, and robustly convex, but its shape is more in line with the upper bound of the critical state stress ratio data presented.

A new stress ratio function for 3-D generalization was proposed which satisfies all the desired conditions with two parameters: R, to control the shape of the curve, with R = 1.51 expected to be sufficient in most cases, and k as the ratio of triaxial extension to compression stress ratios. The proposed equation is differentiable everywhere, remains convex for a reasonable range of k and R values, and captures the existing data better than most other equations. This new stress ratio function has the benefit of being customizable and maintains numerical stability within a large range of its two parameters so that it can be used even as more data becomes available.

The proposed equation was implemented within the critical state constitutive model NorSand in a multi-element drained triaxial compression test. It was demonstrated that the new equation improves both stability and numerical efficiency, compared to NorSand’s default model equation which does not have a continuous derivative at triaxial conditions.

The datasets analysed during the current study are available from the corresponding author upon reasonable request.

The support from the National Science and Engineering Research Council of Canada (NSERC) (RGPIN-2016-05622 and ALLRP 556818-20), WSP, Conetec and Vale Brazil is acknowledged.

Once familiar with the definition, the human brain can easily detect convex or concave shapes, including those of yield surfaces on the π-plane. Mathematically, the sign of the second derivative cannot be used as a criterion the way it is in the Cartesian space, because of the polar nature of the π-plane. Instead, for the common shapes of yield surfaces that have the origin of the π-plane within them, convexity is checked by drawing a tangent line to the function and checking whether the function falls on the side of the line that is closer to the origin in the neighborhood.

On the π-plane shown in Figure A1, the radial polar coordinate is the value of the yield surface function M(θ) (⁠$2/3q$ per Figure 1) and the angular polar coordinate is the Lode angle θ. The convexity of the yield surface function M(θ) at the Lode angle θ_o is determined by drawing a tangent line g(θ) at θ_o and then determining whether a small angle dθ away, M(θ + dθ) is closer to the origin than the line g(θ + dθ).

The slope of a function, used to create a tangent line, is expressed as dy/dx in Cartesian coordinates, which can be readily converted to polar coordinates, consisting of a radius and a coordinate angle, using Equation A1.

Equation A2 will then define the tangent line g(θ) passing through M(θ_o) with the slope dy/dx.

The convexity condition is then written as in Equation A3:

as illustrated in Figure A1.