Introduction

We study stress–strain rate constitutive relations describing granular materials like cohesionless soils or broken rocks within the hypoplastic theory proposed by Kolymbas [21] and developed further by [34, 36, 37]. Unlike hyper- and hypoelastic material laws, the hypoplastic response differs for loading and unloading, thus corresponding to inelastic materials [19]. On the other hand, the strain is not decomposed into elastic and plastic parts unlike classical elastoplastic concepts [3, 20]. The hypoplastic model is incrementally nonlinear and of the rate type. See other incrementally nonlinear constitutive equations in [4, 13], the rate-independence concept and hysteresis in [31, 32, 35]. Under cyclic loading, one can naturally obtain limit circles with barodesy [22, 23]. For mathematical modeling of granular and multiphase media, we cite [15,16,17, 28,29,30], and for well-posedness and stability analysis of the relevant differential equations, we refer to [11, 12, 14, 33].

Our study relies on a simplified version of the hypoplastic constitutive relation that was introduced independently by Bauer [5] and Gudehus [18]. In previous works [8, 9, 24, 26, 27], we omitted the pressure and density dependence of soil materials and considered the strain–stress law as a nonlinear differential equation for the stress under a given proportional strain rate, the so-called strain control. The existence of an analytical solution made it possible to describe explicitly stress paths obtained for the specific cases of monotonic compression, extension, and isochoric deformations. Recently, the case of unknown strain rate that should be derived from a given proportional stress, called stress control, was recovered within implicit differential equations in [10, 25] together with the pressure and density dependence of model parameters in the hypoplastic law. In such studies, the density of a granular material, which represents the distribution of pores and solid grains, is expressed with the help of the corresponding void ratio.

Our current mathematical modeling is based on the degradation of granular hardness suggested by Bauer [6, 7], which is suitable for the description of mechanical weathering of soils and granular materials under dry and wet states. This can be caused by long-time phenomena in various environmental events. We investigate the dynamic of the hypoplastic model subjected to the exponential degeneration for special cases of strain and stress control. Under the strain control, the spherical stress is constant. Whereas in the stress control case, the spherical components of stress and strain rate are described by implicit differential equations, which yield to two different solutions, and consequently, two degradation scenarios. It is shown that both solutions converge asymptotically to an attractor implying sparsification of the material state. When cyclic loading-unloading was applied, then finite ratcheting occurred, that is, a shift of the hysteresis loops was observed in a numerical simulation carried out by the standard 4th order Runge–Kutta method. The corresponding hysteresis curve of the deviatoric strain rate factor has, in this case, the form of a square spiral.

Modeling

For the reader’s convenience, the notations used throughout the paper are collected alphabetically in Table 1.

Table 1 Alphabetical list of notations
Full size table

In this study, we consider a deformable granular body with constant volume and whose granular hardness, \(h_\textrm{s}\), can be subjected to degradation according to the following phenomenological equation found in [7]:

$$\begin{aligned} \dot{h}_\textrm{s} =\frac{1}{c} \Bigl ( h_\textrm{s}^\infty -h_\textrm{s}\Bigr ),\quad c\ge 0, \end{aligned}$$
(1)

where the parameter c has the dimension of time and \(h_\textrm{s}^\infty\) denotes the final value of the degraded solid hardness. Therefore, the dependence of the granular hardness on the time t is expressed for prescribed \(0\le h_\textrm{s}^\infty <h_\textrm{s}^0\) by an exponential function

$$\begin{aligned} h_\textrm{s}(t) =h_\textrm{s}^\infty +\bigl ( h_\textrm{s}^0 -h_\textrm{s}^\infty \bigr ) \exp \left( \frac{t_0-t}{c}\right) , \end{aligned}$$
(2)

where \(h_\textrm{s}^0 =h_\textrm{s}(t_0)\) for some \(t_0\ge 0\), and \(h_\textrm{s}(t)\rightarrow h_\textrm{s}^\infty\) as \(t\rightarrow \infty\). In particular, if \(h_\textrm{s}^\infty =0\), then Eq. 1 follows a constant proportion

$$\begin{aligned} \frac{\dot{h}_\textrm{s}}{h_\textrm{s}} =-\frac{1}{c}. \end{aligned}$$

It is worth mentioning that the degradation of the solid hardness as described in Eq. 1 is reasonable for modeling the cases when the material is subjected, for example, to variable weather conditions. Indeed, in these cases, even in the absence of applied external forces, we may still observe a degradation of the material caused by progressive or instantaneous weathering; see [7].

To describe the constitutive stress–strain relation for hypoplastic granular materials, we consider an extension of the model from [10] that takes into account the degradation of the solid hardness. More precisely, considering the Cauchy stress \(\varvec{\sigma }\) to be a three-by-three matrix and assuming that its orthogonal decomposition after normalization is of the form

$$\begin{aligned} \hat{\varvec{\sigma }} =\frac{\varvec{\sigma }}{\textrm{tr}\varvec{\sigma }} =\hat{\varvec{\sigma }}^* +\frac{1}{3} \textbf{I}, \end{aligned}$$
(3)

the constitutive response between the stress and the three-by-three matrix of the strain rate \(\dot{\varvec{\varepsilon }}\) can be expressed as follows

$$\begin{aligned} \dot{\varvec{\sigma }} =f_\textrm{s} \textrm{tr}\varvec{\sigma } \left\{ a^2 \dot{\varvec{\varepsilon }} +(\dot{\varvec{\varepsilon }}: \hat{\varvec{\sigma }}) \hat{\varvec{\sigma }} +a f_\textrm{d} (\hat{\varvec{\sigma }} +\hat{\varvec{\sigma }}^*) \Vert \dot{\varvec{\varepsilon }}\Vert \right\} +\textrm{tr}\varvec{\sigma } \frac{\dot{h}_\textrm{s}}{h_\textrm{s}} \left( \kappa \hat{\varvec{\sigma }}^* +\frac{1}{3} \textbf{I} \right) , \end{aligned}$$
(4)

using the dyadic product of tensors, where \(a>0\) is the yield strength, weight \(\kappa \in (0,1)\), \(\textbf{I}\) is the three-by-three identity matrix, and \(\Vert \dot{\varvec{\varepsilon }}\Vert =\sqrt{\dot{\varvec{\varepsilon }}:\dot{\varvec{\varepsilon }}}\) stands for the Frobenius norm. Concerning the physical interpretation of the right-hand side of the constitutive relation Eq. 4, the two terms which are linear with respect to the strain rate \(\dot{\varvec{\varepsilon }}\) account for the hypoelastic behavior, whereas the term with the norm \(\Vert \dot{\varvec{\varepsilon }}\Vert\) describes the inelastic behavior and allows for different signs of the strain rate, therefore distinguishing the loading and unloading cycles. Note that the inelastic term is weighted by the density factor \(f_\textrm{d}>0\), while the stiffness factor, \(f_\textrm{s}<0\), influences both the hypoelastic and the inelastic processes of the system. Finally, a non-constant granular hardness \(h_\textrm{s}\) entering Eq. 4 allows degradation phenomena.

For simplicity, the governing relations are given in time \(t\ge t_0\ge 0\) and we omit the spatial dependence for all physical variables. Note that, unlike the model in [7], the constitutive relation Eq. 4 is rate dependent. This is due to the perturbation originated by the term accounting for the degradation which, as seen in Eq. 2, shows a strong dependence on the time. There is a number of studies devoted to the investigation of rate-dependent processes, and in the past few years, this theory became more popular in the engineering community for applications in modeling hysteresis in multi-functional materials, see, e.g., [1, 2].

For cohesionless granular materials only non-positive principal stresses and not all equal to zero

$$\begin{aligned} \sigma _1 \le 0,\quad \sigma _2 \le 0,\quad \sigma _3 \le 0 \end{aligned}$$
(5)

are physically relevant; this yields, in particular, that the pressure \(p =-(\sigma _1 +\sigma _2 +\sigma _3)/3\) is positive. Like in [7], the factors \(f_\textrm{d}\) and \(f_\textrm{s}\) depend on the void ratio e. Indeed, in our investigation, while the volume of the grains in the body is assumed to be a constant \(V_\textrm{s}\), the void space has variable volume \(V_\textrm{v}\), hence the void ratio e is considered as a state quantity

$$\begin{aligned} e =\frac{V_\textrm{v}}{V_\textrm{s}}\in (e_\textrm{d}, e_\textrm{i}). \end{aligned}$$
(6)

The bounds \(e_\textrm{d}>0\) and \(e_\textrm{i}>0\) depend proportionally on the mechanical pressure according to the following identities found in [6] (which are related to Eq. 15)

$$\begin{aligned} e_\textrm{d} =e_\textrm{min} \exp (- y),\quad e_\textrm{i} =e_\textrm{max} \exp (- y), \end{aligned}$$
(7)

for fixed values \(0< e_\textrm{min} <e_\textrm{max}\), by means of the auxiliary variable

$$\begin{aligned} y =\left( \frac{3 p}{h_\textrm{s}}\right) ^n=\left( \frac{- \textrm{tr}\varvec{\sigma }}{h_\textrm{s}}\right) ^n. \end{aligned}$$
(8)

Notably, the variable y appears in the compression law proposed by Bauer with typically \(n\in (0,1)\) and is related to the point of inflection observed in experiments in the compression curve, see [7].

Now, observe that \(\textrm{tr} \dot{\varvec{\varepsilon }}\) is the relative volume increment and, as mentioned before, the solid volume \(V_\textrm{s}\) is constant during the evolution. Denoting by V the total volume, the physical consistency of the system relies on a volume balance equation which reads as follows:

$$\begin{aligned} \dot{V_\textrm{v}}=V \textrm{tr} \dot{\varvec{\varepsilon }}= (V_\textrm{v}+V_\textrm{s}) \textrm{tr} \dot{\varvec{\varepsilon }} \end{aligned}$$

which in terms of the void ratio e introduced in Eq. 6 can be interpreted as

$$\begin{aligned} \dot{e} =(1 +e) \dot{\delta }, \end{aligned}$$
(9)

where \(\delta =\textrm{tr}\varvec{\varepsilon }\) is the volume strain (dilatation). Further, the solution of Eq. 9 for prescribed \(\delta _0\) and \(e_0\) is

$$\begin{aligned} \delta -\delta _0 =\ln \frac{1 +e}{1 +e_0}. \end{aligned}$$
(10)

Note that Eq. 9 can be also understood as a mass balance equation. Indeed, recalling that the material’s density is expressed in terms of the mass m of the granular body by

$$\begin{aligned} \rho =\frac{m}{V_\textrm{s}(1 +e)} \end{aligned}$$

we get the continuity equation

$$\begin{aligned} \dot{\rho } +\rho \, \dot{\delta } =0. \end{aligned}$$

In order to obtain a more approachable formulation of the problem, let us consider the decomposition of the strain rate into deviatoric and spherical parts

$$\begin{aligned} \dot{\varvec{\varepsilon }} =\dot{\varvec{\varepsilon }}^* +\frac{1}{3} \dot{\delta } \textbf{I},\quad \dot{\delta } =\textrm{tr}\dot{\varvec{\varepsilon }}. \end{aligned}$$
(11)

By differentiating Eq. 3, we get

$$\begin{aligned} \dot{\hat{\varvec{\sigma }}}^* =\frac{\dot{\varvec{\sigma }}}{\textrm{tr}\varvec{\sigma }} -\frac{\textrm{tr}\dot{\varvec{\sigma }}}{\textrm{tr}\varvec{\sigma }} \hat{\varvec{\sigma }}. \end{aligned}$$

Recalling that \(p =-\textrm{tr}\varvec{\sigma }/3\), from Eq. 4, it follows the following scalar differential equation for the trace

$$\begin{aligned} \frac{\dot{p}}{p}=\frac{\textrm{tr}\dot{\varvec{\sigma }}}{\textrm{tr}\varvec{\sigma }} =f_\textrm{s} \left\{ a^2 \dot{\delta } +\dot{\varvec{\varepsilon }}: \hat{\varvec{\sigma }} +a f_\textrm{d} \Vert \dot{\varvec{\varepsilon }}\Vert \right\} +\frac{\dot{h}_\textrm{s}}{h_\textrm{s}}. \end{aligned}$$

Using the identities above, we rewrite Eq. 4 in terms of the deviatoric part

$$\begin{aligned} \dot{\hat{\varvec{\sigma }}}^* =f_\textrm{s} \left\{ a^2 (\dot{\varvec{\varepsilon }}^* -\dot{\delta } \hat{\varvec{\sigma }}^*) +a f_\textrm{d} \hat{\varvec{\sigma }}^* \Vert \dot{\varvec{\varepsilon }}\Vert \right\} -(1-\kappa ) \frac{\dot{h}_\textrm{s}}{h_\textrm{s}} \hat{\varvec{\sigma }}^*, \end{aligned}$$
(12)

while by differentiating Eq. 8, we obtain

$$\begin{aligned} \frac{\dot{y}}{n y} =\frac{\dot{p}}{p} -\frac{\dot{h}_\textrm{s}}{h_\textrm{s}} =f_\textrm{s} \left\{ \left( a^2 +\frac{1}{3} \right) \dot{\delta } +\dot{\varvec{\varepsilon }}^*: \hat{\varvec{\sigma }}^* +a f_\textrm{d} \Vert \dot{\varvec{\varepsilon }}\Vert \right\} . \end{aligned}$$
(13)

Besides accounting for the effects of degradation, in our study of the constitutive equation for stress–strain relation, we assume that the void ratio changes in time according to the equation

$$\begin{aligned} \dot{e} =-e n y \left( \frac{\dot{p}}{p} -\frac{\dot{h}_\textrm{s}}{h_\textrm{s}} \right) =-e \dot{y}. \end{aligned}$$
(14)

The identity above has been proposed in [7] for modeling under isotropic compression and it allows us to simplify the system by restoring constant density and stiffness factors in Eq. 13. Similar technique has been used in [25] for studying the long-term dynamic of such type of system. More precisely, for prescribed \(y_0\ge 0\) and \(e_0\in (e_\textrm{min}, e_\textrm{max}) \exp (- y_0)\), the positive solution to Eq. 14 is expressed in terms of the auxiliary variable y as

$$\begin{aligned} e =e_0\exp (y_0-y). \end{aligned}$$
(15)

In turn, the density and stiffness factors are functions of the current void ratio

$$\begin{aligned} f_\textrm{d}(e) =\left( \frac{e -e_\textrm{d}}{e_\textrm{c} -e_\textrm{d}} \right) ^\alpha , \quad f_\textrm{s}(e) =-b\left( \frac{e_\textrm{i}}{e} \right) ^\beta ,\quad \alpha \in (0,0.5),\quad \beta >1, \end{aligned}$$
(16)

for a model parameter \(b >0\), and the critical void ratio \(e_c\) is described, as counterpart of Eq. 7, by the identity

$$\begin{aligned} e_\textrm{c} =e_\textrm{crt} \exp (- y) \end{aligned}$$
(17)

with constant \(e_\textrm{crt}\in (e_\textrm{min}, e_\textrm{max})\). As a consequence of Eqs. 15 and 16 the dependence of the void ratio can be suppressed and, thanks to Eq. 7, both the density and the stiffness correspond to the following constant factors

$$\begin{aligned} f_\textrm{d} =\left( \frac{\exp (y_0) -e_\textrm{min}}{e_\textrm{crt} -e_\textrm{min}} \right) ^\alpha ,\quad f_\textrm{s} =-b \left( \frac{e_\textrm{max}}{e_0 \exp (y_0)} \right) ^\beta . \end{aligned}$$
(18)

In summary, the coupled system we want to investigate comprehends the algebraic Eq. 8, and the differential Eqs. 12 and 13, for unknowns \(\varvec{\varepsilon }^*, p, \hat{\varvec{\sigma }}^*\). Our investigation of system Eqs. 12 and 13 focuses on its behavior when subjected to periodically oscillating loading-unloading cycles, therefore its analysis needs to be carried out in the corresponding monotonicity intervals of the cycle. Since the system is positive 1-homogeneous, by choosing a loading parameter s(t) to be a monotone differentiable function in some interval \(t\in (t_0,t_1)\), assuming that the prime denotes the derivative with respect to s, from Eq. 12, we get by the chain rule

$$\begin{aligned} (\hat{\varvec{\sigma }}^*)^\prime \dot{s} =f_\textrm{s} \left\{ a^2 \bigl ( (\varvec{\varepsilon }^*)^\prime \dot{s}-\delta ^\prime \dot{s} \hat{\varvec{\sigma }}^* \bigr ) + a f_\textrm{d} \hat{\varvec{\sigma }}^* \Vert \varvec{\varepsilon }^\prime \dot{s}\Vert \right\} -(1-\kappa ) \frac{h^\prime _\textrm{s}\dot{s}}{h_\textrm{s}} \hat{\varvec{\sigma }}^*, \end{aligned}$$

where \(\dot{s}\) stands for the derivative of s with respect to t; and analogously for Eq. 13. Thus, dividing by this term, we finally obtain the corresponding system of differential Eqs. 12 and 13 in terms of s. More precisely:

$$\begin{aligned} (\hat{\varvec{\sigma }}^*)^\prime =f_\textrm{s} \left\{ a^2 \bigl ( (\varvec{\varepsilon }^*)^\prime -\delta ^\prime \hat{\varvec{\sigma }}^* \bigr ) +s_\pm a f_\textrm{d} \hat{\varvec{\sigma }}^* \Vert \varvec{\varepsilon }^\prime \Vert \right\} -(1-\kappa ) \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \hat{\varvec{\sigma }}^*, \end{aligned}$$
(19)
$$\begin{aligned} \frac{y^\prime }{n y} =f_\textrm{s} \left\{ \left( a^2 +\frac{1}{3} \right) \delta ^\prime +(\varvec{\varepsilon }^*)^\prime : \hat{\varvec{\sigma }}^* +s_\pm a f_\textrm{d} \Vert \varvec{\varepsilon }^\prime \Vert \right\} , \end{aligned}$$
(20)

where \(s_\pm =1\) if \(\dot{s}>0\), and \(s_\pm =-1\) if \(\dot{s}<0\). With this in mind, we identify a loading cycle for s increasing, i.e., \(\dot{s}>0\), while an unloading cycle refers to s decreasing, \(\dot{s}<0\).

The governing system comprehends two equations: one tensor equation for the deviator, Eq. 19, and one scalar equation for the trace, Eq. 20, which is implicitly given in terms of the auxiliary variable y introduced in Eq. 8, where the granular hardness \(h_\textrm{s}\) represents an inhomogeneous coefficient given by Eq. 2. Note that both equations are implicit through the norm \(\Vert \varvec{\varepsilon }^\prime \Vert\). Besides, the system is under-determined for unknown \(p, (\varvec{\varepsilon }^*)^\prime , \hat{\varvec{\sigma }}^*\) and it can be interpreted in terms of the stress if the strain is given, and vice versa. More precisely, the strain-controlled case consists in choosing a fixed tensor \(\textbf{D}\) normalized by \(\textrm{tr} \textbf{D} =1\), and assuming that the strain tensor has the form \(\varvec{\varepsilon } =\delta (s(t)) \textbf{D}\) as presented in the “Strain control” section. Analogously, in the stress-controlled case, we assume that a fixed tensor \(\textbf{T}\) is normalized by \(\textrm{tr} \textbf{T} =1\), and consider the stress given by \(\varvec{\sigma } =-3 p(s(t)) \textbf{T}\) in the “Stress control” section. The study of the long-time behavior for periodical cycling loading-unloading under stress control is contained in the “Long-time behavior of the stress-controlled case” section.

As mentioned before, in our investigation, we will also take into account equations Eqs. 9 and 14 which describe the changes of the void ratio. Combining these two identities, together with Eq. 15, and considering the time transformation, we get

$$\begin{aligned} \delta ^\prime =-\frac{e}{1+e} y^\prime =-\frac{y^\prime }{1+\frac{\exp (y-y_0)}{e_0}}. \end{aligned}$$
(21)

We can see that the monotonicity intervals of y and \(\delta\) coincide; therefore, it is natural to choose the loading parameter s so that \(\dot{s}>0\) for \(\dot{\delta }>0\), while \(\dot{s}<0\) for \(\dot{\delta }<0\). The opposite relation holds for the variable y, that is, \(\dot{s}>0\) for \(\dot{y}<0\), and \(\dot{s}<0\) for \(\dot{y}>0\). Since y is associated with the pressure \(p=-\textrm{tr} \varvec{\sigma }\) via identity Eq. 8, the rate \(\dot{y}<0\) corresponds to compressive loading, and \(\dot{y}>0\) to unloading related to the extension phase.

Throughout the paper, we consider the initial conditions at \(t=t_0\) corresponding to the loading parameter value \(s(t_0)\) prescribed by given \(e_0, p_0, \delta _0, (\varvec{\varepsilon }^*)^\prime _0, \hat{\varvec{\sigma }}^*_0\) and \(y_0 =(3 p_0/h_\textrm{s}^0)^n\).

Analysis

The solvability of the initial value problem in the general case is open. We study the special cases under strain and stress control conditions.

Strain control

For a prescribed three-by-three symmetric second-order tensor \(\textbf{D}\) such that \(\textrm{tr} \textbf{D} =1\), let us assume that

$$\begin{aligned} \varvec{\varepsilon }^\prime =\delta ^\prime \textbf{D},\quad \textbf{D} =\textbf{D}^* +\frac{1}{3} \textbf{I}, \end{aligned}$$
(22)

hence its deviatoric part reads \((\varvec{\varepsilon }^*)^\prime =\delta '\textbf{D}^*\), and

$$\begin{aligned} \Vert \varvec{\varepsilon }^\prime \Vert =|\delta ^\prime | \sqrt{\Vert \textbf{D}^*\Vert ^2 +\frac{1}{3}}. \end{aligned}$$

Recalling that the loading parameter s follows the monotonicity intervals of \(\delta\), from the identities above and Eq. 19, we get the following differential equation with respect to \(\hat{\varvec{\sigma }}^*\) coupled with \(\delta ^\prime\):

$$\begin{aligned} (\hat{\varvec{\sigma }}^*)^\prime =f_\textrm{s} \left\{ a^2 (\textbf{D}^* -\hat{\varvec{\sigma }}^*) +s_\pm a f_\textrm{d} \hat{\varvec{\sigma }}^* \sqrt{\Vert \textbf{D}^*\Vert ^2 +\frac{1}{3}} \right\} \delta ^\prime -(1-\kappa ) \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \hat{\varvec{\sigma }}^*. \end{aligned}$$
(23)

Similarly, from Eq. 20, we obtain

$$\frac{y^\prime }{ny} =f_\textrm{s}\left\{ a^2 +\frac{1}{3} +\textbf{D}^*: \hat{\varvec{\sigma }}^* +s_\pm a f_\textrm{d} \sqrt{\Vert \textbf{D}^*\Vert ^2 +\frac{1}{3}}\right\} \delta ^\prime ,$$

this together with Eq. 21 implies

$$\begin{aligned} \left( \frac{1}{f_\textrm{s} n y} +\frac{\lambda _{\pm } +\textbf{D}^*: \hat{\varvec{\sigma }}^*}{1+\frac{\exp (y-y_0)}{e_0}} \right) y^\prime =0,\quad \lambda _{\pm } =a^2 +\frac{1}{3} +s_\pm a f_\textrm{d} \sqrt{\Vert \textbf{D}^*\Vert ^2 +\frac{1}{3}}. \end{aligned}$$
(24)

Equation 24 is trivially satisfied if \(y^\prime =0\), that is, \(y=y_0\) is constant and

$$\begin{aligned} e=e_0,\quad \delta =\delta _0,\quad p =\frac{p_0}{h_\textrm{s}^0} h_\textrm{s}. \end{aligned}$$
(25)

In this case, by Eq. 21, we have \(\delta ^\prime =0\) and consequently the differential Eq. 23 for the deviatoric stress corresponds to

$$\begin{aligned} (\hat{\varvec{\sigma }}^*)^\prime =-(1-\kappa ) \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \hat{\varvec{\sigma }}^*, \end{aligned}$$

hence its solution can be obtained explicitly and has the form

$$\begin{aligned} \hat{\varvec{\sigma }}^* =\hat{\varvec{\sigma }}^*_0 \left( \frac{h_\textrm{s}}{h_\textrm{s}^0} \right) ^{-(1-\kappa )}. \end{aligned}$$

Note that, for \(\delta '=y'=0\), Eq. 22 indicates that there is no change of strain; therefore, in accordance with the identity above, the stress forces are due only to the degradation of the granular material.

The other possible solution of Eq. 24 is expressed by the transcendental equation for y coupled with \(\hat{\varvec{\sigma }}^*\), namely,

$$\begin{aligned} \frac{1}{f_\textrm{s} n y} +\frac{\lambda _{\pm } +\textbf{D}^*: \hat{\varvec{\sigma }}^*}{1+\frac{\exp (y-y_0)}{e_0}} =0. \end{aligned}$$
(26)

Although the solvability of Eq. 26 follows from the general theory of ordinary differential equations, there is no explicit formula for its solution. The investigation of such a problem can thus rely on numerical techniques and it is a question for a future work. The theoretical and numerical analysis of a related hypoplastic problem like Eq. 4 for the strain-controlled case has been carried out in [8, 9, 24, 26, 27], however, disregarding the degradation effects.

Stress control

For a prescribed three-by-three symmetric second-order tensor \(\textbf{T}\) with non-positive eigenvalues, and \(\textrm{tr} \textbf{T} =1\), let

$$\begin{aligned} \frac{\varvec{\sigma }}{\textrm{tr}\varvec{\sigma }} =\hat{\varvec{\sigma }} =\textbf{T} =\textbf{T}^* +\frac{1}{3} \textbf{I}, \end{aligned}$$
(27)

such that Eq. 5 holds and \(\textrm{tr} \varvec{\sigma }<0\). Inserting Eq. 27 into Eqs. 19 and 20, dividing by \(f_\textrm{s}\) and \(a^2\), using \((\hat{\varvec{\sigma }}^*)^\prime =\textbf{0}\), we get the coupled system for \((\varvec{\varepsilon }^*)^\prime\) and \(\delta ^\prime\):

$$\begin{aligned} (\varvec{\varepsilon }^*)^\prime =k_\textrm{d} \textbf{T}^*,\quad k_\textrm{d} =\delta ^\prime +\frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} -s_\pm \frac{f_\textrm{d}}{a} \Vert \varvec{\varepsilon }^\prime \Vert , \end{aligned}$$
(28)
$$\begin{aligned} \frac{y^\prime }{f_\textrm{s} n y} =\left( a^2 +\frac{1}{3} \right) \delta ^\prime +(\varvec{\varepsilon }^*)^\prime : \textbf{T}^* +s_\pm a f_\textrm{d} \Vert \varvec{\varepsilon }^\prime \Vert . \end{aligned}$$
(29)

Note that, the first equation is implicit as the norm \(\Vert \varvec{\varepsilon }^\prime \Vert =\sqrt{\Vert (\varvec{\varepsilon }^*)^\prime \Vert ^2 +(\delta ^\prime )^2/3}\). On the other hand, recalling the definition of y, we can see that Eq. 29 gives an equation for the pressure p, while the assumption of Eq. 27 indicates the relation between the stress components and the pressure, \(\varvec{\sigma } =-3 p \textbf{T}\).

We will show that, under certain conditions, we can find an explicit solution for the system. The first step towards such an analytical solution is to derive an expression for the norm \(\Vert \varvec{\varepsilon }^\prime \Vert\).

Theorem 1

Let Eq. 27 hold and assume that

$$\begin{aligned} 1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \ge 0. \end{aligned}$$
(30)

Then

$$\begin{aligned} \Vert \varvec{\varepsilon }^\prime \Vert =\frac{-s_\pm \frac{f_\textrm{d}}{a} \Vert \textbf{T}^*\Vert ^2 \delta ^\prime _\textrm{s} +\sqrt{\textrm{Disc}}}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2},\quad \delta ^\prime _\textrm{s} =\delta ^\prime +\frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}}, \end{aligned}$$
(31)

with the positive discriminant

$$\begin{aligned} \textrm{Disc} =\Vert \textbf{T}^*\Vert ^2 (\delta ^\prime _\textrm{s})^2 +\frac{1}{3} \left( 1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \right) (\delta ^\prime )^2. \end{aligned}$$
(32)

Proof

Recalling the decomposition of the strain rate Eq. 11 (with the corresponding time transformation), from Eq. 28 it follows the identity

$$\begin{aligned} \varvec{\varepsilon }^\prime =\left( \delta ^\prime _\textrm{s} -s_\pm \frac{f_\textrm{d}}{a} \Vert \varvec{\varepsilon }^\prime \Vert \right) \textbf{T}^* +\frac{1}{3} \delta ^\prime \textbf{I}. \end{aligned}$$
(33)

This leads to a quadratic equation for \(\Vert \varvec{\varepsilon }^\prime \Vert\), more precisely:

$$\begin{aligned} \left( 1 -\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \right) \Vert \varvec{\varepsilon }^\prime \Vert ^2 +2 s_\pm \frac{f_\textrm{d}}{a} \Vert \textbf{T}^*\Vert ^2 \delta ^\prime _\textrm{s} \Vert \varvec{\varepsilon }^\prime \Vert -\Vert \textbf{T}^*\Vert ^2 (\delta ^\prime _\textrm{s})^2 -\frac{1}{3} (\delta ^\prime )^2 =0, \end{aligned}$$

whose discriminant is given by Eq. 32 and is positive thanks to Eq. 30. To see that Eq. 31 is indeed the positive solution of the equation above it is enough to observe that Eq. 30 implies \(\frac{f_\textrm{d}}{a} \Vert \textbf{T}^*\Vert ^2 \le \Vert \textbf{T}^*\Vert\), while \(\sqrt{\textrm{Disc}}\ge \Vert \textbf{T}^*\Vert |\delta ^\prime _\textrm{s}|\). \(\square\)

Theorem 1 allows us to write the deviatoric part of the strain rate, Eq. 28, in terms of the function \(\delta\) as follows

$$\begin{aligned} (\varvec{\varepsilon }^*)^\prime =k_\textrm{d} \textbf{T}^*,\quad k_\textrm{d} =\frac{\delta ^\prime _\textrm{s} -s_\pm \frac{f_\textrm{d}}{a} \sqrt{\textrm{Disc}}}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2}. \end{aligned}$$
(34)

Using this and Eq. 21, Eq. 29 can also be reduced to a quadratic differential equation in \(y^\prime\) whose analysis is the content of the following theorem.

Theorem 2

(Analytical solution) If Eqs. 27 and 30 hold, then the auxiliary variable y in Eq. 8 satisfies the quadratic differential equation

$$\begin{aligned} c_2(y) (y^\prime )^2 -2 c_1(y) \frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} y^\prime +c_0 \left( \frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \right) ^2 =0 \end{aligned}$$
(35)

with the inhomogeneous coefficients depending on y as follows

$$\begin{aligned} {\begin{matrix} c_2(y) &{}=\left( A(y) +\frac{B}{1+\frac{\exp (y-y_0)}{e_0}} \right) ^2 -D \left( \frac{1}{1+\frac{\exp (y-y_0)}{e_0}} \right) ^2 C^2,\\ c_1(y) &{}=\left( A(y) +\frac{B}{1+\frac{\exp (y-y_0)}{e_0}} \right) B -\frac{1}{1+\frac{\exp (y-y_0)}{e_0}} \Vert \textbf{T}^*\Vert ^2 C^2, \end{matrix}} \end{aligned}$$
(36)

where for brevity

$$\begin{aligned} A(y) =\frac{1}{f_\textrm{s} n y} +\frac{a^2 +\frac{1}{3}}{1+\frac{\exp (y-y_0)}{e_0}}, \end{aligned}$$

and \(c_0, B, C, D\) are constant factors given by

$$\begin{aligned} {\begin{matrix} c_0 =B^2 -\Vert \textbf{T}^*\Vert ^2 C^2, \quad &{}B =\frac{(1-f_\textrm{d}^2) \Vert \textbf{T}^*\Vert ^2}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2}, \\ C =\frac{f_\textrm{d}}{a} \frac{a^2 -\Vert \textbf{T}^*\Vert ^2}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2}, \quad &{}D =\Vert \textbf{T}^*\Vert ^2 +\frac{1}{3} \left( 1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \right) . \end{matrix}} \end{aligned}$$
(37)

Assume that

$$\begin{aligned} c_0 \left( 1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \right) \ge 0, \end{aligned}$$
(38)

then Eq. 35 has two solutions, \(y_+\) and \(y_-\), which satisfy respectively the following two differential equations:

$$\begin{aligned} y^\prime =\frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \frac{c_1(y) +\sqrt{c_1^2(y)-c_2(y) c_0}}{c_2(y)}, \end{aligned}$$
(39)
$$\begin{aligned} y^\prime =\frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \frac{c_1(y) -\sqrt{c_1^2(y)-c_2(y) c_0}}{c_2(y)}. \end{aligned}$$
(40)

Proof

Using Eqs. 31 and 34, from Eq. 29, we get the equation

$$\begin{aligned} \frac{y^\prime }{f_\textrm{s} n y} -\left( a^2 +\frac{1}{3} +\frac{(1-f_\textrm{d}^2) \Vert \textbf{T}^*\Vert ^2}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2} \right) \delta ^\prime -\frac{(1-f_\textrm{d}^2) \Vert \textbf{T}^*\Vert ^2}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2} \frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}}=s_\pm \frac{f_\textrm{d}}{a} \frac{a^2 -\Vert \textbf{T}^*\Vert ^2}{1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2} \sqrt{\textrm{Disc}}. \end{aligned}$$

Recalling that by Eq. 21\(\delta ^\prime =-y^\prime /(1+\frac{\exp (y-y_0)}{e_0})\), we rewrite the identity above using the notation in Eq. 37 to obtain

$$\begin{aligned} \left( A(y) +\frac{B}{1+\frac{\exp (y-y_0)}{e_0}} \right) y^\prime -B \frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} =s_\pm C \sqrt{\textrm{Disc}}, \end{aligned}$$
(41)

where \(\textrm{Disc}\), given by Eq. 32, corresponds to

$$\begin{aligned} \textrm{Disc} =D (\delta ^\prime )^2 +2 \Vert \textbf{T}^*\Vert ^2 \frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \delta ^\prime +\Vert \textbf{T}^*\Vert ^2 \left( \frac{1-\kappa }{f_\textrm{s} a^2} \frac{h^\prime _\textrm{s}}{h_\textrm{s}} \right) ^2. \end{aligned}$$
(42)

Squaring Eq. 41, after gathering alike terms, we get the quadratic Eq. 35 for \(y^\prime\) with the corresponding coefficients \(c_2(y), c_1(y), c_0\) from Eqs. 36-37. Note that

$$\begin{aligned} c_1^2(y) -c_2(y) c_0&=\left[ \left( A(y) +\frac{B}{1+\frac{\exp (y-y_0)}{e_0}} \right) B -\frac{1}{1+\frac{\exp (y-y_0)}{e_0}} \Vert \textbf{T}^*\Vert ^2 C^2 \right] ^2\\&-\left( A(y) +\frac{B}{1+\frac{\exp (y-y_0)}{e_0}} \right) ^2 (B^2 -\Vert \textbf{T}^*\Vert ^2 C^2)\\&-\left[ \Vert \textbf{T}^*\Vert ^2 +\frac{1}{3} \left( 1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \right) \right] \left( \frac{1}{1+\frac{\exp (y-y_0)}{e_0}} \right) ^2 C^2 (B^2 -\Vert \textbf{T}^*\Vert ^2 C^2) \end{aligned}$$

after some simplification yields to

$$\begin{aligned} c_1^2(y) -c_2(y) c_0 =C^2 \left[ \Vert \textbf{T}^*\Vert ^2 A^2(y) +\frac{c_0}{3} \left( 1-\frac{f_\textrm{d}^2}{a^2} \Vert \textbf{T}^*\Vert ^2 \right) \left( \frac{1}{1+\frac{\exp (y-y_0)}{e_0}} \right) ^2 \right] . \end{aligned}$$
(43)

Therefore, assumption Eq. 38 ensures that the discriminant of Eq. 35 is positive. Further, Eq. 38 together with Eq. 30 implies \(c_0\ge 0\), hence the quadratic Eq. 35 has two solutions, Eqs. 39 and 40, and the proof is complete. \(\square\)

Finally, by solving the differential Eqs. 39 or 40 with respect to y and recalling the relation \(\delta ^\prime =-y^\prime /(1+\frac{\exp (y-y_0)}{e_0})\), we obtain an explicit formula for the strain rate via Eq. 34.

Long-time behavior of the stress-controlled case

In the stress-controlled case analyzed in Theorem 2, we can make the following statements about the long-time behavior of the model:

Theorem 3

(Attractor) Under the stress control Eq. 27 and assumption Eq. 30, as \(t\rightarrow \infty\) the system Eqs. 9, 14, 28, and 29 tends exponentially to a stationary state given by constant

$$\begin{aligned} h_\textrm{s} =h_\textrm{s}^\infty ,\quad p =p_\infty ,\quad e=e_\infty ,\quad \delta =\delta _\infty . \end{aligned}$$
(44)

Proof

Indeed, as \(t\rightarrow \infty\) in Eq. 2, we have \(h_\textrm{s}(t)\rightarrow h_\textrm{s}^\infty\) and \(\dot{h}_\textrm{s} =0\) in the limit. Then \(h^\prime _\textrm{s} =0\) and, by Theorem 1, the norm \(\Vert \varvec{\varepsilon }^\prime \Vert\) is given by Eq. 31 with \(\delta ^\prime _\textrm{s} =\delta ^\prime\). Moreover, for \(h^\prime _\textrm{s} =0\), the discriminant in Eq. 42 corresponds to \(\sqrt{\textrm{Disc}} =\sqrt{D} |\delta ^\prime |\) and we can factorize Eq. 41 as follows

$$\begin{aligned} \left( A(y) +\frac{B +s_\pm C \sqrt{D}}{1+\frac{\exp (y-y_0)}{e_0}} \right) y^\prime =0. \end{aligned}$$
(45)

This implies that, in the limit, y is constant. We claim that this convergence is exponential, and as a consequence of Eqs. 10 and 15, we conclude Eq. 44. Indeed, note that equations Eqs. 39 and 40 are of the form \(y^\prime =(h^\prime _\textrm{s}/ h_\textrm{s}) F(y)\) for some bounded function F. Observing that \(h^\prime _\textrm{s}/ h_\textrm{s}\) behaves as the exponential function \(\frac{\exp (-t)}{\gamma +\exp (-t)}\), with \(\gamma >0\) constant, we need to distinguish two cases: \(F(y_0)<0\) and \(F(y_0)>0\), for a prescribed \(y(0)=y_0\). The case \(F(y_0)=0\) is trivial. If \(F(y_0)<0\) the exponential decay is clear. On the other hand, if \(F(y_0)>0\), then the convergence is driven by the exponential term which multiplies F. In summary, in either case, one can show that y converges to a constant at an exponential rate which proves the result. \(\square\)

For an illustrative example, we consider the normalized stress matrix in Eq. 27 of the shear form:

$$\begin{aligned} {\textbf {T}} =\begin{pmatrix} \frac{2}{5}&{}-\frac{2}{5}&{}0\\[.5ex]-\frac{2}{5}&{}\frac{2}{5}&{}0\\[.5ex]0&{}0&{}\frac{1}{5} \end{pmatrix},\quad {\textbf {T}}^* =\begin{pmatrix} \frac{1}{15}&{}-\frac{2}{5}&{}0\\[.5ex]-\frac{2}{5}&{}\frac{1}{15}&{}0\\[.5ex]0&{}0&{}-\frac{2}{15} \end{pmatrix} \end{aligned}$$

such that the principal eigenvalues of \({\textbf {T}}\) are \(\{0.8, 0.2, 0\}\), and \(\Vert {\textbf {T}}^*\Vert \approx 0.588\). The material parameters are taken from [6, 7] as follows: \(h_\textrm{s}^\infty =78.5\) MPa, \(h_\textrm{s}^0=120\) MPa and \(c=4\) hours for the degradation in Eq. 1; the exponents \(n =0.82\) in Eq. 8, and \(\alpha =0.18\), \(\beta =1.05\), \(b=1\) in Eq. 16; the void ratio bounds \(e_\textrm{min}=0.1\) and \(e_\textrm{max}=0.3\) in Eq. 7, \(e_\textrm{crt}=0.24\) in Eq. 17; \(a=0.6\) and \(\kappa =0.6\) in the hypoplastic law of Eq. 4. The initial values at \(t_0=0\) are prescribed as \(e_0=0.2\), \(p_0=10\) MPa, and \(\delta _0=0\) such that the factors \(f_\textrm{d}\approx 1.561\) and \(f_\textrm{s}\approx -0.844\) in Eq. 16. These data yield the coefficient \(c_0\approx 0.136\) in Eq. 36 and satisfy the solvability condition of Eq. 38 of Theorem 2.

We discretize the system of equations based on equidistant time points. To compute solutions \(y=y_+\) to the first-order differential Eq. 39 and \(y=y_-\) to Eq. 40, the standard fourth-order Runge–Kutta method is applied with the uniform step of 1 (hour).

A monotone behavior as \(t\rightarrow \infty\) in the hypoplastic model is described by the unified factor \(s_\pm =\pm 1\) implying the loading parameter \(s(t)=s_\pm t\) versus time t. This means that \(s_\pm =1\) when \(\delta\) is increasing, and \(s_\pm =-1\) when \(\delta\) decreases. In this case, a numerical solution to the system Eqs. 30–40 is plotted in Fig. 1.

Fig. 1
figure 1

An example solution under monotone loading versus time t (hours)

Full size image

The plots from (a) to (f) show respectively: the granular hardness \(h_\textrm{s}\), auxiliary variable y, pressure p, void ratio e, volume strain \(\delta\), and the deviatoric strain rate factor \(k_\textrm{d}\) for \(t\in [0,50]\). As seen in the figures, the spherical components \(y, p, e, \delta\) coincide for the loading/unloading, and they demonstrate two different scenarios of degradation corresponding either \(y_+\) or \(y_-\). In contrast, the deviatoric component \(k_\textrm{d}\) is the same for \(y=y_+\) and \(y=y_-\), whereas two curves in the plot (f) distinguish either loading as \(s_\pm =1\) or unloading as \(s_\pm =-1\) in Eq. 34.

Cyclic loading-unloading

In this section, we aim to investigate the model behavior in the stress-controlled case when the loading-unloading cycles are given by a continuous function s(t) such that

$$\begin{aligned} {\left\{ \begin{array}{ll} \dot{s}(t)>0 &{}\text {for}\;t\in (t_{2j}, t_{2j+1}),\\ \dot{s}(t)<0 &{}\text {for}\;t\in (t_{2j+1}, t_{2j+2}), \end{array}\right. }\quad j=0,1,2,\ldots \end{aligned}$$

with an increasing sequence of turning points \(t_{k+1}>t_k\) for \(k =0,1,2,\ldots\). For instance, starting with \(s(t_0)=0\), we set \(t_{k+1}=2k-1\) and the piecewise-linear periodic function

$$\begin{aligned} {\left\{ \begin{array}{ll} s(t) =s(t_{2j}) +t-t_{2j} &{}\text {for}\;t\in (t_{2j}, t_{2j+1}),\\ s(t) =s(t_{2j+1}) +t_{2j+1} -t &{}\text {for}\;t\in (t_{2j+1}, t_{2j+2}) \end{array}\right. } \end{aligned}$$
(46)

such that \(\dot{s}=1\) for loading, and \(\dot{s}=-1\) for unloading, which is illustrated in the plot (a) of Fig. 2 during 10 cycles.

Fig. 2
figure 2

Spherical components under cyclic loading-unloading

Full size image

After solving the hypoplastic system Eqs. 30–40 numerically as described before, in Fig.  2b, we depict the corresponding void ratio e versus pressure p, while the stress–strain curve presenting the trace of the stress tensor \(\textrm{tr}\varvec{\sigma } =-3p\) versus volume strain \(\delta\) is plotted in Fig. 2c. Recalling that assumption of Eq. 38 yields two solutions, if \(y_+\not = y_-\), we get distinct \(\delta (y_+)\not =\delta (y_-)\) in Eq. 21. As a consequence Fig. 2b and c describe two possible scenarios of the granular body sparsification corresponding to the two roots \(y_+\) and \(y_-\) of the quadratic differential Eq. 35. Indeed, in Fig. 2b, we can see that the curves identifying the void ratio in either case, \(e(y_+)\) and \(e(y_-)\), approach one another as the pressure increases. Similarly, Fig. 2c shows that at higher pressure (therefore, smaller values of \(-3p\)) the volume strain has basically the same curvature for both \(\delta (y_+)\) and \(\delta (y_-)\).

To study the long-time behavior of the deviatoric strain rate factor \(k_\textrm{d}\), we consider the process of loading-unloading with an increasing number of cycles. The deviatoric strain rate factor \(k_\textrm{d}\) during 2, 3, and 10 loading-unloading cycles is shown within the plots (a), (b), and (c) in Fig. 3, respectively. The dynamic presents a square-spiral curve, in which hysteresis loops are not closed and obey ratcheting. This limit cyclic behavior can be explained by the fact that \(h_s^\infty >0\), which is in accordance with the expected limiting behavior of the granular degradation of Eq. 2 which is rate dependent. Despite that, when the number of cycles increases, the loops converge to an attractor set in accordance with Theorem 3.

Fig. 3
figure 3

The deviatoric component under cyclic loading-unloading with creep time \(c=4\)

Full size image

In order to highlight the role of degradation, we vary the creep time c in the exponential function Eq. 2 for the granular hardness \(h_\textrm{s}\). Compared to \(c =4\) hours from Fig. 3, increasing the creep time to \(c =20\) hours decreases the deviation of the solution and suppresses its hysteresis ratcheting as presented in Fig. 4.

Fig. 4
figure 4

The deviatoric component under cyclic loading-unloading with creep time \(c=20\)

Full size image

From our numerical tests, we found that an increase in the value of c delays the degradation, leading to a kind of limiting cycle which can be observed, while the convergence to a stationary constant state becomes visible just with a higher number of cycles.