Abstract
We study stress paths that are obtained under proportional deformations within the rateindependent hypoplasticity theory of Kolymbas type describing granular materials like soil and broken rock. For a particular simplified hypoplastic constitutive model by Bauer, a closedform solution of the corresponding system of nonlinear ordinary differential equations is available. Since only negative principal stresses are relevant for the granular body, the feasibility of the solution consistent with physics is investigated in dependence of the direction of a proportional strain path and constitutive parameters of the model.
Highlights
 Explicit solution to the hypoplastic model of Kolymbas type in simplified version by Bauer is obtained in a closed form
 Mathematically confirmed that, under monotonic deformation, stress paths tend asymptotically towards an independent stress state
 Location of principal stresses in the negative octant is proved, which is a physical requirement for cohesionless granular materials
 Admissible range of the values of constitutive constants to model the asymptotic behavior of stress paths is represented in a diagram
 This study provides reasonable analytic solutions for numerical validation tests
1. Introduction
The constitutive stressstrain relation for hypoplastic granular materials like cohesionless soil or broken rock is under our consideration. The respective constitutive law is of the rate type, incrementally nonlinear, and it is based on the concept of hypoplasticity proposed in 90th by Kolymbas [1]. Compared to incrementally linear constitutive equations, e.g. for hyperelastic and hypoelastic material laws, the hypoplastic constitutive equations are not differentiable at zero strain rate. This results in different stiffnesses for loading and unloading which is typical for inelastic materials. In contrast to the classical elastoplastic concept, the strain in the theory of hypoplasticity is not decomposed into the elastic and plastic parts. For modelling of soil behavior with hypoplasticity we cite [24], and refer to [5, 6] for other mathematical approaches. The underlying relations lie within the general context of implicit constitutive theories, see e.g. introduction in [7] and references therein.
An intrinsic property of granular materials observed in experiments is that under monotonic deformation the corresponding stress paths tends asymptotically towards a stress state, which is independent on the initial stress state [24]. As asymptotic states are well defined, such states are preferably qualified for the calibration of the constitutive constants. The present paper presents a novel investigation of the possible range of the values of constitutive constants to model the asymptotic behavior of stress paths under proportional deformation in hypoplasticity. Particular attention is paid on the restriction of stress states located in the negative octant of the space of principal stresses which is, from the physical point of view, a general requirement for cohesionless granular materials. Such behavior of granular materials is of practical importance for constitutive modelling not only in hypoplasticity, but, more widely, in soil mechanics and geotechnical engineering. This study also provides reasonable analytical solutions for numerical validation tests, thus it has a significant impact on efficient numerical simulations.
In the present paper, a simplified version of the hypoplastic model described in Bauer [8] and Gudehus [9] is considered. For the sake of simplicity, the pressure and density dependent properties of granular materials are omitted so that only two constitutive parameters remain in the model. Further relying on proportional strain paths, the rate model is reduced to a system of nonlinear ordinary differential equations (ODE). As the result of these simplifications, the stress path is determined by the direction of a proportional strain path, by the initial stress state and the two constitutive parameters involved in the hypoplastic equation.
The main challenge in developing proper mathematical tools is the strongly nonlinear character of the underlying ODE system. Following the rateindependent technique developed in [10], for the simplified hypoplastic model, the asymptotic behavior was proved in our previous work [11]. Then in [12] we constructed the solution of the corresponding nonlinear problem in a closed form, and extended it to a modified model in [13]. The exact solution allows us to describe analytically the evolution of stress paths for various loading scenarios obtained from monotonic compression, extension, and isochoric deformations.
Since only negative principal stresses are relevant for the granular body, it is important to discuss the feasibility of a solution which is consistent with physics. In the present contribution, on the basis of a closedform solution we identify the domain of the constitutive parameters which guarantee that the corresponding stress path starting from feasible initial stress state remains feasible under proportional deformation paths.
2. Theory
In the space of symmetric 3by3 tensors of second order, which are undersigned with the upper tilde, we assume that the strain tensor $\stackrel{~}{\epsilon}$ and the corresponding constant quantity ${\stackrel{~}{\epsilon}}^{\text{'}}=d\stackrel{~}{\epsilon}/ds$ are given by proportional paths with respect to the monotonic increasing parameter $s\ge 0$ as:
The tensor $\stackrel{~}{U}$ is normalized by the Frobenius norm and determines a prescribed constant direction of the proportional strain path. Based on Eq. (1), the rateindependent hypoplastic relations take the form of the following nonlinear ODE system with respect to the unknown symmetric Cauchy stress tensor $\stackrel{~}{\sigma}\left(s\right)$ (see [1113] for more modelling issues):
where $\stackrel{~}{\sigma}:\stackrel{~}{U}$ stands for the scalar product of tensors, and the trace implies $tr\left(\stackrel{~}{\sigma}\right)=\stackrel{~}{\sigma}:\stackrel{~}{I}$. The dimensionless tensor $\stackrel{~}{V}$ is spanned by $\stackrel{~}{U}$ and the identity tensor $\stackrel{~}{I}$ as follows:
The initial condition for Eq. (2) prescribes the initial stress tensor ${\stackrel{~}{\sigma}}^{0}$ at $s=0$:
From a physical point of view, the constitutive parameter $c$ in Eq. (2) scales the amount of ${\stackrel{~}{\sigma}}^{\text{'}}$ and the second constitutive parameter $a$ characterizes the shape of the conical limit stress surface or socalled critical stress state surface in the principal stress space [8]. By this we have assumed deformations for which the objective time derivative coincides with the material time derivative, that holds for example for coaxial deformations.
A closedform solution to the initial value problem given by Eqs. (24) was found in [1213]:
where the values $E$, $D$, and $h\left(s\right)$ are defined as follows:
$\begin{array}{cc}h\left(s\right):=\frac{C}{atr\left(\stackrel{~}{V}\right)}\left({e}^{catr\left(\stackrel{~}{V}\right)s}1\right),& C:=\frac{{\stackrel{~}{\sigma}}^{0}:\stackrel{~}{U}}{tr\left({\stackrel{~}{\sigma}}^{0}\right)}\frac{a\frac{1}{3}tr\left(\stackrel{~}{U}\right)}{tr\left(\stackrel{~}{V}\right)}.\end{array}$
Since only negative principal stresses ${\sigma}_{1}\left(s\right)$, ${\sigma}_{2}\left(s\right)$, ${\sigma}_{3}\left(s\right)$ are relevant for the underlying physical model of granular solids, it is important to discuss their feasibility.
2.1. Computational simulation
We present a computational example of coaxial deformation according to Eqs. (56). Set the parameter values:
such that $\Vert \stackrel{~}{U}\Vert =1$ and $tr\left(\stackrel{~}{U}\right)=7/\sqrt{105}\approx 0.683$. Eq. (3) determines the tensor:
with $tr\left(\stackrel{~}{V}\right)=2.45/\sqrt{105}1\approx 1.239$. For the value of the initial stress tensor we choose ${\stackrel{~}{\sigma}}^{0}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(1,0.5,0.5\right)<0$. Then the parameters from Eq. (6) are calculated as:
As the corresponding material behavior is coaxial, the stress tensor has zero elements outside the diagonal, thus $\stackrel{~}{\sigma}\left(s\right)=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\sigma}_{1}\right(s),{\sigma}_{2}(s),{\sigma}_{3}(s\left)\right)$, and from Eq. (5) we get:
In Fig. 1 (a) we show for coaxial deformation the diagonal of $\stackrel{~}{U}$ as a direction vector of the proportional deformation prescribed in the coordinate system ${\epsilon}_{1}=$${\epsilon}_{11}$, ${\epsilon}_{2}=$${\epsilon}_{22}$, ${\epsilon}_{3}=$${\epsilon}_{33}$. Then the proportional strain path $\stackrel{~}{\epsilon}\left(s\right)=s\stackrel{~}{U}$ from Eq. (1) forms a halfline along $\stackrel{~}{U}$ for $s\ge 0$. The dashed line here is the isotropic axis. Intersection of a perpendicular plane with the negative octant constitutes an equilateral triangle drawn for visualization reason. In Fig. 1(b) the diagonals of tensors $\stackrel{~}{V}$, ${\stackrel{~}{\sigma}}^{0}$, $\stackrel{~}{\sigma}$ are depicted in the coordinates ${\sigma}_{1}=$${\sigma}_{11}$, ${\sigma}_{2}=$${\sigma}_{22}$, ${\sigma}_{3}=$${\sigma}_{33}$. Starting as $s=0$ at ${\stackrel{~}{\sigma}}^{0}$ marked by the cross, the stress path $\stackrel{~}{\sigma}\left(s\right)$ (marked by a solid curve) approaches the halfline along a vector $\stackrel{~}{V}$ as $s\to \infty $. The asymptotic behavior was investigated in our previous works [1113]. In the present context, it is worth noting that $\stackrel{~}{V}$, ${\stackrel{~}{\sigma}}^{0}$ and $\stackrel{~}{\sigma}\left(s\right)$ for all $s>0$ lie in the negative octant in Fig. 1(b) (although not $\stackrel{~}{U}$ in Fig. 1(a)), thus implying a feasible solution consistent with physics.
Fig. 1Computational simulation example: a) orientation of the proportional strain path prescribed and b) asymptotic behavior of the corresponding stress path
a) Principal strain space
b) Principal stress space
3. Feasibility conditions
According to the analytic solution from Eqs. (56), the stress path $\stackrel{~}{\sigma}\left(s\right)$ lies in the hyperplane spanned by ${\stackrel{~}{\sigma}}^{0}$ and $\stackrel{~}{V}$. Therefore, we define a cone of feasible directions depending on $a$ and $\stackrel{~}{U}$ such that:
where ${v}_{1}$, ${v}_{2}$, ${v}_{3}$ are eigenvalues of $\stackrel{~}{V}$. That means that ${K}_{V}$ is the set of all directions $\stackrel{~}{V}$ that lie in the negative octant and that can be reached by Eq. (3). On the other hand, we define a cone of feasible stresses such that:
where ${\sigma}_{1}$, ${\sigma}_{2}$, ${\sigma}_{3}$ imply the principal stresses corresponding to $\stackrel{~}{\sigma}$. We draw both ${K}_{V}$ and ${K}_{\sigma}$ in the same Fig. 2.
In the space of principal stresses the apex of the cone ${K}_{\sigma}$ lies in the origin of the principal stress axes. In the stress space the socalled deviator plane is defined as the plane for $tr\left(\stackrel{~}{\sigma}\right)={\sigma}_{1}+{\sigma}_{2}+{\sigma}_{3}=$ –1, and the intersection of this plane with the negative octant constitutes an equilateral triangle with the side length of $\sqrt{2}\approx \text{1.414}$ as shown in Fig. 2. In the figure we use marks ${\widehat{\sigma}}_{ii}^{*}={\sigma}_{ii}/tr\left(\stackrel{~}{\sigma}\right)1/3$ for the axes $i=\mathrm{1,2},3$. Furthermore, the intersection of the deviator plane with a circular cone with apex at the origin and main axis perpendicular to the deviator plane is represented by a disc of some radius $r\ge 0$. We can draw the circumcircle with radius ${r}_{c}=\sqrt{2/3}\approx \text{0.8165}$ and the incircle with radius ${r}_{i}=\sqrt{1/6}\approx \text{0.408}$ of this triangle. The intersection of the above triangle with a disc of radius ${r}_{a}\ge 0$ forms the cone ${K}_{\sigma}$. The value of ${r}_{a}$ depends on $a$, see [14].
Note that $tr\left(\stackrel{~}{V}\right)={v}_{1}+{v}_{2}+{v}_{3}<0$ for feasible $\stackrel{~}{V}$ from Eq. (7). Therefore, if $E\ge D$ in Eq. (6), then $\mathrm{e}\mathrm{x}\mathrm{p}\left(cEs\right)\le \mathrm{e}\mathrm{x}\mathrm{p}\left(cDs\right)$ (recalling that $c<0$), so that the factors in front of ${\stackrel{~}{\sigma}}^{0}$ and $\stackrel{~}{V}$ in Eq. (5) are positive. Choosing both $\stackrel{~}{V}\in {K}_{V}$ and the initial stress state ${\stackrel{~}{\sigma}}^{0}\in {K}_{\sigma}$ thus guarantees that the whole stress path $\stackrel{~}{\sigma}\left(s\right)$ is contained in ${K}_{\sigma}$ for all $s\ge 0$.
We now investigate necessary and sufficient conditions on $a$ and $\stackrel{~}{U}$ to ensure that $\stackrel{~}{V}\in {K}_{V}$.
Fig. 2Representation of the intersection of the limit cones with the deviator plane in the principal stress space
3.1. Lower boundary of admissible parameters
Let $\stackrel{~}{V}$ be as in (3) and let us assume that $tr\left(\stackrel{~}{V}\right)=atr\left(\stackrel{~}{U}\right)1<0$. Then ${\stackrel{~}{V}}^{*}=\stackrel{~}{V}/tr\left(\stackrel{~}{V}\right)$ is the axial projection of $\stackrel{~}{V}$ onto the deviator plane, which for instance is marked with the cross in Fig. 2. Let $r(a,\stackrel{~}{U})$ denote its distance from the center of the triangle. By definition of ${\stackrel{~}{V}}^{*}$ we have:
The norm in the numerator can be evaluated using the condition $\Vert \stackrel{~}{U}\Vert =1$ (which implies that $\lefttr\left(\stackrel{~}{U}\right)\right\le \sqrt{3}$) as follows:
Finally, we define ${r}_{a}$ to be the minimal value such that:
for all $\stackrel{~}{U}$ such that $\Vert \stackrel{~}{U}\Vert =1$. The maximum in Eq. (9) over all admissible values of $\stackrel{~}{U}$ is reached for $tr\left(\stackrel{~}{U}\right)=3a$, that is:
Hence, a sufficient condition for $\stackrel{~}{V}$ to belong to ${K}_{V}$ independently of the strain rate direction $\stackrel{~}{U}$ is given by Eq. (10) for ${{r}_{a}=r}_{i}=\sqrt{1/6}$, that is:
For $a>1/3$, the condition $\stackrel{~}{V}\in {K}_{V}$ is fulfilled only for certain values of $tr\left(\stackrel{~}{U}\right)$, namely those for which $r\left(a,\stackrel{~}{U}\right)\le {r}_{i}$ according to Eq. (9). Taking into account the constraint $atr\left(\stackrel{~}{U}\right)<1$, we rewrite Eq. (9) as:
The critical value of $tr\left(\stackrel{~}{U}\right)$ is such that:
that is $tr\left(\stackrel{~}{U}\right)=\sqrt{2}.$ If $tr\left(\stackrel{~}{U}\right)\le \sqrt{2}$, then Eq. (11) holds for all $a>0$. Otherwise, for $tr\left(\stackrel{~}{U}\right)\in (\sqrt{2},\sqrt{3}]$, only the values:
are admissible, and the lower boundary of the admissible parameter domain is represented by the dotdashed curve in Fig. 3.
3.2. Upper boundary of admissible parameters
The question is more delicate if the pair $(tr\left(\stackrel{~}{U}\right),a)$ lies above the dotdashed line in Fig. 3. As long as ${r}_{a}$ in Eq. (10) is between ${r}_{i}$ and ${r}_{c}$ as in Fig. 2, the circle of radius ${r}_{a}$ intersects the triangle in the negative octant. In other words, some points ${\stackrel{~}{V}}^{*}$ at the distance $r(a,\stackrel{~}{U})$ from the center of the triangle fall out of ${K}_{V}$. The points ${\stackrel{~}{V}}^{*}$ with $r\left(a,\stackrel{~}{U}\right)>{r}_{c}=\sqrt{2/3}$ do not belong to ${K}_{V}$ at all. The case of ${r}_{a}>{r}_{c}$ can be treated similarly as above. We replace Eq. (9) with:
still under the constraint $1atr\left(\stackrel{~}{U}\right)>0$, and obtain:
The critical value of $tr\left(\stackrel{~}{U}\right)$ is such that:
that is $tr\left(\stackrel{~}{U}\right)=1.$ If $tr\left(\stackrel{~}{U}\right)\le 1$, then Eq. (14) is never satisfied for any $a>0$. Otherwise, for $tr\left(\stackrel{~}{U}\right)\in (1,\sqrt{3}]$, only the values:
exclude the possible existence of points ${\stackrel{~}{V}}^{*}$at the distance $r(a,\stackrel{~}{U})$ belonging to ${K}_{V}$. The dashed curve in Fig. 3 represents the upper boundary of this parameter domain.
Fig. 3Visualization of the domain of parameters with feasibility conditions
4. Conclusions
The simplified hypoplastic constitutive Eq. (2) considered in this paper is incrementally nonlinear and includes only two constitutive parameters: a stiffness parameter $c$ and a parameter $a$ related to the critical friction angle of a cohesionless granular material, which is defined for the axisymmetric triaxial compression in the limit state. Although the influence of the direction of the stress deviator on the critical stress state is neglected, i.e. the constitutive parameter $a$ is assumed to be constant, it does not mean a restriction of the general results drawn. A relevant relation between the value of parameter $a$ for the direction of the stress deviator of triaxial compression and the one for another direction can be obtained by interpolation.
For future work, we propose to study the simplified hypoplastic model under cyclic loading, when the granular solid exhibits the socalled ratchetting effects of fading and shifting hysteresis loops.
Acknowledgements
This work was supported by the OeAD Scientific and Technological Cooperation (WTZ CZ 01/2016 and 18/2020) financed by the Austrian Federal Ministry of Science, Research and Economy (BMWFW) and by the Czech Ministry of Education, Youth and Sports (MŠMT).
Further support by the Russian Foundation for Basic Research (RFBR) and JSPS research project 195150004 (V.A.K.) and by the European Regional Development Fund, Project No. CZ.02.1.01/0.0/0.0/16_019/0000778 (P.K.) is gratefully acknowledged.
References

Kolymbas D. An outline of hypoplasticity. Archive of Applied Mechanics, Vol. 61, 1991, p. 143151.

Gudehus G. Physical Soil Mechanics. Springer, Berlin, Heidelberg, 2011.

Kolymbas D., Medicus G. Genealogy of hypoplasticity and barodesy. International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 40, 2016, p. 25322550.

Mašín D. Modelling of Soil Behaviour with Hypoplasticity. Springer Nature, Switzerland, 2019.

Annin B. D., Kovtunenko V. A., Sadovskii V. M. Variational and hemivariational inequalities in mechanics of elastoplastic, granular media, and quasibrittle cracks. Analysis, Modelling, Optimization and Numerical Techniques. Springer Proceedings in Mathematics and Statistics, Vol. 121, 2015, p. 4956.

Khludnev A. M., Kovtunenko V. A. Analysis of Cracks in Solids. WITPress, Southampton, Boston, 2000.

Itou H., Kovtunenko V. A., Rajagopal K. R. Wellposedness of the problem of nonpenetrating cracks in elastic bodies whose material moduli depend on the mean normal stress. International Journal of Engineering Science, Vol. 136, 2019, p. 1725.

Bauer E. Calibration of a comprehensive hypoplastic model for granular materials. Soils and Foundations, Vol. 36, 1996, p. 1326.

Gudehus G. A comprehensive constitutive equation for granular materials. Soils and Foundations, Vol. 36, 1996, p. 112.

Brokate M., Krejčí P. Wellposedness of kinematic hardening models in elastoplasticity. RAIRO Modélisation Mathématique et Analyse Numérique, Vol. 32, 1998, p. 177209.

Kovtunenko V. A., Krejčí P., Bauer E., Siváková L., Zubkova A. V. On Lyapunov stability in hypoplasticity. Proceedings of Equadiff Conference, Bratislava, 2017, p. 107116.

Bauer E., Kovtunenko V. A., Krejčí P., Krenn N., Siváková L., Zubkova A. V. On proportional deformation paths in hypoplasticity. Acta Mechanica, 2020, https://doi.org/10.1007/s00707019025973.

Bauer E., Kovtunenko V. A., Krejčí P., Krenn N., Siváková L., Zubkova A. V. Modified model for proportional loading and unloading hypoplastic materials. Extended Abstracts Spring 2018. Trends in Mathematics, Vol. 11, 2019, p. 201210.

Bauer E. Conditions for embedding Casagrande’s critical state into hypoplasticity. Mechanics of CohesiveFrictional Materials, Vol. 5, 2000, p. 125148.