Stress-controlled ratchetting in hypoplasticity: a study of periodically proportional loading cycles

We investigate rate-independent strain paths in a granular material generated by periodically oscillating stress cycles using a particular constitutive model within the hypoplasticity theory of Kolymbas type. It is assumed that the irreversible hypoplastic effects decay to zero when the void ratio reaches its theoretical minimum, while the void ratio is in turn related to the evolution of the volumetric strain through the mass conservation principle. We show that under natural assumptions on material parameters, both isotropic and anisotropic stress cycles are described by a differential equation whose solution converges asymptotically to a limiting periodic process taking place in the shakedown state when the number of loading cycles tends to infinity. Furthermore, an estimation of how fast, in terms of the number of cycles, the system approaches the limit state is derived in explicit form. It is shown how it depends on the parameters of the model, on the initial void ratio, and on the prescribed stress interval.


Introduction
This is a continuation of the paper [12], where we studied some mathematical consequences of the constitutive stress-strain relation for hypoplastic granular materials like cohesionless soil or broken rock. The constitutive law is of the rate type, incrementally non-linear, and it is based on the hypoplastic concept proposed by Kolymbas [35]. Compared to other types of nonlinear constitutive equations, e.g. for hyperelastic and hypoelastic material laws, the hypoplastic constitutive relations are not differentiable at zero strain rate. This permits the constitutive modelling of 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, and there is no need to introduce additional functions to model the transition from elasticity to plasticity on monotone loading or unloading paths. Detailed discussions about physical aspects of hypoplastic models can be found for instance in [25,37,43], their response to cycling loading is studied in [14,38,47], shear localization in [4,11,15,22,29], critical states in [5,7,53], behaviour under undrained condition in [44], creep and stress relaxation in [6,9,40,46], grain fragmentation in [8,13], extension to a micropolar continuum in [29,30], and thermodynamic aspects in [32,48,49].
Other typical representatives for incrementally non-linear constitutive equations are for instance the Armstrong-Frederick model [2] as a hypoplastic extension of the von Mises model for plasticity, the endochronic model by Valanis [51], the octolinear model by Darve [20,21], the CLoE model by Chambon et al. [16], and the barodesy model by Kolymbas [36]. Variational approaches to modelling granular and multiphase media can be found for instance in [1,33,39].
There is a large number of investigations devoted to the properties of the stress paths under a monotonic proportional deformation path. Goldscheider [23] formulated two rules saying (a) that a proportional deformation path starting from a nearly stress-free state results in a proportional stress path, and (b) 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 are an intrinsic feature of frictional granular materials also confirmed by further experiments, e.g., [19,50], and numerical simulations using the discrete element method, e.g., [31,41]. This feature can be interpreted by fading memory of the material on the initial state [26]. The asymptotic behaviour of the stress path for proportional deformation has a great importance for constitutive modelling of frictional granular materials. This is, for instance, the basic concept of the framework of barodesy [45] and it is also an intrinsic feature in hypoplasticity [25,27,34,42,43,46]. In hypoplasticity, asymptotic properties were numerically investigated by Kolymbas [34] and a mathematical criterion was proposed first by Niemunis [46].
In [12], we have proposed a different strategy by considering the strain-stress law under a given proportional strain trajectory as a nonlinear differential equation for the stress. To this end, a simplified version of the hypoplastic model by Gudehus [24] and Bauer [3] was considered. For the sake of simplicity the pressure dependent properties of granular materials described in [3,24] were omitted, so that only six material parameters remain in the present model. Note that the same simplified version can also be obtained from the hypoplastic model by von Wolffersdorff [52] as shown in [10]. The solution of the corresponding non-linear problem was found in closed form. In this way we made also a close link to barodesy models [36]. The exact solution provided in [12] made it possible to describe analytically various scenarios of the behaviour of stress paths obtained from monotonic proportional 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 have identified the domain of the constitutive parameters which guarantee that for proportional deformation paths the corresponding stress path starting from an arbitrary initial stress state are asymptotically stable.
The goal of the present paper is to study the phenomenon of ratchetting (that is, the shift of the strain cycles under periodic loading and unloading stress cycles) resulting from the simplified hypoplastic model interpreted as a stress-controlled element test. In its original simplified form [12], the theoretical ratchetting would be infinite, which is not in agreement with experimental observations. This is due to the fact that the influence of the change of the void ratio of the granular material on the incremental stiffness is neglected there. Similarly to the hypoplastic version proposed in [3,24], we therefore assume here that the non-linear part of the constitutive law contains a multiplier depending on the current void ratio with the property that it vanishes if the void ratio reaches the most condensed admissible state. Herein, the current void ratio, which is defined as the ratio of the volume of the void space to the volume of the solid material, is considered as a state quantity in the constitutive model. In addition, the stress rate obtained from the constitutive equation is scaled by a factor which, for densification of the material, tends to infinity and thus the material becomes asymptotically rigid.
For granular materials with incompressible grains the evolution of the void ratio itself can be expressed as a function of the strain rate by virtue of the mass conservation principle. Since we are in the stress-controlled case, the strain rate is the unknown of the problem and has to be found as a solution of a system of differential equations. Unlike in the strain-controlled case, the solution is not found in closed form and has to be computed numerically. On the other hand, we prove that the strain paths subjected to periodic loading-unloading stress cycles converge to an equilibrium independently of whether the proportional stress paths are isotropic or not.
The paper is organized as follows. In Sect. 1, we present the constitutive model used and derive the differential equations which are studied in the sequel. The case of isotropic stress cycles is investigated in detail in Sect. 2. The convergence of the strain trajectory toward the equilibrium when the number of the loadingunloading cycles tends to infinity as well as an estimate of the speed of convergence are proved in Sect. 3. Sect. 4 is devoted to the case of anisotropic stress cycles and we show that the differential equation describing the process has the same asymptotical properties observed in the isotropic case. Some numerical tests are shown in Sect. 5.

The constitutive model
In what follows we study the extension of the constitutive model from [12], accounting only for coaxial deformations on the element level, i.e. the case when the direction of principal stresses coincides with the direction of the corresponding principal strains. In this case, the objective stress rate equals the material time derivative of the Cauchy stress, herein denoted byσ (t). Then the constitutive model that we consider to describe the relation betweenσ (t) and the strain rate tensorε(t) is governed by the following equatioṅ where c 1 (t) < 0 and f (t) ≥ 0 are given in terms of the current void ratio of the granular material in a sense that will be made more precise below. We further assume that L(σ ) is a fourth order tensor and N (σ ) is a second order tensor of the form under the condition that the trace tr σ of σ does not vanish. We denote by I 4 the unit fourth order tensor, and a is assumed to be a material constant. The incrementally non-linear character of the constitutive equation (1.1) ensures the modelling of inelastic and dissipative material properties. For the case of coaxial deformations, the constitutive equation (1.1) is not restricted to problems of small deformations. Besides, the constitutive equation (1.1) is rate independent and, consequently, invariant with respect to an arbitrary time transformation. Therefore, the variable t need not be necessarily interpreted here as physical time, but rather as an abstract monotonically increasing parameter within an interval t ∈ [b, b + T ]. In our case, as we shall see in (2.6), it will be measured in terms of the stress variation. For our purposes we denote by ·, · the canonical scalar product ":" in R 3×3 and, with the help of (1.2)-(1.4), the evolution equation (1.1) can be written in the forṁ where, in accordance to [5], f is given as a function of the void ratio e(t) by the function c 1 is assumed to depend on e(t) in the form Moreover, for granular materials with incompressible particle density the evolution of the void ratio can be computed from the mass balance equatioṅ In the density function (1.6) the current value of the void ratio, e(t), is related to the critical void ratio e c and minimum void ratio e d in such a way that for e(t) = e d the value of function f (e(t)) becomes zero. The parameters e c > e d > 0 and α > 0 in (1.6) are assumed to be constant. The idea is to derive from the identities (1.5)-(1.9) a differential equation for the unknown function e(t). We need to check that we have e(t) > e d during the whole process to make formulas (1.6) and (1.7) meaningful. Besides, we assume the constantsc > 0 and β > 1, meaning that the material becomes rigid when e(t) asymptotically converges to e d .
The stress-controlled proportional loading consists in choosing a fixed tensor satisfying the condition tr S < 0, and assuming σ (t) in the form . This is what we call a proportional loading if σ is increasing, and proportional unloading if σ is decreasing. Let us note that for a cohesionless granular material, only negative principal stresses are physically admissible, i.e., all elements in the diagonal of the tensor S must be negative. These assumptions ensure that we haveσ 11 < 0, σ 22 < 0,σ 33 < 0 for the compression phase andσ 11 > 0,σ 22 > 0,σ 33 > 0 for the extension phase. In view of (1.11), we rewrite Eq. (1.5) aṡ In Sect. 2 we first consider the isotropic case, while a more general case will be studied in Sect. 4.

The case of isotropic stress cycles
Assume for the moment that the condition The isotropic case is characterized by the choice S = −I; therefore, equation (1.12) corresponds toσ Taking the scalar product with I we get wherefrom we obtain Inserting the identity above in (2.2) we can see thaṫ Recalling that c 1 (t) < 0, for the loading asσ (t) > 0, from (2.3) together with (2.4) and (2.1) it follows: whereas for the unloading asσ (t) < 0 respectively Therefore, by using (1.9) for ε(t), I we obtain the following differential equation for e(t) where sign denotes the sign function.
Let us consider now a periodically oscillating function σ between some fixed values 0 < σ < σ . More specifically, we put T = log(σ ) − log(σ ), t k = kT for k = 0, 1, 2, . . . , and define Note that with this choice of a periodic function the right-hand side of (2.5) no longer depends on t which is rather advantageous for the analysis to be carried out in the sequel. Indeed, for σ as in (2.6) we have (2.10) or, equivalently, putting problem (2.9) and (2.10) is of the form (2.12)

Convergence
We first note that, since h in (2.7) and g in (2.8) are positive for u > 0, both H and G are increasing in (0, ∞), We shall see below that the convergence is guaranteed as long as the function H − G is increasing. This is, however, not the case on the whole half-line (0, ∞). Indeed, we have which is positive if and only if (2.11) and (2.12) we get the recurrent relations for j = 0, 1, 2, . . .
with an initial condition u 0 > 0. We prove the following result.

Summing (3.3) and (3.4) we obtain
Recalling that H + G is increasing and lim u→0+ (H + G)(u) = −∞, we must have 0 < u 2 j+2 < u 2 j . The sequence {u 2 j } is thus decreasing and contained in the interval (0, U 0 ). There exists thereforeū ∈ [0, U 0 ) such that lim j→∞ u 2 j =ū. We now rewrite (3.4) as and adding this identity to (3.3) yields by virtue of (3.5) that Since H − G is increasing in (0, U 0 ), we have 0 < u 2 j+1 < u 2 j−1 , in other words, the sequence {u 2 j+1 } is decreasing. Hence, using also (3.5), we see that it converges to some 0 ≤ u ≤ū. Assume first that u > 0. We can pass to the limit in (3.3) and (3.4) and obtain Recall that the system (3.3)-(3.4) has been derived under the condition (2.1) provided u 0 < U 0 . It remains to check a posteriori that its solution satisfies (2.1). Indeed, we have 0 < u(t) < U 0 during the whole evolution, so that by (1.6) and (3.2) we have which is precisely (2.1). The condition u 0 < U 0 is not a real restriction. It was shown for example in [12] that physically meaningful values of the parameter a are 1/(2 In terms of the initial void ratio e(0), this can be interpreted by virtue of (3.2) as (3.10) Typical values of α for sands lie within the range of 0.1 ≤ α ≤ 0.3, e.g. [3,28]. Therefore, in practice, the factor on right-hand side of (3.10) can never be exceeded and does not impose a real restriction in the modelling of materials like sand. For other granular materials the relevant range for values of α may be different depending on particle shape and particle size distribution.
Summing up the two equations and using (3.12) we get where the last equality holds true for some u ∈ [u 2 j−1 , u 2 j ] thanks to the Mean Value Theorem. By (2.7), (2.8), and (3.2) we have (3.14) We can restrict ourselves to the case u < U 0 and obtain from (3.14) the upper bound Let us estimate the decay rate of H (u) as u → 0+. By (2.7) we have with q =c(1 + 3a 2 ), r = 1 + e d . For x ∈ (0, U 0 ) we thus obtain The two integrals in (3.16) have to be treated separately. For u ∈ (0, U 0 ), we have r + u < r + U 0 , while for u ≥ U 0 we use the formula r + u ≤ ( r U 0 + 1)u. We thus get In particular, for We now prove for x > 0 and ξ < 0 the implication .

(3.23)
Then (H − G)(x) > ξ, which contradicts the implication (3.21). Using (3.22) we obtain from (3.20) the inequality The function (ξ) := 1 1 + ρξ 1+ρ is increasing and convex on (0, ∞) thanks to the assumption that β > 1. Hence, (ξ j+1 ) − (ξ j ) ≥ (ξ j )(ξ j+1 −ξ j ), that is, for every j ∈ N. We conclude that for all j ∈ N with some constantĈ > 0. Using (3.19) and (3.26) we thus get On the other hand, from (3.15) we deduce the following lower bound for H . (3.27) Hence, using the fact that G(u 2 j+1 ) < G(U 0 ) which is a constant independent of j we obtain that (note that β > 1 by hypothesis)μ For j ∈ N sufficiently large, we thus deduce the convergence behaviour with a constant C > 0 independent of j. Note that the constant C is related to all material parameters involved in the constitutive equation as well as to the initial void ratio and to the prescribed cyclic stress interval. We see that the number of cycles necessary to reach a given neighborhood of the shakedown state gets large for small values of α and large values of β, see Figs. 5 and 6, respectively.

The case of anisotropic stress cycles
In this section we consider the case when the tensor S is not necessarily a multiple of the unit tensor, i.e., S 11 , S 22 and S 33 can possibly have different (negative) values. This is what characterizes the anisotropy in the stress cycles. Our goal is to show that even in the anisotropic case the problem can be reduced to an equation like (2.5) which exhibits similar asymptotically periodic character as in the isotropic case. The main difference is that in order to obtain an expression for the norm of the strain rateε some further intermediate steps are needed. Assumingσ (t) = 0, we introduce the following notation In the case of loading, that is,σ (t) > 0, (1.12) can be rewritten as (note that c 1 (t) < 0) Similarly, for unloadingσ (t) < 0 we get We rewrite (4.2) and (4.3) in a more convenient form. Taking the scalar product of (4.2) and (4.3) with Q and noting that Q, I = 1 yields so that (4.2) and (4.3) can be reformulated as we can simplify (4.5) and get In summary, the problem consists in solving the ODEs Note that by (4.1) and (4.6) we have Q 2 ≥ 1/3 and The counterpart of condition (2.1) reads here Assuming that (4.9) holds we get positiveness for the discriminant of the quadratic equation following from (4.8), and its positive roots are: Hence, The reader is invited to verify that in the isotropic case, i.e. S = −I, we have Q = 1 3 I, A = B = 1 1+3a 2 I, and consequently the identity above yields equation (2.5) studied in Sect. 2. We now rewrite (4.11) in a more concise form with (4.14) Then we have We now use (4.1) and (1.9) to derive the identity Recalling that σ (t) σ (t) = 1 forσ (t) > 0, and σ (t) σ (t) = −1 forσ (t) < 0 and using also (1.6), we obtain differential equations for the unknown function e(t) in the forṁ

Numerical tests
The focus of the numerical simulations presented in this section is on the influence of some prescribed quantities on the evolution of the volume change with increasing number of loading cycles. More precisely, the influence of the given non-zero components of the stress direction S 11 , S 22 , S 33 , the initial void ratio e(0), and the value of parameters α and β. The material parameters involved in the present constitutive model, i.e. {a, e d , e c , α,c, β}, have the following meaning: parameter a is related to the intergranular friction angle in a stationary state or the so-called critical state [5,7,10]; e d and e c denote the minimum void ratio and critical void ratio, respectively; parameter α triggers the amount of the density function f (e) depending on the current void ratio e; parameter c scales the amount of the incremental stiffness (i.e. a higher value ofc results in a higher stiffness); and parameter β > 1 triggers the influence of the state dependent quantity u = e − e d on the stiffness factor c.
More details about the physical meaning of the material parameters and the procedure for calibration can be found for instance in [3,10]. In connection to that, it is worth recalling that for highly non-linear constitutive equations, like the one we investigate here, certain restrictions on the possible range of the parameter values may apply. For the present numerical simulations the material parameters are adapted to experimental data typical for sands, i.e. a = 0.4, e c = 0.8, e d = 0.4,c = 2, α = 0.1, and assuming β = 1.03. For the interval of cyclic loading along the given stress direction, the interval limits of the scalar stress factor in (1.11) are σ = 10 and σ = 12.
For prescribed isotropic cyclic stress states with S = − I, starting from initial conditions at t 0 = 0 chosen as e(0) = 0.7, σ = 10, ε(0) = 0 · I, and for n = 25 cycles the results of the numerical simulation are shown in Fig. 1. From Fig. 1a it is clearly visible that for a full loading-unloading-reloading cycle the compaction, i.e. the decrease of the volumetric strain ε v = ε 11 + ε 22 + ε 33 , is more pronounced at the beginning of cyclic loading and becomes smaller with number of stress cycles. This behaviour is in agreement with experimental observations and can be explained by an increase of the incremental stiffness caused by a reduction of the void ratio e. As the number of cycles increases, the void ratio tends towards the minimum void ratio e d , i.e. the state quantity u = e − e d asymptotically approaches zero, as shown in Fig. 1b. It can be noted that experiments with stiff grains usually show after reloading a certain loop of the the stress-strain curve. This is caused by a higher stiffness at the beginning of the reloading path, which is not captured by the present model, but it can be modelled with additional state quantities using more refined hypoplastic material descriptions, e.g. [14,47]. However, the essential shakedown feature of granular materials can already be simulated with the present simplified model.
The behaviour under axisymmetric cyclic stress states, i.e. S 11 = −3, S 22 = −1, S 33 = −1, is shown in Fig. 2. From Fig. 2a it can be seen that the components ε 22 = ε 33 have a tendency to increase, while the component ε 11 decreases. Within a full stress cycle, the accumulated volumetric strain exhibits compaction, which means a decrease of the void ratio towards the minimum void ratio as shown in Fig. 2b. Numerical simulations show that the sign of the evolution of the individual strain components strongly depends on the   The results of numerical simulation in the case of triaxial anisotropy with the choice S 11 = −1.5, S 22 = −0.5, S 33 = −1 is shown in Fig. 3. From Fig. 3a it can be seen that the components ε 11 and ε 33 have a tendency to decrease, while the component ε 22 increases. Again the volumetric strain exhibits compaction, which is also related to an decrease of the void ratio towards the minimum void ratio as shown in Fig. 3b. It cannot be expected, however, that in this kind of problems numerical tests would give any useful validation of the convergence rate computed theoretically in Sect. 3. Note that by formula (3.28) the convergence is very slow, therefore any numerical analysis in this matter would require a large number of steps subjected to accumulated numerical error.
Further numerical simulations show that for a full loading-unloading-reloading cycle the changes in the volume depend not only on the number of loading cycles, but also on the material parameters involved in the constitutive model, the given stress direction, the initial void ratio and to the prescribed stress interval. For example, the influence of the initial void ratio e(t = 0) = e(0) on the change of the density factor f (e) is demonstrated in Fig. 4. In particular, Fig. 4a shows the behaviour under isotropic cyclic stress states with S = − I, and Fig. 4b shows the behaviour under axisymmetric anisotropic cyclic stress states corresponding to S with the components S 11 = −3, S 22 = −1, S 33 = −1. As the density factor f (e) is a function of the current void ratio, a reduction of f (e) means a densification. It is clearly visible that in both cases the cyclic  behaviour shows additional compaction with number of cycles. The compaction is more pronounced for the initially higher void ratio, in particular, at the beginning of cyclic loading. After same number of cycles the compaction is higher under isotropic cyclic stress states than for the case of axisymmetric anisotropic cyclic stress states.
For an initial void ratio of e(0) = 0.7 the influence of parameter α on the evolution of the density function f (e) is shown in Fig. 5a for the case of isotropic cyclic stress states and in Fig. 5b for the case of axisymmetric anisotropic cyclic stress states (using, in each case, the same components for the tensors S as in Fig. 4a and Fig. 4b, respectively). In both cases the evolution of f (e) shows again compaction with an increase in the number of cycles. For same stress cycles and an initial void ratio of e(0) = 0.7 the influence of parameter β on the evolution of the density function f (e) is shown in Fig. 6. For a higher value of β the compaction is lower which can be explained by a higher incremental stiffness, i.e., the value of c 1 obtained from (1.7) scales the incremental stiffness and it is higher for a higher value of β.
In accordance with the analytical results, the numerical simulations also show that the cyclic behaviour under the particular conditions (1.10) and (1.11) leads to an additional compaction of the material in each loading-unloading-reloading cycle. Independently of the given stress direction and the initial void ratio, the amount of the additional compaction decreases with number of cycles. In this context, it should be noted that in the present constitutive model compaction is not an inherent property of formula (1.6) for the density factor f (e). Indeed, it is shown in [3] that an increase of the volume can be observed for the same choice of density function (1.6) under stress cycles different from (1.10) and (1.11).

Conclusion
We study a special case of the extended hypoplastic constitutive model for granular materials subjected to proportional periodically oscillating stresses under the assumption that the evolution of the strain components depends on the prescribed stress interval and the current void ratio, which is defined as the ratio of the volume of the void space to the volume of the solid material. For a prescribed stress cycle a part of the corresponding deformations is irreversible, indicating to the dissipative property of the incrementally non-linear constitutive model used. In the analysis it is proved that within each full stress cycle the resulting volumetric strain exhibits a compaction, i.e., a reduction of the void ratio. The amount of additional compaction reduces non-linearly with the number of cycles in both isotropic and anisotropic loading scenarios. Under technical assumptions on the parameters of the model, we find an upper bound for the initial void ratio which guarantees that the current void ratio converges asymptotically to its minimum value, which corresponds to the full compaction of the material. Numerical results indicate that the ratchetting trajectories very strongly depend on the degree of stress anisotropy. This question, however, deserves to be studied in more detail in a future work.