On feasibility of rate-independent stress paths under proportional deformations within hypoplastic constitutive model for granular materials

1. Introduction

The constitutive stress-strain 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 non-linear, 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 [2-4], 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 [2-4]. 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 non-linear 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 non-linear character of the underlying ODE system. Following the rate-independent 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 non-linear 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 closed-form 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 3-by-3 tensors of second order, which are undersigned with the upper tilde, we assume that the strain tensor $\tilde{\epsilon }$ and the corresponding constant quantity ${\tilde{\epsilon }}^{\text{'}}=d\tilde{\epsilon }/ds$ are given by proportional paths with respect to the monotonic increasing parameter $s\ge 0$ as:

The tensor $\tilde{U}$ is normalized by the Frobenius norm and determines a prescribed constant direction of the proportional strain path. Based on Eq. (1), the rate-independent hypoplastic relations take the form of the following non-linear ODE system with respect to the unknown symmetric Cauchy stress tensor $\tilde{\sigma }\left(s\right)$ (see [11-13] for more modelling issues):

where $\tilde{\sigma }:\tilde{U}$ stands for the scalar product of tensors, and the trace implies $tr\left(\tilde{\sigma }\right)=\tilde{\sigma }:\tilde{I}$. The dimensionless tensor $\tilde{V}$ is spanned by $\tilde{U}$ and the identity tensor $\tilde{I}$ as follows:

The initial condition for Eq. (2) prescribes the initial stress tensor ${\tilde{\sigma }}^0$ at $s=0$:

From a physical point of view, the constitutive parameter $c$ in Eq. (2) scales the amount of ${\tilde{\sigma }}^{\text{'}}$ and the second constitutive parameter $a$ characterizes the shape of the conical limit stress surface or so-called 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 closed-form solution to the initial value problem given by Eqs. (2-4) was found in [12-13]:

where the values $E$, $D$, and $h\left(s\right)$ are defined as follows:

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. (5-6). Set the parameter values:

$\begin{array}{cc}\begin{array}{cc}c=-1,& a=0.35,\end{array}& \tilde{U}=\frac{\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(1,2,-10\right)}{\sqrt{105}}\approx \mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}(0.098,0.195,-0.976),\end{array}$

such that $‖\tilde{U}‖=1$ and $tr\left(\tilde{U}\right)=-7/\sqrt{105}\approx -0.683$. Eq. (3) determines the tensor:

$\tilde{V}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\frac{0.35}{\sqrt{105}}-\frac{1}{3},\frac{0.7}{\sqrt{105}}-\frac{1}{3},-\frac{3.5}{\sqrt{105}}-\frac{1}{3}\right)\approx -\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(0.299,0.265,0.675\right)<0,$

with $tr\left(\tilde{V}\right)=-2.45/\sqrt{105}-1\approx -1.239$. For the value of the initial stress tensor we choose ${\tilde{\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:

$\begin{array}{cc}\begin{array}{cc}E\approx 0.234,& D\approx -0.2,\end{array}& \begin{array}{cc}C\approx 0.32,& h\left(s\right)\approx 0.7375\left(e^{-0.434s}-1\right).\end{array}\end{array}$

As the corresponding material behavior is coaxial, the stress tensor has zero elements outside the diagonal, thus $\tilde{\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:

$\tilde{\sigma }\left(s\right)\approx \left(e^{-0.234s}{\tilde{\sigma }}^0+1.614\left(e^{0.2s}-e^{-0.234s}\right)\tilde{V}\right)e^{0.7375\left(\mathrm{e}\mathrm{x}\mathrm{p}\left(-0.434s\right)-1\right)}<0.$

In Fig. 1 (a) we show for coaxial deformation the diagonal of $\tilde{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 $\tilde{\epsilon }\left(s\right)=s\tilde{U}$ from Eq. (1) forms a half-line along $\tilde{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 $\tilde{V}$, ${\tilde{\sigma }}^0$, $\tilde{\sigma }$ are depicted in the coordinates $-{\sigma }_1=$$-{\sigma }_{11}$, $-{\sigma }_2=$$-{\sigma }_{22}$, $-{\sigma }_3=$$-{\sigma }_{33}$. Starting as $s=0$ at ${\tilde{\sigma }}^0$ marked by the cross, the stress path $\tilde{\sigma }\left(s\right)$ (marked by a solid curve) approaches the half-line along a vector $\tilde{V}$ as $s\to \infty$. The asymptotic behavior was investigated in our previous works [11-13]. In the present context, it is worth noting that $\tilde{V}$, ${\tilde{\sigma }}^0$ and $\tilde{\sigma }\left(s\right)$ for all $s>0$ lie in the negative octant in Fig. 1(b) (although not $\tilde{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. (5-6), the stress path $\tilde{\sigma }\left(s\right)$ lies in the hyperplane spanned by ${\tilde{\sigma }}^0$ and $\tilde{V}$. Therefore, we define a cone of feasible directions depending on $a$ and $\tilde{U}$ such that:

where $v_1$, $v_2$, $v_3$ are eigenvalues of $\tilde{V}$. That means that $K_V$ is the set of all directions $\tilde{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 $\tilde{\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 so-called deviator plane is defined as the plane for $tr\left(\tilde{\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(\tilde{\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(\tilde{V}\right)=v_1+v_2+v_3<0$ for feasible $\tilde{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 ${\tilde{\sigma }}^0$ and $\tilde{V}$ in Eq. (5) are positive. Choosing both $\tilde{V}\in K_V$ and the initial stress state ${\tilde{\sigma }}^0\in K_{\sigma }$ thus guarantees that the whole stress path $\tilde{\sigma }\left(s\right)$ is contained in $K_{\sigma }$ for all $s\ge 0$.

We now investigate necessary and sufficient conditions on $a$ and $\tilde{U}$ to ensure that $\tilde{V}\in K_V$.

Fig. 2Representation of the intersection of the limit cones with the deviator plane in the principal stress space

3.2. Upper boundary of admissible parameters

The question is more delicate if the pair $(tr\left(\tilde{U}\right),a)$ lies above the dot-dashed 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 ${\tilde{V}}^{*}$ at the distance $r(a,\tilde{U})$ from the center of the triangle fall out of $K_V$. The points ${\tilde{V}}^{*}$ with $r\left(a,\tilde{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:

13

$r\left(a,\tilde{U}\right)=a\frac{\sqrt{1-{tr}^2\left(\tilde{U}\right)/3}}{1-atr\left(\tilde{U}\right)}>r_c,$

still under the constraint $1-atr\left(\tilde{U}\right)>0$, and obtain:

14

$a\left(r_{c}tr\left(\tilde{U}\right)+\sqrt{1-\frac{1}{3}{tr}^2\left(\tilde{U}\right)}\right)>r_c.$

The critical value of $tr\left(\tilde{U}\right)$ is such that:

$r_{c}tr\left(\tilde{U}\right)+\sqrt{1-\frac{1}{3}{tr}^2\left(\tilde{U}\right)}=0,$

that is $tr\left(\tilde{U}\right)=-1.$ If $tr\left(\tilde{U}\right)\le -1$, then Eq. (14) is never satisfied for any $a>0$. Otherwise, for $tr\left(\tilde{U}\right)\in (-1,\sqrt{3}]$, only the values:

15

$a>\frac{r_c}{r_{c}tr\left(\tilde{U}\right)+\sqrt{1-{tr}^2\left(\tilde{U}\right)/3}},$

exclude the possible existence of points ${\tilde{V}}^{*}$at the distance $r(a,\tilde{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 non-linear 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 so-called ratchetting effects of fading and shifting hysteresis loops.