1 Introduction

Coaxiality between principal stresses and principal plastic strain rates was first postulated by Saint–Venant [27]. If we follow the theories of stable equilibriumFootnote 1 phenomena according to Gibbs [11], systems that are in stable equilibrium conditions, when slightly perturbed, have the tendency to maximize dissipation (or conversely to minimize internal energy). For a given stress state and plastic strain increment, the plastic dissipation is higher when they are coaxial. In the absence of constraints, irregularities and non-homogeneities, coaxiality seems a logical assumption. However, maximization of dissipation is subjected to, for example, kinematic constraints in the medium. In granular materials, such constraints may arise due to anisotropy, non-homogeneity and bifurcation. It has been pointed out by Hill [17] that for anisotropic material, generally principal stresses and principal (plastic) strain increments are non-coaxial except for the special case where the principal stress axes coincide with the axes of anisotropy. In fact, in granular materials, non-coaxiality between principal stresses and principal (plastic) strain rates has been observed through various techniques. For example, Drescher and De Josselin de Jong [9] studied the deformation behavior of a photo-elastic disk assembly to verify the double-sliding free-rotating model [7] and they were able to calculate the degree of non-coaxiality between the axes of principal stresses and strain rates. Roscoe et al. [25], using simple shear tests, observed that principal stresses and principal plastic strain rates can be non-coaxial. Using the directional shear cell apparatus (DSC) [2], Arthur et al. [3] investigated the stress–strain behavior of the Leighton-Buzzard Sand samples due to change of stress path direction and found that principal stresses and principal strain rates are generally non-coaxial. Gutierrez et al. [15] employed the hollow cylinder apparatus to investigate deformation behavior of dense air-pulviated Toyoura Sand subjected to proportional stress path, pure principal stress rotation and loading with increasing deviatoric stress combined with principal stress rotation. Their test results show that loading conditions that involve principal stress rotations are in general non-coaxial. The hollow cylinder has been popularly applied to the investigation of stress–strain behavior of soils under loading conditions that involve principal stress rotation e.g., [3, 5, 29, 35]. Non-coaxiality has also been observed in discrete element model (DEM) setups [1, 3, 29, 35].

1.1 Definition

Consider a stress state in simple shear conditions as shown in Fig. 1. Let the direction of the minor principal stress and the minor principal strain rate direction make angles \(\alpha_{\sigma }\) and \(\alpha_{{\dot{\varepsilon }}}\), respectively, with the x-axis. When there is no deviation between the two angles, i.e., \(\alpha_{{\dot{\varepsilon }}} = \alpha_{\sigma }\), the principal stresses and the principal strain rates are said to be coaxial; otherwise they are non-coaxial. The degree of non-coaxiality may be defined by the deviation angle \(\Delta = \alpha_{{\dot{\varepsilon }}} - \alpha_{\sigma }\). The definition had been extended (between principal stresses and principal plastic strain rates) to a 3D stress–plastic strain rate space by Gutierrez and Ishihara [12]. In general, if a stress tensor \(\sigma_{ij}\) and a plastic strain rate tensor \(\dot{\varepsilon }_{ij}^{p}\) are defined in a Cartesian coordinate system \({\mathbf{x}} = \{ \begin{array}{*{20}c} {x_{1} } & {x_{2} } & {x_{3} } \\ \end{array} \}^{\text{T}}\) as

$$\sigma_{ij} = \left[ {\begin{array}{*{20}c} {\sigma_{11} } & {\sigma_{12} } & {\sigma_{13} } \\ {\sigma_{21} } & {\sigma_{22} } & {\sigma_{33} } \\ {\sigma_{31} } & {\sigma_{32} } & {\sigma_{33} } \\ \end{array} } \right]\quad {\text{and}}\quad \dot{\varepsilon }_{ij}^{p} = \left[ {\begin{array}{*{20}c} {\dot{\varepsilon }_{11}^{p} } & {\dot{\varepsilon }_{12}^{p} } & {\dot{\varepsilon }_{13}^{p} } \\ {\dot{\varepsilon }_{21}^{p} } & {\dot{\varepsilon }_{22}^{p} } & {\dot{\varepsilon }_{23}^{p} } \\ {\dot{\varepsilon }_{31}^{p} } & {\dot{\varepsilon }_{32}^{p} } & {\dot{\varepsilon }_{33}^{p} } \\ \end{array} } \right],$$
(1)

the principal stress tensor and the principal strain rate tensor are obtained by transforming each as

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\sigma }_{ij} = T_{ik}^{\sigma } \sigma_{kl} T_{lj}^{\sigma } = \left[ {\begin{array}{*{20}c} {\sigma_{1} } & 0 & 0 \\ 0 & {\sigma_{2} } & 0 \\ 0 & 0 & {\sigma_{3} } \\ \end{array} } \right]\quad {\text{and}}\quad \underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{\dot{\varepsilon }}_{ij}^{p} = T_{ik}^{{\dot{\varepsilon }}} \dot{\varepsilon }_{kl}^{p} T_{lj}^{{\dot{\varepsilon }}} = \left[ {\begin{array}{*{20}c} {\dot{\varepsilon }_{1}^{p} } & 0 & 0 \\ 0 & {\dot{\varepsilon }_{2}^{p} } & 0 \\ 0 & 0 & {\dot{\varepsilon }_{3}^{p} } \\ \end{array} } \right],$$
(2)

respectively, such that the off-diagonal terms are zero. \(T_{ik}^{\sigma }\) and \(T_{ik}^{{\dot{\varepsilon }}}\) are, respectively, matrices that transform the stress tensor and the plastic strain rate tensors into their respective principals. If \(T_{ik}^{\sigma } \ne T_{ik}^{{\dot{\varepsilon }}}\), then the principal stresses and principal strain rates are said to be non-coaxial, Gutierrez and Ishihara [12]. From now on, the condition is referred to as non-coaxiality.

Fig. 1
figure 1

Non-coaxial behavior in simple shear tests (Roscoe et al. [25])

Full size image

1.2 Trend

Prior to bifurcation, experimental evidences confirm that the degree of non-coaxiality decreases in magnitude with increasing stress ratio, see e.g., [3, 24]. In Fig. 2, test results by Arthur et al. [3] are presented. The test results show that degree of non-coaxiality vanishes with stress ratio.

Fig. 2
figure 2

Evolution of degree of non-coaxiality with stress ratio when the principal stress \(\sigma_{1B}\) is preceded by \(\sigma_{1A}\). a DEM simulations. b DSC experiments (after Arthur et al. [3])

Full size image

The tendency of the degree of non-coaxiality during post-bifurcation deformation is controversial as also pointed out in Tejchman and Wu [28]. For example, Vardoulakis and Georgopoulos [34] presented a biaxial test on Karlsruhe Sand that shows that the degree of non-coaxiality vanishes with shear strain after bifurcation. Similarly, DEM simulations in Thornton and Zhang [29] reveal that the degree of non-coaxiality vanishes toward the critical state regardless of post-peak reduction in stress ratio. On the other hand, Gutierrez and Vardoulakis [13] presented tests on Nevada Sand in which “the post-bifurcation non-coaxiality parameter varies with shear displacement.” The variation in their plot implies that degree of non-coaxiality increases with shear strain during post-bifurcation deformation. This tendency has been attributed mainly to rotation of principal stresses because of “a simple shear loading condition imposed in the shear band.” Tejchman and Wu [28] carried out a numerical investigation of shear localization in dilatant bodies using a micro-polar hypoplastic model with a focus on non-coaxiality and stress–dilatancy behavior of an initially medium dense Karlsruhe Sand. In the same paper, they showed an increase in the degree of non-coaxiality during post-bifurcation deformation. However, they obtained only a small rate of increase (from about 1° at the peak to about 6° at large deformations). Although it is small, the trend does not agree with the observed trends in experimental measurements in Vardoulakis and Georgopoulos [34] on the same sand.

Thornton and Zhang [29] from their DEM simulations pointed out that “at any stage of shearing, during simple shear deformation, the angle of non-coaxiality depends on the mobilized angle of shearing resistance, the rate of dilation, the initial stress state, and the applied loading path.” Furthermore, Thornton and Zhang [29] concluded that when there is no further change in volume, the stress and strain rate directions are coaxial.

1.3 Effect

It is not fully understood yet in what way non-coaxiality affects the deformation behavior of soils. However, one aspect can be clearly envisaged, that is, non-coaxiality affects the energy dissipation mechanism of the medium in some way. Studies also show that dilatancy is strongly influenced by non-coaxiality of principal stresses and principal plastic strain rates. The possible influence of non-coaxiality on stress–dilatancy behavior of a granular medium was first pointed out by De Josselin de Jong [8]. Gutierrez and Ishihara [12] proposed a non-coaxial version of Taylor’s work hypothesis. Later, Gutierrez and Wang [14] introduced a degree of non-coaxiality into Rowe’s stress–dilatancy relation. Closely investigating these theories, some inconsistencies were noticed (Tsegaye [30], Tsegaye et al. [33]). In the next sections, non-coaxial plastic dissipation and stress–dilatancy relations are developed by extending the theory—proposed by Tsegaye and Benz [31].

The following applies thought out the paper:

  1. 1.

    Stress quantities, friction angles and cohesion are always taken to be effective values without any special indication by a prime.

  2. 2.

    Strain rates defined here refer generally to an artificial time increment and can be considered as simultaneous infinitesimal strain increments.

2 Non-coaxial plastic dissipation and stress–dilatancy

For an isothermal condition, the energy variation may be written as

$$\dot{F} + {\mathcal{D}} - \sigma_{ij} \dot{\varepsilon }_{ij} = 0,\quad {\mathcal{D}} \ge 0$$
(3)

where \(\dot{F}\) is the rate of Helmholtz free energy, \({\mathcal{D}}\) is the rate of dissipation, \(\sigma_{ij}\) is Cauchy’s stress tensor, and \(\dot{\varepsilon }_{ij}\) is the strain rate tensor. The strain rate is additively decomposed into elastic and plastic, i.e., \(\dot{\varepsilon }_{ij} = \dot{\varepsilon }_{ij}^{e} + \dot{\varepsilon }_{ij}^{p}\) according to the well-known hypothesis of elastoplasticity. Accordingly, the rate of work may be decomposed into elastic and plastic, i.e., \(\tilde{W} = \tilde{W}^{e} + \tilde{W}^{p} = \sigma_{ij} \dot{\varepsilon }_{ij}^{e} + \sigma_{ij} \dot{\varepsilon }_{ij}^{p}\). The Helmholtz free energy may be decomposed into elastic and plastic under certain assumptions as \(\dot{F} = \dot{F}^{e} + \dot{F}^{p}\) where \(\dot{F}^{e} = \sigma_{ij} \dot{\varepsilon }_{ij}^{e}\) and hence \(\dot{F}^{p} : = \dot{F} - \sigma_{ij} \dot{\varepsilon }_{ij}^{e}\). However, in the following \(\dot{F}^{p} = 0\) is assumed such that \({\mathcal{D}} = {\mathcal{D}}^{p} = \sigma_{ij} \dot{\varepsilon }_{ij}^{p}\). \({\mathcal{D}}^{p}\) is referred to as plastic dissipation.

Next, we will present a theoretical framework that was developed [30] for describing stress–dilatancy relations and plastic dissipations in geomaterials. We begin to lay down the theory first in the plane strain and in the axisymmetric condition, and then we will continue to apply the same approach for establishing plastic dissipation and stress–dilatancy relation considering the full stress and plastic strain rate tensors.

2.1 Plane strain and axisymmetric conditions

Assuming coaxiality, for axisymmetric and plane strain conditions, Tsegaye and Benz [31] derived a plastic dissipation of the form

$${\mathcal{D}}_{N}^{p} = {{r_{1} \sigma_{1} \dot{\varepsilon }_{1}^{p} \left( {C_{N} - 1} \right)} \mathord{\left/ {\vphantom {{r_{1} \sigma_{1} \dot{\varepsilon }_{1}^{p} \left( {C_{N} - 1} \right)} {C_{N} }}} \right. \kern-0pt} {C_{N} }}.$$
(4)

where \(\sigma_{1}\) is the major principal stress and \(\dot{\varepsilon }_{1}^{p}\) is the major principal plastic strain rate (along \(\sigma_{1}\) in this case). \(r_{1}\) are defined such that \(r_{1} = 1\), for plane strain and for triaxial compression and \(r_{1} = 2\) for triaxial extension, \(C_{N}\) is a stress ratio at the critical state and hence a constant.

For \(\sigma_{1} \dot{\varepsilon }_{1}^{p} \ge 0\) and, where \(\varphi_{\mu }\) is the interparticle friction angle, the original stress–dilatancy relationship of Rowe’s [26] and the plastic dissipation thereof is found. However, as discussed in Tsegaye and Benz [31], \(C_{N}\) does not have to be the one given in Rowe [26] and the stress ratio \(N_{\sigma }\) does not have to obey the Mohr–Coulomb criterion. The advantage of this approach has been demonstrated in Tsegaye and Benz [31] by extending existing stress–dilatancy equations and deriving a stress–dilatancy relationship for the Hoek–Brown criterion [18–21] which in its preliminary validation turns out to be in a good agreement with Farmer’s [10] servo controlled triaxial compression tests on sand stone and mudstone. Assuming \(C_{N} = C_{N}^{U}\) when \(\sigma_{1} \dot{\varepsilon }_{1}^{p} < 0\), the condition \(0 < C_{N}^{U} = {1 \mathord{\left/ {\vphantom {1 {C_{N}^{L} }}} \right. \kern-0pt} {C_{N}^{L} }} \le 1\) guarantees non-negative plastic dissipation. The suitability of the corresponding stress–dilatancy equation for soil models intended for the modelling of deformation behaviour of soils under cyclic loading has been demonstrated in Tsegaye [30]. The assumption of coaxiality has further been lifted in Tsegaye [30] and Tsegaye et al. [33] for axis symmetric and plane strain conditions. In Tsegaye [30], the theory is extended such that a full 3D stress–strain condition is considered. To give the full treatment of the theory, the approach in Tsegaye [30] and Tsegaye et al. [33] for axis symmetric and plane strain conditions is presented next.

For non-coaxial principal stresses and principal plastic strain rates, the rate of plastic work in the principal stress space may be written as

$$\bar{W}^{p} = r_{1} \sigma_{1} \dot{\bar{\varepsilon }}_{1}^{p} + r_{3} \sigma_{3} \dot{\bar{\varepsilon }}_{3}^{p} ,$$
(5)

where \(\sigma_{i}\) are principal stresses (\(i = 1 - \,{\text{major}},\,\,i = 3 - \,{\text{minor}}\)) and \(\dot{\bar{\varepsilon }}_{i}^{p}\) are the conjugate plastic strain rate components projected along the principal stress components (coaxial components). \(r_{i}\) are defined such that \(r_{1} = r_{3} = 1\), for plane strain, \(2r_{1} = r_{3} = 2\), for triaxial compression and \(r_{1} = 2r_{3} = 2\) for triaxial extension.

For plane strain and axisymmetric conditions, the relation

$$\left\{ {\begin{array}{*{20}c} {\dot{\bar{\varepsilon }}_{1}^{p} } \\ {\dot{\bar{\varepsilon }}_{3}^{p} } \\ \end{array} } \right\} = \frac{1}{{m_{s} }}\cos^{2} \Delta \left[ {\begin{array}{*{20}c} {m_{s} } & {m_{s}^{2} \tan^{2} \Delta } \\ {\tan^{2} \Delta } & {m_{s} } \\ \end{array} } \right]\left\{ {\begin{array}{*{20}c} {\dot{\varepsilon }_{1}^{p} } \\ {\dot{\varepsilon }_{3}^{p} } \\ \end{array} } \right\},$$
(6)

is proposed for projecting principal plastic strain rate components, \(\dot{\varepsilon }_{i}^{p}\), onto the respective principal stresses. In Eq. (6), \(\Delta\) is degree of non-coaxiality between the principal stress and the principal plastic strain rate directions and \(m_{s} = {{r_{3} } \mathord{\left/ {\vphantom {{r_{3} } {r_{1} }}} \right. \kern-0pt} {r_{1} }}\). The mapping presented in Eq. (6) is actually derived for a plane strain condition, i.e., for \(r_{i} = 1\). The extension to axis symmetric conditions is assuming that axis symmetric conditions can be constructed from superposition of two perpendicular planes strain conditions. From the tests in Gutierrez et al. [15], for monotonic tests with fixed principal stress path, the deviation between the directions of the principal stresses and the principal plastic strain increments is small compared to tests that involve principal stress rotation. Therefore, a coaxial condition may be assumed for axis symmetric conditions.

Next, let the principal stress components and the principal plastic strain rate components obey the relations

$$\sigma_{1} = N_{\sigma } \sigma_{3} \,{\text{and}}\,\dot{\varepsilon }_{3}^{p} = - N_{\psi } \dot{\varepsilon }_{1}^{p} ,$$
(7)

respectively, where \(N_{\sigma }\) is the stress ratio and \(N_{\psi }\) is the dilatancy ratio, respectively, and they are called stress–dilatancy conjugates in Tsegaye [30]. Considering Eqs. (5), (6) and (7), the non-coaxial plastic rate of work can be conveniently written as

$$\tilde{W}^{p} = r_{1} \sigma_{1} \dot{\varepsilon }_{1}^{p} \hat{c}_{N} \tilde{d}_{N} ,$$
(8)

in which \(\hat{c}_{N}\) and \(\tilde{d}_{N}\) are convenient functionals that are established by rearranging the plastic rate of work from Eqs. (5) and (6). The function \(\tilde{d}_{N}\) contains the stress–dilatancy conjugates according to

$$\tilde{d}_{N} \mathop = \limits^{1} 1 - m_{s} \frac{{\tilde{N}_{\psi } }}{{N_{\sigma } }},\quad m_{s} \tilde{N}_{\psi } \mathop = \limits^{2} - m_{s} \frac{{\dot{\bar{\varepsilon }}_{3}^{p} }}{{\dot{\bar{\varepsilon }}_{1}^{p} }} = \frac{{m_{s} N_{\psi } - \tan^{2} \Delta }}{{1 - m_{s} N_{\psi } \tan^{2} \Delta }}$$
(9)

The function \(\hat{c}_{N}\) is dilatancy-coaxiality function given by

$$\hat{c}_{N} = \cos^{2} \Delta - \frac{{m_{s} \tilde{N}_{\psi } + \tan^{2} \Delta }}{{1 + m_{s} \tilde{N}_{\psi } \tan^{2} \Delta }}\sin^{2} \Delta .$$
(10)

Postulating the first variation of the function \(\tilde{d}_{N}\) in Eq. (9)1 to vanish, i.e., the variation

$$\delta \tilde{d}_{N} = - \frac{{\delta (m_{s} \tilde{N}_{\psi } )}}{{N_{\sigma } }} + m_{s} \tilde{N}_{\psi } \frac{{\delta N_{\sigma } }}{{N_{\sigma }^{2} }} = 0$$
(11)

yields the relation,

$$C_{N} m_{s} \tilde{N}_{\psi } = N_{\sigma } .$$
(12)

Note that \(m_{s} \tilde{N}_{\psi }\) contains the degree of non-coaxiality according to Eq. (9)2. The non-coaxial plastic dissipation is obtained by substituting Eq. (11) into Eq. (8)

$${\mathcal{D}}_{N}^{p} = \tilde{W}^{p} = r_{1} \sigma_{1} \dot{\varepsilon }_{1}^{p} \hat{c}_{N} \left( {1 - \frac{1}{{C_{{\tilde{N}}} }}} \right) \ge 0,$$
(13)

Assuming \(C_{{\tilde{N}}} = C_{{\tilde{N}}}^{L}\) when \(\sigma_{1} \dot{\varepsilon }_{1}^{p} \ge 0\) and \(C_{{\tilde{N}}} = C_{{\tilde{N}}}^{U}\) when \(\sigma_{1} \dot{\varepsilon }_{1}^{p} < 0\), guaranteeing the inequality in Eq. (13) requires a bit more than the inequality \(0 < C_{{\tilde{N}}} = C_{{\tilde{N}}}^{U} = {1 \mathord{\left/ {\vphantom {1 {C_{{\tilde{N}}}^{L} \le 1}}} \right. \kern-0pt} {C_{{\tilde{N}}}^{L} \le 1}}\), it also requires \(\hat{c}_{N} \ge 0\) which limits the range of the degree of non-coaxiality \(\Delta \in [ - {\pi \mathord{\left/ {\vphantom {\pi 4}} \right. \kern-0pt} 4},{\pi \mathord{\left/ {\vphantom {\pi 4}} \right. \kern-0pt} 4}]\). This range is more relaxed than the range proposed by de Josselin de Jong [7] who obtained, for a given angle of shearing resistance, \(\varphi\), \(\Delta \in [ - {\varphi \mathord{\left/ {\vphantom {\varphi 2}} \right. \kern-0pt} 2},{\varphi \mathord{\left/ {\vphantom {\varphi 2}} \right. \kern-0pt} 2}]\).

The integration constant, \(C_{N}\), may be established by considering the phase transformation point, i.e., \(m_{s} N_{\psi } \to 1\) and hence \(m_{s} \tilde{N}_{\psi } \to 1\) in the sense of loading.

Substituting Eq. (9)2 into Eq. (12), the non-coaxial stress–dilatancy relationship for axisymmetric and plane strain conditions is written as [30]

$$m_{s} N_{\psi } = \frac{{N_{\sigma } + C_{N} \tan^{2} \Delta }}{{C_{{\tilde{N}}} + N_{\sigma } \tan^{2} \Delta }}.$$
(14)

Note that for a coaxial condition, i.e., when \(\Delta = 0\), the relationship in Eq. (14) simplifies to \(C_{N} m_{s} N_{\psi } = N_{\sigma }\) and \(\hat{c}_{N} = 1\), leading to the plastic dissipation in Eq. (4).

2.2 General stress–strain condition

Next, the theoretical framework is extended to the general stress–strain conditions by considering the rate of plastic work given by

$$\tilde{W}^{p} = \sigma_{ij} \dot{\varepsilon }_{ij}^{p} ,$$
(15)

where the full stress tensor \(\sigma_{ij}\) and the full plastic strain rate tensor \(\dot{\varepsilon }_{ij}^{p}\) are defined with respect to a common reference axes. The stress and the plastic strain rate tensors can be equivalently written as

$$\sigma_{ij} = p\delta_{ij} + s_{ij} \quad {\text{and}}\quad \dot{\varepsilon }_{ij}^{p} = \frac{1}{3}\dot{\varepsilon }_{v}^{p} \delta_{ij} + \dot{e}_{ij}^{p},$$
(16)

respectively, where \(p\) is the effective confining pressure defined by one-third of the trace of the effective stress tensor, i.e., \(p = {{\sigma_{ij} \delta_{ij} } \mathord{\left/ {\vphantom {{\sigma_{ij} \delta_{ij} } 3}} \right. \kern-0pt} 3}\), \(\dot{\varepsilon }_{v}^{p}\) is the plastic volumetric strain rate tensor defined by the trace of the plastic strain rate tensor, \(\dot{\varepsilon }_{v}^{p} = \dot{\varepsilon }_{ij}^{p} \delta_{ij}\), which is valid for infinitesimal strain assumption, \(s_{ij}\) is the deviatoric stress tensor, and \(\dot{e}_{ij}^{p}\) is the deviatoric plastic strain rate tensor.

In the principal stress space described by a set of eigenvectors, the deviatoric stress tensor can be written as

$$\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{s}_{k} = \frac{2}{3}q\sin \{ \theta_{\sigma } + {{2(3 - k)\pi } \mathord{\left/ {\vphantom {{2(3 - k)\pi } 3}} \right. \kern-0pt} 3}\} ,\,\,k = 1,\,\,2,\,\,3$$
(17)

wherein \(q = \sqrt {{{3s_{ij} s_{ij} } \mathord{\left/ {\vphantom {{3s_{ij} s_{ij} } 2}} \right. \kern-0pt} 2}}\) is a stress invariant called deviatoric stress. Similarly, the principal deviatoric plastic strain rate can be written as

$$\dot{\underset{\raise0.3em\hbox{$\smash{\scriptscriptstyle-}$}}{e} }_{k} = \dot{\varepsilon }_{q}^{p} \sin \{ \theta^{\prime}_{{\dot{\varepsilon }}} + {{2(3 - k)\pi } \mathord{\left/ {\vphantom {{2(3 - k)\pi } 3}} \right. \kern-0pt} 3}\} ,\,\,k = 1,\,\,2,\,\,3,$$
(18)

wherein \(\dot{\varepsilon }_{q}^{p} = \sqrt {{{2\dot{e}_{ij}^{p} \dot{e}_{ij}^{p} } \mathord{\left/ {\vphantom {{2\dot{e}_{ij}^{p} \dot{e}_{ij}^{p} } 3}} \right. \kern-0pt} 3}}\) is the deviatoric plastic strain rate.

Considering the transformation in Eq. (2), the deviatoric stress and the deviatoric plastic strain rate are given by

$$s_{ij} = q\tilde{m}_{ij} \quad {\text{and}}\quad \dot{e}_{ij}^{p} = \dot{\varepsilon }_{q}^{p} \tilde{n}_{ij} ,$$
(19)

wherein \(\tilde{m}_{ij}\) and \(\tilde{n}_{ij}\) are given by

$$\tilde{m}_{ij} = \frac{2}{3}[T_{i1}^{\sigma } T_{j1}^{\sigma } \sin (\theta_{\sigma } + {{4\pi } \mathord{\left/ {\vphantom {{4\pi } {3)}}} \right. \kern-0pt} {3)}} + T_{i2}^{\sigma } T_{j2}^{\sigma } \sin (\theta_{\sigma } + {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {3)}}} \right. \kern-0pt} {3)}} + T_{i3}^{\sigma } T_{j3}^{\sigma } \sin \theta_{\sigma }]$$
(20)
$$\tilde{n}_{ij} = T_{i1}^{{\dot{\varepsilon }}} T_{j1}^{{\dot{\varepsilon }}} \sin (\theta_{{\dot{\varepsilon }}}^{p} + {{4\pi } \mathord{\left/ {\vphantom {{4\pi } {3)}}} \right. \kern-0pt} {3)}} + T_{i2}^{{\dot{\varepsilon }}} T_{j2}^{{\dot{\varepsilon }}} \sin (\theta_{{\dot{\varepsilon }}}^{p} + {{2\pi } \mathord{\left/ {\vphantom {{2\pi } {3)}}} \right. \kern-0pt} {3)}} + T_{i3}^{{\dot{\varepsilon }}} T_{j3}^{{\dot{\varepsilon }}} \sin \theta_{{\dot{\varepsilon }}}^{p},$$
(21)

respectively.

Considering Eqs. (15), (16) and (19), the plastic rate of work can then be conveniently written as

$$\tilde{W}^{p} = p\dot{\varepsilon }_{v}^{p} + q\dot{\varepsilon }_{q}^{p} \tilde{m}_{ij} \tilde{n}_{ij} .$$
(22)

Note that Eq. (22) has already been derived in Gutierrez and Ishihara [12] and they have then called the quantity \(\tilde{m}_{ij} \tilde{n}_{ij}\) the non-coaxiality factor, whose absolute value we call here the degree of coaxiality and denote it  with a \(\tilde{c}\).

We proceed to establish the stress–dilatancy relation and the plastic dissipation employing the same techniques we have used for the plane strain and the axis symmetric conditions. Consider constraints of stress ratio and plastic strain rate ratio as

$$q = M_{\sigma }^{\theta } p\quad {\text{and}}\quad \dot{\varepsilon }_{v}^{p} = - M_{\psi }^{\theta } \dot{\varepsilon }_{q}^{p},$$
(23)

respectively, where \(M_{\sigma }^{\theta }\) is the stress ratio and \(M_{\psi }^{\theta }\) is the conjugate dilatancy ratio.

The rate of plastic work per unit bulk volume can be written as in

$$\tilde{W}^{p} = p\dot{\varepsilon }_{q}^{p} \tilde{c}\tilde{d}_{M} ,$$
(24)

where \(\tilde{c} = |\tilde{m}_{ij} \tilde{n}_{ij} |\) and \(\tilde{d}_{M}\) is a function of the stress–dilatancy conjugates and is given as

$$\tilde{d}_{M} = - \tilde{M}_{\psi }^{\theta } + \tilde{M}_{\sigma }^{\theta } ,$$
(25)

where \(\tilde{M}_{\psi }^{\theta } = {{M_{\psi }^{\theta } } \mathord{\left/ {\vphantom {{M_{\psi }^{\theta } } {\tilde{c}}}} \right. \kern-0pt} {\tilde{c}}}\) and \(\tilde{M}_{\sigma }^{\theta } = sM_{\sigma }^{\theta }\), in which the quantity \(s = \text{sgn} (\tilde{m}_{ij} \tilde{n}_{ij} )\) is defined.

Postulating the variation, \(\delta \tilde{d}_{M} = 0\) and hence \(\tilde{d}_{M} = C_{M}^{\theta } \ge 0\), Tsegaye [30] was led to a plastic dissipation

$${\tilde{\mathcal{D}}}_{M}^{p} = p\dot{\varepsilon }_{q}^{p} \tilde{c}C_{M}^{\theta } \ge 0.$$
(26)

The corresponding stress–dilatancy relationship is then

$$M_{\psi }^{\theta } = - \tilde{c}\left( {sM_{\sigma }^{\theta } - C_{M}^{\theta } } \right).$$
(27)

The constant \(C_{M}^{\theta }\) may be established by considering \(M_{\psi }^{\theta } = 0\) at the phase transformation condition as just discussed in Sect. 2.1. Convenient ad hoc functions of void ratio may then be introduced such that the effects of density are taken into account in both the stress–dilatancy relation and the plastic dissipation. A non-constant ad hoc may overrule the postulate that \(\delta \tilde{d}_{M}\) vanishes.

Note that the scalar quantity \(\tilde{m}_{ij} \tilde{n}_{ij}\) contains the degree of non-coaxiality. However, it contains not only the degree of non-coaxiality between principal stresses and principal plastic strain rates but also the difference in the Lode angle of the stress tensor and the plastic strain rate tensor, as it was also discussed in Gutierrez and Ishihara [12]. When the principal stress and the principal plastic strain rate tensor are coaxial, one obtains \(\tilde{m}_{ij} \tilde{n}_{ij} = s\cos (\theta_{\sigma } - \theta_{{\dot{\varepsilon }}}^{p} )\). Notice also that for loading \(s = 1\), for neutral loading \(s: = 0\) and for unloading \(s = - 1\). Sometimes the deviation in the Lode angle has also been inclusively considered as non-coaxiality in the definition of non-coaxiality as a non-coincidence of directions of stresses and plastic strain increments, e.g., as stated in Tejchman and Wu [28].

Note that the plastic dissipation in Eq. (26) and the stress–dilatancy relationship in Eq. (27) are different from the one proposed by Gutierrez and Ishihara’s [12]. Gutierrez and Ishihara postulated that the plastic dissipation obeys \({\tilde{\mathcal{D}}}_{M}^{p} = p\dot{\varepsilon }_{q}^{p} M_{c}\), where \(M_{c}\) is given by the stress ratio at the critical state. Gutierrez and Ishihara [12] themselves asserted that non-coaxiality decreases the plastic dissipation although their final proposition that the plastic dissipation obeys \({\tilde{\mathcal{D}}}_{M}^{p} = p\dot{\varepsilon }_{q}^{p} M_{c}\) does seem not reflect that. The difference in the proposed plastic dissipation is also reflected in the resulting stress–dilatancy relationships as will be shown later.

3 Tendency of the degree of non-coaxiality

During pre-bifurcation deformation of granular materials, various test results consistently show that the degree of non-coaxiality vanishes with the stress ratio. For the plane strain and axis symmetric cases, this is translated into the inequality \(\delta \hat{c}_{N} \ge 0\) in Tsegaye [30] and hence

$$\delta \hat{c}_{N} = - \frac{2}{\cos 2\Delta \tan \Delta }\delta \Delta - \frac{1}{{1 + m_{s} \tilde{N}_{\psi } }}\delta (m_{s} \tilde{N}_{\psi } ) \ge 0.$$
(28)

Then, from Eq. (28) the inequality

$$\delta \Delta \le - \frac{1}{2}c_{\Delta } \tan \Delta \frac{1}{{1 + m_{s} \tilde{N}_{\psi } }}\delta \left( {m_{s} \tilde{N}_{\psi } } \right),$$
(29)

is obtained, where \(c_{\Delta } = \cos 2\Delta\) is the degree of coaxiality. As long as the quantity \(m_{s} \tilde{N}_{\psi }\) is increasing, the magnitude of the degree of non-coaxiality must be decreasing, i.e., when \(\Delta ( + )\), \(\delta \Delta ( - )\) and \(\Delta ( - )\), \(\delta \Delta ( + )\). These inequalities may help establish the evolution rule of the degree of non-coaxiality. A simple evolution rule may be established as

$$\delta \Delta = - \frac{1}{2}\pi^{{a_{\Delta } }} c_{\Delta } \tan \Delta \frac{1}{{1 + m_{s} \tilde{N}_{\psi } }}\delta \left( {m_{s} \tilde{N}_{\psi } } \right),$$
(30)

where \(\pi^{{a_{\Delta } }}\) is added to control the rate of decay, \(a_{\Delta }\) is a material parameter that controls the rate of decay (\(\pi = 3.14 \ldots\) so that \(a_{\Delta }\) takes a relatively smaller value and has no mathematical origin).

Equation (29) implies that the evolution of the degree of non-coaxiality is influenced by dilatancy ratio \(\tilde{N}_{\psi }\) (thus also the stress ratio) and stress path (\(m_{s}\)) in agreement with Thornton and Zhang [29].

Note that when \(\Delta_{0}\) is zero, the angle of non-coaxiality is identically zero for the rest of the deformation irrespective of changes in the dilatancy ratio. Such may be true if the sample is truly isotropic (both in initial fabric and initial stress state). This is in agreement with the results from discrete element simulations, e.g., Thornton and Zhang [29], Wang et al. [35]. However, if there is non-coaxiality due to a constraint of some sort from the beginning, the dilatancy ratio works on it such that in each subsequent increment it tends toward coaxiality.

For \(\delta m_{s} = 0\), i.e., proportional stress path, and \(C_{{\tilde{N}}}\) a constant, Eq. (29) simplifies to

$$\delta \Delta = - \frac{1}{2}\pi^{{a_{\Delta } }} c_{\Delta } \tan \Delta \frac{1}{{C_{{\tilde{N}}} + N_{\sigma } }}\delta N_{\sigma } ,$$
(31)

Equation (31) implies that the degree of non-coaxiality, \(\Delta\), decreases with an increasing stress ratio, \(N_{\sigma }\). This tendency agrees with observations, for example, with Roscoe et al. [25], Matsuoka et al. [23], and Thornton and Zhang [29]. The tendency of Eq. (31) is such that non-coaxiality decreases with positive increment in stress ratio, thus toward the critical state, see the DSC tests and DEM simulations by Arthur et al. [3].

However, if \(\Delta > 0\) and \(\delta N_{\sigma } < 0\), \(\delta \Delta > 0\), Eq. (31) implies that the degree of non-coaxiality increases. This may occur during post-bifurcation deformation (see Gutierrez and Vardoulakis [13]). However, Vardoulakis and Georgopoulos [34] presented a biaxial test on Karlsruhe Sand (Fig. 3) in which the degree of non-coaxiality, during post-bifurcation deformation, shows a tendency of vanishing with shear strain. This tendency seems to contradict the observation of the former, i.e., Gutierrez and Vardoulakis [13]. In explaining the contradiction, Gutierrez and Vardoulakis [13] suggested that the degree of post-bifurcation principal stress rotation may depend on the soil type. It is also noted that the latter tendency cannot be deduced from Eq. (31) at least in its present form, but it may be deduced from Eq. (29). In the test results of Vardoulakis and Georgopoulos [34], the dilatancy angle was still increasing for a while after the peak friction angle (Fig. 3). This result is in agreement with Eq. (29). The seeming contradiction arises if one assumes the evolution of the dilatancy ratio is governed by the stress ratio alone upon which Eq. (31) is derived from Eq. (29).

Fig. 3
figure 3

Evolution of degree of coaxiality and dilatancy during post-bifurcation for Karlsruhe Sand (after Vardoulakis and Georgopoulos [34])

Full size image

Figure 4 shows the evolution of dilatancy-coaxiality function \(\hat{c}_{N}\) and the degree of non-coaxiality \(\Delta\) with stress ratio considering the evolution rule in Eq. (31). The higher the value of \(a_{\Delta }\), the faster the degree of non-coaxiality vanishes. On the basis of how fast the system reaches its maximum dissipation potential under continuous deformation, granular materials may be distinguished into highly dissipation efficient, medium dissipation efficient and low dissipation efficient. Other suggested evolution rules can be found, for example, in Gutierrez and Wang [14] and Gutierrez et al. [16].

Fig. 4
figure 4

Plots of a \(\hat{c}_{N}\) versus \(N_{\sigma }\) and b \(\Delta\) versus \(N_{\sigma }\) for \(\Delta_{0} \approx 40^{0}\) and \(C_{{\tilde{N}}} = 3\)

Full size image

The DEM simulation results in Thornton and Zhang [29] show a strong relationship between the initial degree of non-coaxiality, \(\Delta_{0}\), and the K0 stress state. A possible form that can capture such a relationship between the initial degree of non-coaxiality and the K0 for the simple shear condition is

$$\tan \Delta_{0} = \kappa_{\Delta } \frac{{1 - K_{0}^{2} }}{{K_{0} }},$$
(32)

where \(\kappa_{\Delta }\) is a constant. For the specific simulation results in Thornton and Zhang [29], Fig. 5, \(\kappa_{\Delta } = 2/3\) gives a good fit. However, generally one may speculate that non-coaxiality may arise due to non-homogeneity of the media, stress discontinuity, fabric anisotropy and mechanisms due to bifurcation as existing experimental evidences suggest, e.g., Gutierrez and Vardoulakis [13].

Fig. 5
figure 5

Evolution of the degree of non-coaxiality with shear strain from DEM simulation of a simple shear condition (after Thornton and Zhang [29])

Full size image

4 Application to selected yield criteria

Here, selected yield criteria are considered for illustrating how the theoretical frameworks so far presented can be used to enhance existing stress–dilatancy relationships [30].

4.1 Non-coaxial stress–dilatancy formalism for a Mohr–Coulomb material

Here, we consider a material that obeys Coulomb’s friction rule. As presented in Tsegaye and Benz [31], the mobilized stress ratio and the stress ratio at constant volume are defined by

$$N_{\sigma }^{{\mathrm{MC}}} = N_{\varphi } + bN_{\varphi } - b\quad {\text{and}}\quad \bar{K}_{N}^{{\mathrm{MC}}} = \bar{N}_{c} + b\bar{N}_{c} - b,$$
(33)

respectively, where

$$N_{\varphi } = \frac{{1 + \sin \varphi_{m} }}{{1 - \sin \varphi_{m} }},\bar{N}_{\varphi } = \frac{{1 + f_{sd} \sin \varphi_{c} }}{{1 - f_{sd} \sin \varphi_{c} }},0 \le f_{sd} \sin \varphi_{c} < 1,$$
(34)

And \(b = a/\sigma_{3}\), wherein \(\varphi_{m}\) is mobilized friction angle and \(\varphi_{c}\) is critical state friction angle here considered a material constant, and a is the so-called attraction [22] given as \(a = c\cot \varphi_{c}\) where \(c\) is cohesion, \(f_{sd}\) is an ad hoc function that is introduced to capture effect of effective confining pressure and void ratio. Next, specific stress–dilatancy relations are developed considering loading and unloading following the theoretical framework established in Sect. 2.1.

4.1.1 Loading

Substituting Eq. (33) into Eq. (14), the dilatancy ratio multiplied by the shear mode constant, \(m_{s} N_{\psi }^{{\mathrm{MC}}}\), is given by

$$m_{s} N_{\psi }^{{\mathrm{MC}}} = \frac{{N_{\sigma }^{{\mathrm{MC}}} + \bar{K}_{N}^{{\mathrm{MC}}} \tan^{2} \Delta }}{{\bar{K}_{N}^{{\mathrm{MC}}} + N_{\sigma }^{{\mathrm{MC}}} \tan^{2} \Delta }}.$$
(35)

Then the sine of the mobilized dilatancy angle is obtained from [32],

$$\sin \hat{\psi }_{m}^{L} \mathop = \limits^{1} - \frac{{m_{s} N_{\psi } - 1}}{{m_{s} N_{\psi } + 1}}\mathop = \limits^{2} - c_{\Delta } \frac{{N_{\varphi } - \bar{N}_{c} }}{{N_{\varphi } + \bar{N}_{c} - f_{b} }},\quad f_{b} = {{2b} \mathord{\left/ {\vphantom {{2b} {(1 + b)}}} \right. \kern-0pt} {(1 + b)}}.$$
(36)

Equation (36)2 is obtained by substituting Eq. (35) into Eq. (36)1. For \(b = 0\) Eq. (36)2 reduces to

$$\sin \hat{\psi }_{m}^{L} = - c_{\Delta } \frac{{\sin \varphi_{m} - f_{sd} \sin \varphi_{c} }}{{1 - f_{sd} \sin \varphi_{m} \sin \varphi_{c} }}.$$
(37)

The influence of non-coaxiality on the stress–dilatancy is demonstrated in Fig. 6. As can be seen in the figure, non-coaxiality reduces the magnitude of the dilatancy angle. However, if the rate at which the degree of non-coaxiality vanishes is very high, the effect in the dilative region (at higher mobilizations) can be low. It should be noted here that the stress–dilatancy relation presented in Eq. (37) does not agree with Gutierrez and Wang’s [14] non-coaxial version of Rowe’s stress–dilatancy relation. As can be seen from Eq. (37) and as also illustrated in Fig. 6, the higher the degree of non-coaxiality, the less the magnitude of the dilatancy angle for a given mobilized friction angle, whereas the Gutierrez and Wang [14] non-coaxial version of Rowe’s stress–dilatancy relation gives the opposite, i.e., the higher the degree of non-coaxiality, the higher the magnitude of the dilatancy angle for a given mobilized friction angle below the phase transformation. In fact, when \(c_{\Delta }\) tends to zero, the Gutierrez and Wang [14] formulation tends to given unlimited volumetric contraction and unlimited plastic dissipation.

Fig. 6
figure 6

Effect of non-coaxiality on stress–dilatancy relation in plane strain conditions for different values of \(a_{\Delta }\)(the evolution rule in Eq. (31) considered), \(f_{sd} = 1\) and \(\varphi_{c} = 30\) degree during plastic loading (and the initial degree of non-coaxiality is calculated after Eq. (32) where Jakys formula is used for determining K0)

Full size image

4.1.2 Unloading

For the case of unloading, \(C_{N} = 1 - \bar{K}_{N}^{{\mathrm{MC}}}\) may be inserted into Eq. (36) such that \(m_{s} N_{\psi }^{{\mathrm{MC}}}\) is given by

$$m_{s} N_{\psi }^{{\mathrm{MC}}} = \frac{{\bar{K}_{N}^{{\mathrm{MC}}} N_{\sigma }^{{\mathrm{MC}}} + \tan^{2} \Delta }}{{1 + N_{\sigma }^{{\mathrm{MC}}} \bar{K}_{N}^{{\mathrm{MC}}} \tan^{2} \Delta }},$$
(38)

which yields

$$\sin \hat{\psi }_{m}^{U} \mathop = \limits^{1} - c_{\Delta } \frac{{\bar{K}_{N}^{{\mathrm{MC}}} N_{\sigma }^{{\mathrm{MC}}} - 1}}{{\bar{K}_{N}^{{\mathrm{MC}}} N_{\sigma }^{{\mathrm{MC}}} + 1}}\mathop = \limits^{2} - c_{\Delta } \frac{{(1 + b^{2} )\bar{N}_{c} N_{\varphi } - b(1 + b)(\bar{N}_{c} + N_{\varphi } ) + b^{2} - 1}}{{(1 + b^{2} )(\bar{N}_{c} N_{\varphi } + 1) - b(1 + b)(\bar{N}_{c} + N_{\varphi } )}}.$$
(39)

Equation (39)2 is obtained considering Eqs. (33)1,2 into Eq. (39)1. Considering \(b = 0\), Eq. (39)2 further simplifies to

$$\sin \hat{\psi }_{m}^{U} \mathop = \limits^{1} - c_{\Delta } \frac{{\bar{N}_{c} N_{\varphi } - 1}}{{\bar{N}_{c} N_{\varphi } + 1}}\mathop = \limits^{2} - c_{\Delta } \frac{{\sin \varphi_{m} + f_{sd} \sin \varphi_{c} }}{{1 + f_{sd} \sin \varphi_{m} \sin \varphi_{c} }}.$$
(40)

Comparing Eqs. (37) and (40), \(- c_{\Delta }\) remaining the same, the minus sign in the former changes into a plus sign in the latter. Therefore, the non-coaxial stress–dilatancy relation for the loading and the unloading conditions can be combined as

$$\sin \hat{\psi }_{m} = - c_{\Delta } \frac{{\sin \varphi_{m} - sf_{sd} \sin \varphi_{c} }}{{1 - sf_{sd} \sin \varphi_{m} \sin \varphi_{c} }},$$
(41)

where \(s = 1\) during loading and \(s = - 1\) during unloading.

4.2 Non-coaxial stress–dilatancy formalism for a Hoek–Brown material

A stress–dilatancy formulation for the Hoek–Brown [18–21] failure criterion for rocks and rock masses has been proposed in Tsegaye and Benz [31]. Here, the same stress–dilatancy relationship is extended considering non-coaxiality following the formalism established in Sect. 2.1.

In Tsegaye and Benz [31], the stress ratio, \(N_{\sigma }^{\mathrm{HB}}\), and the modified residual strength derived from the generalized Hoek–Brown criterion were given as

$$N_{\sigma }^{\mathrm{HB}} = 1 + \tilde{b}\left( {{{m_{b} } \mathord{\left/ {\vphantom {{m_{b} } {\tilde{b}}}} \right. \kern-0pt} {\tilde{b}}} + s} \right)^{a} ,\quad \tilde{b} = {{\sigma_{ci} } \mathord{\left/ {\vphantom {{\sigma_{ci} } {\sigma_{3} }}} \right. \kern-0pt} {\sigma_{3} }}$$
(42)

and

$$\bar{K}_{N}^{\mathrm{HB}} = 1 + f_{sd} \tilde{b}\left( {{{\hat{m}_{b} } \mathord{\left/ {\vphantom {{\hat{m}_{b} } {\tilde{b}}}} \right. \kern-0pt} {\tilde{b}}} + \hat{s}} \right)^{{\hat{a}}},$$
(43)

respectively, assuming that the Hoek–Brown criterion governs the stress ratio since the onset of plastic deformation and the residual state. \(\sigma_{ci}\) is the uniaxial compressive strength of intact rock. The parameters \(m_{b}\), \(s\) and \(a\) are constants which depend upon the rock mass characteristics: the Geological Strength Index (GSI) and disturbance factor (D). The circumflex letters, \(\hat{m}_{b}\), \(\hat{s}\) and \(\hat{a}\), in Eq. (43) represent the Hoek–Brown parameters at the residual state. Then, the corresponding non-coaxial stress–dilatancy equations may be established by considering plastic deformations under loading and under and unloading as follows.

4.2.1 Loading

Assuming loading and substituting Eqs. (42) and (43) into Eq. (14) leads to

$$m_{s} N_{\psi{\mathrm{-peak}}}^{\mathrm{HB},L} = \frac{{N_{\sigma{\mathrm{-peak}}}^{\mathrm{HB}} + \bar{K}_{N}^{\mathrm{HB}} \tan^{2} \Delta }}{{\bar{K}_{N}^{\mathrm{HB}} + N_{\sigma{\mathrm{-peak}}}^{\mathrm{HB}} \tan^{2} \Delta }}.$$
(44)

Substituting Eqs. (42) and (43) into Eq. (44), the dilatancy angle may be obtained as

$$\sin \hat{\psi }_{{\mathrm{peak}}}^{\mathrm{HB},L} \, = - \frac{{m_{s} N_{\psi{\mathrm{-peak}}}^{\mathrm{HB},L} - 1}}{{m_{s} N_{\psi{\mathrm{-peak}}}^{\mathrm{HB},L} + 1}} = - c_{\Delta } \frac{{\varGamma - \tilde{\varGamma }}}{{2 + \varGamma + \tilde{\varGamma }}}.$$
(45)

where \(\varGamma\) and \(\tilde{\varGamma }\) are, respectively,

$$\varGamma = \tilde{b}{{(m_{b} } \mathord{\left/ {\vphantom {{(m_{b} } {\tilde{b}}}} \right. \kern-0pt} {\tilde{b}}} + s)^{a} \quad {\text{and}}\quad \tilde{\varGamma } = f_{sd} \tilde{b}{{(\hat{m}_{b} } \mathord{\left/ {\vphantom {{(\hat{m}_{b} } {\tilde{b}}}} \right. \kern-0pt} {\tilde{b}}} + \hat{s})^{{\hat{a}}} .$$
(46)

4.2.2 Unloading

For unloading, considering \(C_{N} = 1/\bar{K}_{N}^{\mathrm{HB}}\) such that

$$m_{s} N_{\psi{\mathrm{-peak}}}^{\mathrm{HB},U} = \frac{{\bar{K}_{N}^{\mathrm{HB}} N_{\sigma{\mathrm{-peak}}}^{\mathrm{HB}} + \tan^{2} \Delta }}{{1 + N_{\sigma{\mathrm{-peak}}}^{\mathrm{HB}} \bar{K}_{N}^{\mathrm{HB}} \tan^{2} \Delta }},$$
(47)

substituting Eqs. (42) and (43) into Eq. (44) and considering the definition in Eq. (36), one obtains

$$\sin \hat{\psi }_{{\mathrm{peak}}}^{\mathrm{HB},U} \, = - c_{\Delta } \frac{{\varGamma + \tilde{\varGamma } + \varGamma \tilde{\varGamma }}}{{2 + \varGamma + \tilde{\varGamma } + \varGamma \tilde{\varGamma }}}.$$
(48)

4.3 Non-coaxial stress–dilatancy formulation for Lode angle-dependent yield functions

In this section, we are going to demonstrate how specific stress–dilatancy relations and plastic dissipation equations might be established for the general stress–strain condition following the formalism established in Sect. 2.2 when the yield functions are specified.

Let the stress ratio be written in terms of the triaxial compression stress ratio and a convenient Lode angle-dependent function as

$$M_{\sigma }^{\theta } = \ell^{\theta } M_{\sigma }^{C} ,$$
(49)

where \(M_{\sigma }^{C}\) is the stress ratio for triaxial compression and \(\ell^{\theta }\) is the Lode angle-dependent modification.

Accordingly, the modified ‘constant’ may be obtained as a stress ratio at the phase transformation condition as

$$\hat{C}_{M}^{\theta } = \hat{\bar{K}}_{M}^{\theta } = \hat{\bar{\ell }}_{c}^{\theta } \hat{\bar{K}}_{M}^{\theta } ,$$
(50)

where \(\hat{\bar{\ell }}_{c}^{\theta }\) is a modified Lode angle-dependent function and \(\hat{\bar{K}}_{M}^{\theta }\) is the modified phase transformation stress ratio which evolves to the critical state stress ratio as the plastic deformation progresses toward the critical state condition.

Assuming the Mohr–Coulomb yield function, the current stress ratio and the stress ratio at the phase transformation for the triaxial compression and the triaxial extension states are given, respectively, as

$$M_{c}^{{\mathrm{MC}},C/E} = M_{\varphi }^{C/E} (1 + \hat{b})\quad {\text{and}}\quad \bar{K}_{M}^{C/E} = \bar{M}_{c}^{{\mathrm{MC}},C/E} (1 + \hat{b}) ,$$
(51)

wherein

$$M_{\varphi }^{C/E} = 3\frac{{N_{\varphi } - 1}}{{r_{1} N_{\varphi } + r_{3} }},\bar{M}_{c}^{C/E} = 3\frac{{\bar{N}_{c} - 1}}{{r_{1} \bar{N}_{c} + r_{3} }}\,{\text{and}}\,\hat{b} = a/p.$$
(52)

\(N_{\varphi }\) and \(\bar{N}_{c}\) are as defined in Eq. (34), \(r_{1}\) and \(r_{3}\) are as defined in Eq. (5), and \(a\) is attraction. The subscript c indicates critical, and the superscript C/E indicates triaxial compression or extension conditions. The stress ratio is then obtained by multiplying the triaxial compression stress ratio by an appropriate Lode angle-dependent function. The Bardet [4] Lode angle-dependent function may be considered for instance.

Combining Eqs. (27), (49) and (50), we obtain

$$M_{\psi }^{\theta } = - \tilde{c}(s\ell^{\theta } M_{\sigma }^{C} - \bar{\ell }_{c}^{\theta } \bar{M}_{\sigma }^{C} ).$$
(53)

Inserting Eqs. (51) into Eq. (53), we obtain a stress–dilatancy relationship as

$$M_{\psi }^{\theta } = - \tilde{c}(s\ell^{\theta } M_{\sigma }^{C} - \bar{\ell }_{c}^{\theta } \bar{M}_{c}^{C} )(1 + \hat{b}).$$
(54)

5 Preliminary validation of the proposed  theory

5.1 Dependence of plastic dissipation on non-coaxiality

Let us consider the experimental investigations due to Gutierrez and Ishihara [12] and their implication interms of effects of non-coaxiality. The original data from Gutierrez and Ishihara [12], Fig. 7 (Left), show a plot of the dilatancy ratio against the stress ratio multiplied by the degree of coaxiality, \(\tilde{c} = |\tilde{m}_{ij} \tilde{n}_{ij} |\), for monotonic tests and for tests subjected to pure principal stress rotation. The following can be inferred from the plot.

Fig. 7
figure 7

Left: Stress ratio (q/p) multiplied by degree of coaxiality, \(\tilde{c}\), versus dilatancy ratio (\(- {{M}}_{{{\psi}}}\)) for tests with pure rotation of principal stresses and with monotonic loading, data from Gutierrez and Ishihara [12], Right: Normalized plastic dissipation versus dilatancy ratio interpreted from the data

Full size image

  • At lower mobilizations of the stress ratio, the stress–dilatancy relation is less contractive for the tests with pure stress rotation than those for monotonic loading.

  • On average, the plots tend to have a unique stress ratio at the phase transformation state.

  • On average, the differences in the dilatancy ratio between the two data sets vanish with stress ratio

Further, the data are interpreted in terms of the normalized plastic dissipation, Fig. 7 (right). From the trend of the plastic dissipation for the two sets of data, it can be observed that

Fig. 8
figure 8

Predicted versus measured degrees of coaxiality (\(\cos 2\Delta\)). The predicted degree of coaxiality (\(\cos 2\Delta\)) is calculated using the measured dilatancy ratio as an input into Eq. (55)

Full size image

  • the tests with pure rotation dissipate less than those from monotonic loading.

  • on average the difference in the normalized plastic dissipation between the two sets of data vanishes with stress ratio.

It is known that tests with principal stress rotation induce non-coaxiality between the directions of principal stresses and principal plastic strains increments. The data is therefore in agreement with the non-coaxial plastic dissipation and stress–dilatancy theory presented in this paper.

5.2 Evolution of the degree of non-coaxiality

Equation (30) can be analytically integrated such that it yields the relationship

$$\tan \Delta = \sqrt {\frac{{\chi_{\psi } }}{{\chi_{\psi } + (1 + m_{s} N_{\psi } )^{{\pi^{{a_{\Delta } }} }} }}} ;\quad \chi_{\psi } = \frac{{\tan^{2} \Delta_{0} }}{{1 - \tan^{2} \Delta_{0} }}\left( {1 + m_{s} N_{\psi ,0} } \right)^{{\pi^{{a_{\Delta } }} }} .$$
(55)

The dilatancy ratio calculated from Vardoulakis and Georgopoulos [34], in Fig. 3 is then fed into Eq. (55). A good agreement is obtained between the trend of the experimentally measured degree of coaxiality, \(\cos 2\Delta ,\) and that of the trend of the degree of coaxiality calculated using Eq. (55), Fig. 8.

Integrating Eq. (31) yields the relationship

$$\tan \Delta = \sqrt {\frac{{\chi_{\sigma } }}{{\chi_{\sigma } + (C_{{\tilde{N}}} + N_{\sigma } )^{{\pi^{{a_{\Delta } }} }} }}} ;\quad \chi_{\sigma } = \frac{{\tan^{2} \Delta_{0} }}{{1 - \tan^{2} \Delta_{0} }}\left( {C_{{\tilde{N}}} + N_{\sigma 0} } \right)^{{\pi^{{a_{\Delta } }} }} .$$
(56)

The trend compares reasonably well with data from DSC tests in Arthur et al. [3], Fig. 9.

Fig. 9
figure 9

The tendencey of the degree of non-coaxiality with increasing stress ratio, data from DSC tests in Arthur et al. [3] in comparison with theory, Eq. (56)

Full size image

6 Conclusion

Inconsistencies were noted in some of the theoretical frameworks proposed in literature for the description of stress-dilatancy relations and plastic dissipations in geomaterials when principal stresses and principal plastic strain rates were non-coaxial, which added to the motivation of the work presented in this paper. The paper deals with the development of, a consistent and unifying theoretical framework that describes non-coaxial plastic dissipation and stress–dilatancy relations for geomaterials. For a possible use in models that intend the modelling of deformation behaviour of geomaterials under cyclic loading, both loading and unloading were explicitly considered. The framework accomodates and extends the well-known stress-dilatancy theories that have so far assumed a central place in constitutive modelling of geomaterials. Selected experimental results from literature that demonstrated the validity of elements of the theoretical framework have also been presented. Through theoretical arguments and looking at published experimental results, it is concluded that

  • plastic dissipation decreases with increasing of degree of non-coaxiality between axes of principal stresses and principal plastic strain increments.

  • when the degree of non-coaxiality tends to increase, the volumetric changes due to dilatancy tend to decrease.

  • the degree of coaxiality increases with increasing dilatancy ratio and increasing stress ratio prior to bifurcation.

In addition, limits of the degree of non-coaxiality have been obtained for the case of a plane strain deformation requiring that the plastic dissipation be non-negative. The limits of the degree of non-coaxiality so obtained were found to be more relaxed than previously proposed and they accommodated the maximum values of degrees of non-coaxiality previously observed in DEM simulations. Furthermore, the formulation presented in this paper for describing the evolution of the degree of coaxiality conforms to empirical observations in the pre-bifurcation region. After bifurcation, some inconsistencies can be explained by the proposed framework. However, available experimental evidences are insufficient to draw conclusions.