On proportional deformation paths in hypoplasticity

We investigate rate-independent stress paths under constant rate of strain within the hypoplasticity theory of Kolymbas type. For a particular simplified hypoplastic constitutive model, the exact solution of the corresponding system of nonlinear ordinary differential equations is obtained in analytical form. On its basis, the behaviour of stress paths is examined in dependence of the direction of the proportional strain paths and material parameters of the model.


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 nonlinear, and it is based on the hypoplastic concept proposed by Kolymbas [29]. 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 is due to different stiffnesses for loading and unloading typical for inelastic materials. In contrast to the classical elastoplastic concept, the strain in the theory of hypoplasticity is not decomposed into elastic and plastic parts. Detailed discussions about physical aspects of hypoplastic models can be found, for instance, in [20,31,36]; their response to cycling loading is studied in [9,41], shear localization in [7,10,16,23], critical states in [4,5,49], behaviour under undrained condition in [37], extension to a micropolar continuum in [23,24], and thermodynamic aspects in [26,42,43].
Experiments with granular materials show a typical behaviour of the stress path obtained under a monotonic proportional deformation path which was formulated in two rules by Goldscheider [18]. The first rule says that a proportional deformation path starting from a nearly stress-free state results in a proportional stress path. The second rule says that a proportional deformation path starting from an arbitrary stress state asymptotically converges to the same proportional stress path as the one obtained for the initially stress-free state. These asymptotic properties of granular materials can be interpreted by fading memory of the material [21]. This feature is also confirmed by further experiments, e.g. [13,44], and numerical simulations by the discrete element method, e.g. [25,34].
As the asymptotic behaviour of the stress path for proportional deformation is an intrinsic property of granular materials, it has also great importance for constitutive models relevant to frictional granular materials. It can be noted that in many constitutive models the asymptotic behaviour is only fulfilled for particular deformation paths. The lack of this property can lead to large deviation of the prediction of stress paths not used for calibration.
In constitutive modelling, the asymptotic property is the starting point for barodesy modelling in [38], and it is also an intrinsic feature in hypoplasticity [20,22,28,35,36,40]. With respect to a particular hypoplastic model by Kolymbas [28], stress paths were numerically investigated for initially axisymmetric stress states and axisymmetric proportional deformations. However, no general requirements for the asymptotic behaviour were formulated. A mathematical criterion providing a necessary condition for convergent asymptotic behaviour in incremental models was given first by Niemunis [40]. It requires that the normal distance of the generated stress path to the asymptote decreases with an increase in monotonic deformation, where the asymptote is obtained for the corresponding fixed strain rate starting from the initially stress-free state. The criterion by Niemunis is applicable to all rate-type models; however, it is restricted to proportional deformation paths with contractant volume strain behaviour.
In the present paper, a different strategy is proposed to examine the asymptotic behaviour using the analytical solution for the stress path, which depends on the direction of the proportional strain path, the initial stress state, and the material parameters involved in the constitutive equation. The main challenge in developing proper mathematical tools is the strongly nonlinear character of the underlying differential equations. Following the rate-independent technique developed in [11], we gave a rigorous mathematical proof of the existence of asymptotic states in [8,32]. This idea is further developed in the present work.
To this end, a simplified version of the hypoplastic model by Gudehus [19] and Bauer [3] is considered. For the sake of simplicity, the pressure-and density-dependent properties of granular materials described in [3,19] are omitted so that only two material parameters remain in the model. It can be noted that the same simplified version can also be obtained from the hypoplastic model by von Wolffersdorff [46] as shown in [6].
For the simplified hypoplastic model, we construct the solution of the corresponding nonlinear problem in a closed form. In this way, we make also a close link to barodesy models [30]. The exact solution allows us to describe analytically various scenarios of the behaviour of stress paths obtained from monotonic compression, extension, and isochoric deformations. The latter leads to stress limit states or the so-called critical stress states, which can be represented by a conical surface in the space of principal stresses. In particular, we identify the domain of the constitutive parameters which guarantee that for proportional deformation paths the corresponding stress paths starting from an arbitrary initial stress state are asymptotically stable.

The model
We first fix some basic tensor notation. By R 3×3 sym , we denote the space of symmetric 3-by-3 tensors of second- , which are endowed with the usual double dot product, the associated norm, and the trace, respectively: Here, {I } := {δ i j } 3 i, j=1 stands for the 3-by-3 identity matrix with the Kronecker symbol δ i j = 1 for i = j, and zero otherwise. In physical interpretation, σ corresponds to the Cauchy stress tensor, andε is the strain rate tensor. For time t ≥ 0, we interpret σ (t) andε (t) as time-dependent tensor-valued functions.
With respect to the normalized stress tensorσ = σ /tr(σ ), the general representation of the hypoplastic constitutive equation of the Kolymbas type can be written in the factorized form as the following tensor equation for the objective stress rate: is a symmetric fourth-order tensor, and the double dot product is to be interpreted as The dimensionless parameter c < 0 scales the incremental stiffness and can be calibrated, for instance, based on an isotropic compression test. The right-hand side of (1.1) is a homogeneous function of degree one in σ . Note that dry granular materials are cohesionless, so that only negative principal stresses are relevant to the constitutive equation (1.1). Furthermore, we remark that particular representations of the tensor functions in (1.1) are based on terms from the general representation theorem of isotropic tensor-valued functions [47]. Various explicit versions are proposed in the literature (e.g. [3,19,29,43,48,49]). In this paper, we consider a particular version of (1.1) proposed by Bauer in [5] in a simplified manner: with the normalized stress deviatorσ * =σ − I /3, where the symbol 4 I stands for the fourth-order identity tensor, the symbol ⊗ denotes the dyadic product of tensors, and the term in (1.1) which is linear inε can also be represented as: (1. 3) The constitutive constant a > 0 is called limit stress state parameter and characterizes the shape of the conical limit stress surface or the so-called critical stress state surface in the principal stress space [5]. Critical states are defined for a vanishing stress rate under continuous isochoric deformation. For critical stress states, parameter a equals the norm of the normalized stress deviator, i.e. a = σ * , and it can be related to the so-called critical friction angle [4]. While in the model by Gudehus [19] and Bauer [3] the value of a also depends on the orientation of the stress deviator, parameter a is assumed to be a constant in the present paper. For the granular friction angle φ ∈ (0, π/2) such that a = 2 √ 2/3 sin φ/(3 − sin φ), we get the physical restriction a < a phys = √ 2/3 ≈ 0.8165 as sin φ < 1.

Hypoplastic model under proportional deformation
Further we restrict ourselves to strain paths pointing in one fixed direction. Namely, we call a strain path proportional if there exists a time-independent symmetric second-order tensor U ∈ R 3×3 sym normalized by U = 1 and a scalar time-dependent function s(t) with s(0) = 0 andṡ(t) > 0 such that it holds (2.1) In (2.1), U determines a prescribed direction in the space of symmetric second-order tensors, and s represents a monotonic increasing parameter. Moreover, we assume deformations for which the objective time derivative and the material time derivative coincide, i.e.
• σ =σ , which, for example, holds for fixed directions of principal stresses. As the material behaviour described by Eq. (1) is isotropic, the strain rate tensor with fixed principal values, i.e. the strain increments are proportional and coaxial, corresponds to a stress rate tensorσ with zero elements outside the diagonal. Thus, the relevant system of ordinary differential equations (ODE) reduces to three.
In this case, inserting the chain rule d/dt =ṡd/ds and ε =ṡ according to (2.1) and dividing withṡ = 0, we can rewrite the rate-independent relations (1) with respect to s > 0 in the form: or, equivalently, without factorization: where the formulaσ +σ * = 2σ − I /3 was used, and .
In the principal stress space, the representation of the normalized quantityV indicates the direction of the asymptotic stress state for s → ∞. This observation will be specified below after the representation formula (5.6). Indeed, formula (5.6.1) shows that the evolution of σ takes place in the 2D plane generated by constant tensorsV and σ (0). We shall see that the scalar coefficient ofV is dominant and determines the asymptotic convergence or divergence of the stress path as s → ∞.
In system (2.3) of ordinary differential equations (ODE) for σ (s) ∈ R 3×3 sym , the last term in the right-hand side is nonlinear in σ . For (2.3), we prescribe the initial condition at s = 0: where σ 0 ∈ R 3×3 sym is a given initial stress tensor. In the following, we study the behaviour as s → ∞ of solutions σ (s) to the nonlinear initial value problem (2) in dependence of the direction U of proportional strain paths, the material parameter a, and the initial stress σ 0 .

Motivating example: isotropic compression/extension
In this Section an assembly of grains with a grain skeleton is considered under compressive stresses. In mathematical terms, it means that all principal stress components are negative. Here, by compression or extension we mean the dynamics of the process defined by (2.1), where the strain rate can be, respectively, negative or positive.
Since the model is valid only for compressive stresses, we necessarily have tr(σ (0)) < 0. Furthermore, c < 0 and D − < 0, so that the negative stress trace exponentially grows when the compression linearly increases along the direction of −I for all a > 0. Equation (2.3) turns into the form In view of (3.2), the solution of (3.3) reads If the deviator σ * (0) of the initial stress is zero, that is, the initial stress is located on the isotropic axis generated by I , then the whole stress path lies on the isotropic axis. Otherwise, if the initial stress deviator does not vanish, there are three possible scenarios for the evolution of σ (s) depending on whether

Scenario (a)
If E − > 0, then we have cE − < 0, and the stress path converges exponentially to the isotropic trajectory with increasing time, so that the model is exponentially stable with respect to small perturbations of the initial stress.

Scenario (b)
The situation is totally different when E − < 0. Then cE − > 0, and small perturbations of the initial stress produce exponentially large deviation of the stress path from the isotropic trajectory.
Scenario (c) In the limit case E − = 0, the distance of the solution trajectory from the isotropic axis remains constant. Thus, scenarios (b) and (c) do not fulfil the second law by Goldscheider. According to the sign of E − in (3.3) we see that for the stability of the model with respect to variations of the initial stress (scenario (a)) it needs a > a min , and the critical parameter value is a min = 1/(2 √ 3) ≈ 0.2887. For smaller values of a (scenario (b)), we are in the classical philosophical situation 1 described by Frank Brentano.
For illustration, in Fig. 1 we show the projection of the stress paths onto the so-called Rendulic plane for particular values a = 0.1, a = a min , and a = 1, and an initial stress σ 0 = −diag(0.2, 0.01, 0.2). We observe that: for a = 1 > a min (solid curve) according to scenario (a) the stress paths tend towards the isotropic axis (dot dashed line), which is, in this case, the asymptote; for a = 0.1 < a min (dashed curve), there is no asymptotic behaviour implying scenario (b); and for a = a min (dotted curve), the stress path is parallel to the isotropic axis implying scenario (c).
(ii) Isotropic extension: For an initially prestressed state of monotonic isotropic extension, we have U + = −U − = I / √ 3 resulting in the equations for tr(σ ): and, respectively, for σ : As given above, its solution can be written in the form In this case, we have D + > 0 and E + > 0, and the stress σ (s) in (4.3) decays exponentially to zero as s → ∞ independently of the choice of the initial condition. We remark that according to the derived formulas (3) and (4) (as well as in the general case (2.3)) the parameter c < 0 does affect neither the asymptotic states nor the shape of stress paths; rather, it influences how quick a proportional stress path is approached.
Motivated by this special example, it is our aim to investigate in the following the stress path under general proportional deformations.

Analytical solution of the hypoplastic equation for general proportional deformations
The principal difficulty of solving the ODE (2.3) in the general form is its non-linearity in the stress. Here, we apply the following procedure consisting of three steps: Step 1 We derive from (2.3) an ODE for the auxiliary scalar variable (σ : U )/tr(σ ) and find its solution; Step 2 Inserting the solution to Eq. (2.3) projected in the isotropic direction, we obtain a linear equation for tr(σ ) and solve it; Step 3 Substituting this solution in the constitutive relation (2.3), we find an expression of σ in the closed form.
A rigorous derivation of the analytical solution is given in the "Appendix." Below we summarize the resulting formulas and discuss similarities between hypoplasticity and barodesy models. Let V be as in (2.4) and assume tr(V ) = 0. To simplify the formulas, we introduce the constants (see (17.1) and (17.6)) and the function h(s) with h(0) = 0 and depending on the initial stress tensor σ 0 , defined by the formula (see (16.4) and (17.3)): In (17.1), we prove that the differential equation (2.3) can be transformed into the following equation: It is worth noting that (5.3) has the structure (up to a factorization) of a constitutive relation adopted in barodesy; see, for example, [38]: where R is a given direction of the proportional stress path and f and g are model parameters. The exponential expression of R is formal and admits the expansion in −α: The first two linear asymptotic terms in (5.5) with α = −3a yield exactly V in (2.4), and the nonlinear behaviour of R with respect to U is substituted in the Kolymbas-type hypoplastic model by the norm U (which is normalized to one here). In this context, it can be mentioned that with constant a Eq. (1) models the stress limit condition by Drucker-Prager. For a refined modelling of the stress limit condition for granular materials, factor a should depend on the orientation of the stress deviator, i.e. it should be a function of the so-called Lode angle, which allows the adaptation of arbitrary conically shaped stress limit conditions in the principal stress space as outlined in [3,4]. Taking into account the contribution O(α 2 ) in Eq. (5.5), the barodesy model describes a stress limit condition close to the one by Matsuoka-Nakai [39] as shown in [17,38]. In hypoplasticity, the stress limit condition by Matsuoka-Nakai can be predefined as shown, for instance, in [4,46]. In both models (barodesy and hypoplasticity), parameters α and a are related to the granular friction angle. Based on the representation (5.3) which is linear in σ , the solution of the initial value problem (2) can be found in the closed form (see (18.3)): with E and D defined as in (5.1), h(s) as in (5.2), and V as in (2.4), whereasV = V /tr(V ). For the special case tr(V ) = 0 , that is, tr(U ) = 1/a for a = 0, see formulas (19). Due to tr(V ) = 1, from (5.6) we compute tr(σ )(s) = tr(σ 0 ) exp(cDs + h(s)) and the normalized stress tensor Therefore, D − E < 0 in (5.7) proves the asymptotic convergence exponentially as s → ∞.
Since only negative principal stresses are relevant for the underlying model, it is important to discuss their feasibility. For this reason, we define the feasible cone: where σ 1 , σ 2 , and σ 3 are eigenvalues of σ . Note that tr( In the next Sections, we investigate the asymptotic behaviour as s → ∞ of stress paths (5.6) in dependence of specific deformations U .

Asymptotic behaviour of stress paths under various proportional deformations
Similarly to the motivation example, we should distinguish contractant (volume-decreasing) from dilatant (volume-increasing) states. For this task, we call by contractant the compression corresponding to tr(U ) < 0, and by dilatant the extension when tr(U ) > 0, while tr(U ) = 0 corresponds to the volume-preserving deformation.
To this end, we note that paths with tr(U ) > 0 will eventually approach the origin as they are proportional volume-increasing deformation paths. However, there are examples of axial loading paths, e.g.
From (5.1), we infer D − < 0; hence, the second term in the right-hand side of (5.6) grows exponentially along V − . The asymptotic behaviour of the first term depends on the sign of E − : it grows to plus or minus infinity if E − < 0 and decays to a finite number if E − = 0 and to zero if E − > 0. The latter case corresponds to the proportional contractant deformation such that: This condition is equal to and describes the asymptotically stable stress path σ (s) in (5.6) attracting the direction of V − as s → ∞ by the mean of descent distance: We remark that the decay of the distance in (7.3) holds only for sufficiently large values of s. For instance, the distance may increase in some interval s ∈ (0, s 0 ) with s 0 > 0 satisfying cE − s 0 + h − (s 0 ) = 0, when cE − < 0 and h − (s 0 ) > 0 for C − < 0 in (7.1).

the left plot) obtain
Therefore, the stable asymptotic behaviour under the proportional dilatant deformation is guaranteed by D + > 0 and E + > 0, that is: or, equivalently, then the stress path σ (s) in (5.6) decays asymptotically to zero: By this, if E + > D + , then the stress path is closer to V + ; otherwise, it is closer to σ 0 − tr(σ 0 )V + /tr(V + ) when the opposite inequality E + < D + holds.
In the special case of tr(V + ) = 0, that is, a = 1/tr(U + ), from (19) we have The asymptotic stability should be considered separately.
If tr(V + ) > 0, then −ca tr(V + ) > 0 and follows the limit in (5.2) (see Fig. 2 the right plot): . (10.1) In this case, exp(cE + s + h + (s)) and exp(cD + s + h + (s)) entering the stress path σ (s) in (5.6) may obey various asymptotic behaviour in dependence of the sign of parameter C + in (10.1) and the values of E + , D + in (5.1). The asymptotically stable decay to zero of stress paths is described by the following two cases: either C + = 0 (hence, h + (s) ≡ 0) and E + > 0, D + > 0, that is: or h + (s) → −∞ provided by C + > 0: . The non-monotone behaviour of σ (s) is admissible for dilatant, too.
(iii) Volume-preserving deformation Consider a proportional deformation U * with tr(U * ) = 0; then, tr(V * ) = −1, E = a, and D = 0 in (5.1); hence, with h(s) = (σ 0 : U * )/(a tr(σ 0 )) + 1 (exp(cas) − 1), and the limit: We note that the limit states σ ∞ in (11.1) for all admissible directions U * with tr(U * ) = 0 and initial stresses σ 0 form the cone of the limit stress states with the apex in the origin of the space of principal stress axes. From (11.1), we find that σ * ∞ = aU * tr(σ ∞ ); thus, the interaction of this cone with the deviator plane, i.e. the plane with tr(σ ) = 1 (see [4]), is represented by the circle of radius a: The analysis of asymptotic behaviour of stress paths under various proportional deformations from Sect. 4 in dependence of the parameters is summarized for convenience in Table 1.  To prove this, we considered U = 1; therefore, tr(U ) becomes extreme if U i j = 0, i = j, so only the entries on the diagonal of U are distinct from zero. With the arithmetic quadratic mean inequality, we get: and equality holds for U 11 = U 22 = U 33 = 1/ √ 3. As we see in Table 1, the relations of D and E, one to each other and to zero, are crucial for the asymptotic behaviour of σ (s) as s → ∞. Depending on tr(U ) (x-axis) and a (y-axis), these relations are plotted in Fig. 3.
Here, we consider the formal description of the plot, which means the equations of the curves separating the different greyscale areas, starting in the down left corner: between area 1 and area 2: , tr(U ) ≤ 0; between area 2 and area 3: between area 3 and area 4: , tr(U ) ≥ 0; between area 4 and area 5: , tr(U ) ≥ 0; between area 5 and area 6: The meaning of the corresponding areas is the following: in the light grey areas 2 and 4 (where parameters D < 0 < E), stress paths are going to infinity attracting V ; in the medium grey areas 1 and 5 (where D < E < 0), stress paths do not have any asymptote and diverge; in the dark grey area 3 (where 0 < D < E) stress paths go to zero attracting V ; in the white area 6 (where 0 < E < D), the behaviour of the stress paths depends also on the initial stress state σ 0 , while in the other areas the asymptotic behaviour is independent of the initial stress state.
Further we illustrate the usage of Table 1 for the cases of contractant, dilatant, and volume-preserving deformations.

(i) Contractant straining
If tr(U ) < 0, the stress path is asymptotically stable in area 2. We get the lower bound for tr(U ) = − √ 3, which means isotropic compression, by For a > a min , we can be sure to have asymptotically stable behaviour for arbitrary proportional deformation (see Fig. 1 from Sect. 3.1).

(ii) Dilatant straining
If tr(U ) > 0, the stress path is asymptotically stable in areas 3 and 4: in area 3, the path tends to zero, while in area 4, it tends to infinity. We get the upper bound for tr(U ) = √ 3, which means isotropic extension, by For a < a max , we have asymptotically stable behaviour for arbitrary proportional extension; in this case, the stress path will also tend to zero.

(iii) Volume-preserving deformation
For the limiting case of tr(U ) = 0, the stress path will tend to a certain stress state at the critical cone. To illustrate the matter, we plot few results of numerical simulation.
In the left plot of Fig. 4, the stress paths are projected in the plane spanned by −σ 33 and the first median of the − √ 2σ 11 = − √ 2σ 22 -plane. The dot dashed line is the isotropic axis; the outer solid lines are the intersection with the critical cone. In the right plot, the paths are projected onto the deviator stress planeσ * = σ /tr(σ )− I /3. In the deviator plane, the cutting line with the critical cone forms a circle of radius a according to (11.2). In Fig. 4, we see three example stress paths under volume-preserving deformation U = diag(1, 1, −2)/ √ 6 with the parameter a = 0.35 and three initial stress states We see that for different deformations U we get different behaviour of the stress paths. For U 1 (solid curve), we have isotropic compression, and the corresponding stress path tends to the isotropic axis. For U 2 (dashed curve), we have uniaxial compression, and we can also observe an asymptotic behaviour. For U 3 (dotted curve), we have volume-preserving deformation and the stress path goes to the critical cone.

Case studies
In this Section, we will discuss the case of isotropic deformation and uniaxial deformation in more detail.   (9), and (10.2), and from (4.3), we conclude unconditional asymptotic decay to zero: By this, if E + > D + , which holds for a < 1/ √ 3, then σ (s) attracts closer the direction of V + . Otherwise, (σ 0 ) * is more attractive. The result in dependence on the parameter a is gathered for convenience in Table 2.
We remark that this particular case of isotropic compression/extension implies μ = D ∓ being the single eigenvalue and μ = E ∓ the double eigenvalue for the following eigenvalue problem written with respect to μ, respectively: where 1 stands for the 3-by-3 matrix of ones and V ∓ is the eigenvector corresponding to D ∓ ; see [32] for details.

Conclusions
In the paper, for a constitutive model based on the concept of hypoplasticity of the Kolymbas type the stress paths obtained under proportional deformations are investigated. In particular, a simplified hypoplastic constitutive equation is considered where the objective stress rate is a function of the current stress and strain rate. The model is obtained by omitting the influence of the change in the pressure-dependent relative density in the hypoplastic model originally proposed by Bauer and Gudehus. For arbitrary proportional deformations starting from the stress-free state, the corresponding stress paths are linear which is in accordance with the first law by Goldscheider and also an important property for constitutive models relevant to frictional granular materials. From experiments, it is known that significant changes in the density of the granular material can lead to a slightly curved stress path which can also be influenced by grain crushing under higher stresses. Such properties can be taken into account using more enhanced hypoplastic models, which are not considered in the present paper.
The hypoplastic equation considered only includes two constitutive constants: a stiffness parameter and a limit stress state parameter. The former can be related to an isotropic compression test and the latter to the so-called critical friction angle defined in the steady state of a cohesionless granular material under triaxial compression. Although the influence of the direction of the stress deviator on the limit stress state parameter is neglected, it does not mean a restriction of the general results drawn. A relevant relation between the values for the direction of the stress deviator of triaxial compression and the one for another direction can, for instance, be obtained by interpolation.
For the hypoplastic constitutive equation considered, the analytical solution of the stress paths depending on monotonous linear deformation paths is obtained in a closed form. On its basis, the course of the stress paths is examined in dependence of the material parameters and prescribed proportional contractant, dilatant, and volume-preserving deformations. It is shown that stress paths starting from arbitrary initial stress states are usually nonlinear, and their convergence to a proportional stress path (second law by Goldscheider) is asymptotically stable only for a certain domain of the limit stress state parameter, which is related to the critical friction angle.