It is not allowed to uniformly apply the criteria of Tresca, the Mohr-Coulomb, the Mises, the Lade-Duncan and the Matsuoka-Nakai to different soil types when the stresses are different. Refer to the form of strength criterion given in literatures (Liu et al., 2010; Xiao et al., 2010; Xiao et al., 2011a; Xiao et al., 2011b; Xiao et al., 2012a; Xiao et al., 2012b; Xiao et al., 2012c; Sun et al., 2013), a simple UNS criterion is given. The formulation can be expressed as: q = M g ( θ ) p a ( p + σ 0 p a ) n (1a) or ln ( q p a ) = n ⁡ ln ( p + σ 0 p a ) + ln [ M g ( θ ) ] (1b) where: p = σ 1 + σ 2 + σ 3 3 (2) q = 2 2 ( σ 1 − σ 2 ) 2 + ( σ 2 − σ 3 ) 2 + ( σ 1 − σ 3 ) 2 (3) g ( θ ) = 3 2 L 1 cos ⁡ δ θ + 6 α ⁡ sin ⁡ φ 0 3 − sin ⁡ φ 0 3 2 L 1 cos ⁡ δ 0 + 6 α ⁡ sin ⁡ φ 0 3 − sin ⁡ φ 0 (4) δ θ = 1 3 arccos ( − 3 3 L 2 ⁡ cos ( 3 θ ) 2 L 1 3 / 2 ) δ 0 = 1 3 arccos ( − 3 3 L 2 2 L 1 3 / 2 ) } (5) L 1 = ( 3 − sin ⁡ φ 0 ) 3 3 ( 3 − sin ⁡ φ 0 ) 3 − 81 ( 1 − sin ⁡ φ 0 − sin 2 φ 0 + sin 3 φ 0 ) L 2 = 2 ( 3 − sin ⁡ φ 0 ) 3 27 ( 3 − sin ⁡ φ 0 ) 3 − 729 ( 1 − sin ⁡ φ 0 − sin 2 φ 0 + sin 3 φ 0 ) } (6)

σ ₀ denotes the bonding stress, reflecting the cohesive property exhibited by the cohesive-frictional material, σ ₀ = ccotφ ₀ (where c is cohesion, φ ₀ is friction angle). n denotes the slope in the ln ( q p a ) − ln ( p + σ 0 p a ) plane. ln [ M g ( θ ) ] denotes the intercept, the effect of geomaterial frictional characteristic. The exponent n serves for describing how hydrostatic pressure affects material failure. α, denotes the combination coefficient, and it is able to confirm its value based on triaxial compression strength and any strength when b ≠ 0. Figure 1 illustrates how φ ₀ and α affect the curve change of the linear strength criterion in the deviatoric plane when other conditions are the same. It can be seen from the figure that, in theory, the criterion shape is a circular firstly when the friction angle is 20°, and then changes into a curved triangle when it gradually changes to 50°. In an increasing range of −0.3 to 10, the criterion is firstly a SMP curved triangle and then becomes the Mises circle, which approximates all the curves between the SMP curved triangle and the Mises circle. Specifically, with n = 1 and α = 0, we adopt the Lade-Duncan criterion as the UNS criterion. With n = 1 and α = −0.3, we adopt the SMP criterion as the UNS criterion. With n = 1 and α → +∞, it is the Extended Mises criterion. Figure 2 illustrates how σ ₀ and n affect the failure curves. Accordingly, parallel failure curves can be obtained if other parameters are equal. Figure 2B illustrates how n affects the failure curves in meridian plane and octahedral plane. We can see expanding failure planes with n increase.

In the consolidated drained true triaxial strength experiments performed by Toyota et al. (2004) on saturated silty sand, they set the hydrostatic pressures to values of 100, 200, and 300 kPa. In the process of loading, in a given experiment, b held a constant value of 0, 0.25, 0.5, 0.75, and 1.0. Figure 3 shows the experimental results, and as calculated by the proposed criterion, σ ₀ = 0 kPa, M = 1.509, and n = 1, α = 0.229. Figure 3 also shows the curves under theoretical prediction, which well fit experimental data, showing that the developed criterion is capable of well describing silty sand strength behavior, namely the variations of b and p. In the deviatoric plane, when the hydrostatic pressure is the same, the shear strength shows a downtrend when b increases, showing how the intermediate principal stress affects the shear strength (Figure 3B).

In general, there are many factors that decide soil yielding and soil strength, namely the current stress state, the loading history, the internal structure and the loading direction regarding internal structure as well. Also, the yield function shall be objective and not affect or be affected by the coordinate system. In this paper, the yield function is proposed following for describing soil yielding when the intermediate principal stress is taken into account (Gao, 2012): f = q − H g ( θ ) p = 0 (7)

The above yield function is modified from the UNS criterion proposed this paper, then a hardening parameter H is used to replace the originally constant frictional coefficient M. The other parameters (p, g(θ), q) are shown in Eq. 2 and 3 and Eq. 4, respectively. In the true triaxial tests, on the deviatoric plane, θ denotes the angle between the current stress state and the vertical stress axes, Figure 4A partitions the deviatoric plane; Figure 4B is the yield surface; Figure 4C is the yield loci with different hardening parameter values.

The H evolution law is proposed as follows: d H = 〈 d L 〉 r H = 〈 d L 〉 G c h H p a ( M f − H ) (8)

In the equation, r _H shows the H evolution direction, and r _H ≥ 0; dL is a loading index. <x> = 0 is the Macauley bracket, and x ≤ 0 when x > 0, and <x> = x when x > 0; c _h denotes a positive constant.

Dilatancy relation supports the soil constitutive model. The dilatancy is described by the dilatancy relation found in the literature of Li and Dafalias (2004): D = d ε v p 2 / 3 d e i j p d e i j p = d 1 exp ( ∫ 〈 d L 〉 ) ( M p − H ) (9)

In the equation, d ε v p shows the plastic volumetric strain increment; the plastic deviatoric strain increment is denoted as d e i j p ( = d ε i j p − d ε v p δ i j / 3 ) ; d ₁ denotes a positive parameter of model; The phase transformation stress ratio is expressed as M _p which is obtained from the conventional triaxial compression tests about remolded samples. The denominator in Eq. 9 serves for controlling the change of volume, particularly in the case of high strain level. With the sample being sheared to critical state, the plastic deviatoric strain increase can be unlimited. Hence, when the flow is unlimited, the denominator term can be infinite, thus the value of D goes to 0.

Besides, considering the yield function of Eq. 7, we propose the associated flow rule: d e i j p = 〈 d L 〉 n i j (10)

In Equation 10, n _ij denotes a unit tensor, which is defined by: n i j = 1 C ( ∂ f ∂ σ i j − 1 3 ∂ f ∂ σ m n δ m n δ i j ) (11) where C can be explained as the quantity norm in the parentheses in Eq. 11. Of particular note is that Pietruszczak (1999) also has adopted a similar flow rule.

By applying the consistency condition to the general yield function form in Eq. 7, we can obtain: d f = ∂ f ∂ σ i j d σ i j + ∂ f ∂ H d H = 0 (12a) or: d f = ∂ f ∂ σ i j d σ i j + 〈 d L 〉 ∂ f ∂ H r H = ∂ f ∂ σ i j d σ i j − 〈 d L 〉 K p = 0 (12b)

K _p is the plastic modulus, K p = − ∂ f ∂ H r H (13)

The constitutive relations can be derived by following the classical plasticity theory. Regarding the deviatoric and volumetric elastic strain increments, the present model makes a postulation on the d e i j e and d ε v e , an isotropic hypoelastic relation: d e i j e = d s i j G d ε v e = d p K } (14) wherein the expressions regarding the elastic moduli K and G are employed according to: K = E 3 ( 1 − 2 ν ) G = E 2 ( 1 + ν ) } (15)

Regarding the relevant plastic strain increments, the dilatancy equation together with the flow rule assist in obtaining the relations as follows: d e i j p = 〈 d L 〉 n i j (16) and d ε v p = D 2 / 3 d e i j p d e i j p = 〈 d L 〉 2 / 3 D (17)

Equation 11 already gives the definition of n _ij . Calculation of the incremental stress–strain relation can be achieved with the help of equations mentioned above. Then, d σ i j = d s i j + d p δ i j = 2 G d e i j e + K d ε v e δ i j = 2 G ( d e i j − d e i j p ) + K ( d ε v − d ε v p ) δ i j = 2 G ( d e i j − 〈 d L 〉 n i j ) + K ( d ε v − 〈 d L 〉 2 / 3 D ) δ i j (18)

Eq. 18 assumes an additive decomposition regarding the total strain increment d ε i j = d ε i j e + d ε i j p . Eq. 19 gives the loading index expression on the basis of Eqs 12–18, 〈 d L 〉 = 2 G ∂ f ∂ σ i j + ( K − 2 3 G ) ∂ f ∂ σ k l δ k l δ i j ∂ f ∂ σ i j ( 2 G n i j + 2 3 K D δ i j ) + K p d ε i j 2 / 3 D = Θ i j d ε i j (19)

The chain rule serves for obtaining ∂ f ∂ σ i j = ∂ f ∂ q ∂ q ∂ σ i j + ∂ f ∂ p ∂ p ∂ σ i j + ∂ f ∂ g ( θ ) ∂ g ( θ ) ∂ θ ∂ θ ∂ σ i j (20) where ∂ f ∂ q = 1 (21) ∂ q ∂ σ i j = ( q 2 − I 1 2 ) ∂ I 1 ∂ σ i j + 9 ∂ I 3 ∂ σ i j 2 q ( q − I 1 ) (22) ∂ f ∂ p = − H g ( θ ) (23) ∂ p ∂ σ i j = 1 3 δ i j (24) ∂ f ∂ g ( θ ) = − H p (25) ∂ g ( θ ) ∂ θ = − 3 3 2 m l 3 r 0 1 − n 2 sin ( δ θ ) sin ( 3 θ ) l 2 3 / 2 ⁡ cos ( δ θ ) (26) m = 3 2 l 1 (27) n = − 3 3 l 2 2 l 1 3 / 2 cos ( 3 θ ) (28)

Combining Eqs 18, 19, the obtained constitutive relation regarding incremental form serves for the following numerical computations. d σ i j = Λ i j k l d ε k l (29) where Λ i j k l = G ( δ i k δ j l + δ i l δ j k ) + ( K − 2 3 G ) δ i j δ k l − h ( d L ) ( 2 G n i j + 2 3 K D δ i j ) Θ k l (30)

The Heaviside step function is expressed as h (dL), and h (dL > 0) = 1, h (dL ≤ 0) = 0.