# Fluctuations and failure in granular materials: theory and numerical simulations

Luigi La Ragione, Giuseppina Recchia, Felix Darve, Francois Nicot, Antoine Wautier

► To cite this version:

Luigi La Ragione, Giuseppina Recchia, Felix Darve, Francois Nicot, Antoine Wautier. Fluctuations and failure in granular materials: theory and numerical simulations. Granular Matter, 2024, 26 (3), pp.68. 10.1007/s10035-024-01431-5. http://dx.doi.org/10.4179/2-02-001431-5

HAL Id: hal-04689851

https://hal.inrae.fr/hal-04689851v1

Submitted on 10 Sep 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Fluctuations and failure in granular materials: theory and numerical simulations

Luigi La Ragione<fnref n="1" />* · Giuseppina Recchia<fnref n="1" /> · Felix Darve<fnref n="2" /> · Francois Nicot<fnref n="3" /> · Antoine Wautier<fnref n="4" />

Received: 15 January 2024 / Accepted: 31 March 2024 / Published online: 4 June 2024
© The Author(s) 2024

## Abstract

We consider a dense aggregate of elastic, frictional particles isotropically compressed and next uniaxial strained at constant pressure. We show how failure can be predicted if fluctuations in the kinematics of contacting particles are introduced. We focus on the second order work and the possibility that at some stressed states it becomes negative under proper perturbations. Our analysis involves both a theoretical model and numerical simulations based upon the distinct element method (DEM). The theoretical model deals with contacting particles with incremental relative displacements that deviate from the average deformation in order to ensure their equilibrium. Because of this, the macroscopic stiffness tensor of the aggregate, that relates increments in stress with increments in strain, does not have the major symmetry. Consequently, in the hardening regime, we predict stressed states in which the second order work vanishes. The model seems transparent, and it makes clear and illustrative the role played by the fluctuations introduced in the kinematics of contacting particles in relation to the vanishing of second order work in an aggregate of compressed particles. The comparison with numerical simulations data supports the model.

**Keywords** Granular materials · Micromechanics · Fluctuations · Discrete element method · Failure

## 1 Introduction

Among complex systems, granular materials have attracted much interest in the past few years in materials sciences. Granular materials can be encountered in a variety of engineering contexts, such as pharmaceutical engineering, food

particle storage, geological mass-driven hazards, and civil engineering. It is well recognized that granular materials can exhibit a wide spectrum of emergent properties that have received much interest [1–6]. Such properties appear at the specimen scale, whereas they cannot be observed at smaller scales. As an illustration, granular materials are dissipative structures [7, 8]. Thanks to their disordered microstructure, such materials can adapt themselves through microstructural reorganizations to resist against a given external loading on the specimen scale. Such reorganizations mostly entail dissipative mechanisms through grain sliding and rolling at the inter-granular contacts. They act as elementary, microscopic events that are responsible for the whole system to be in mechanical equilibrium. It is worth emphasizing that the notion of mechanical equilibrium should not be confused with that of thermodynamic equilibrium. Indeed, thermodynamic equilibrium means that any fluctuation that takes place stays bounded and eventually vanishes. On the other hand, mechanical equilibrium imposes that the external tractions applied on the boundary of the system are balanced by the internal stress field developing within the system. When a system is in thermodynamic equilibrium, it is automatically in mechanical equilibrium which is a less restrictive

<figure><img src="image_2.png" /></figure>

Luigi La Ragione
luigi.laragione@poliba.it

Felix Darve
felix.darve1@univ-grenoble-alpes.fr

Francois Nicot
francois.nicot@univ-smb.fr

Antoine Wautier
antoine.wautier@inrae.fr

<footnote n="1">¹ Politecnico di Bari, Via Edoardo Orabona, 4, 70126 Bari,
BA, Italy</footnote>

<footnote n="2">² University Grenoble Alpes, Grenoble INP, CNRS, Lab 3SR,
Grenoble, France</footnote>

<footnote n="3">³ Université Savoie Mont Blanc, IISTerre, Chambéry, France</footnote>

<footnote n="4">⁴ INRAE, Aix-Marseille University, UMR RECOVER,
Marseille, France</footnote>
equilibrium type. Thus, even at mechanical equilibrium, a given material system can be out of thermodynamic equilibrium. In that case, some fluctuations are expected to develop within the system, making it evolve in an irreversible way from a given state toward another state. Irreversibility means that the system cannot return spontaneously to the previous states without any external exchange, mainly because of internal dissipative processes. Fluctuations at the micro scale, instead of vanishing, can be amplified. In some circumstances, the system can then bifurcate towards a new response mode, with a new topological organization of the energy dissipation. This is what happens on dense granular assemblies under a shear loading (as a drained triaxial loading, for example). It is repeatedly observed that the deviatoric stress reaches a peak value after a monotonous increase (so-called hardening regime), followed by a decrease (so-called softening regime) until the critical state regime takes place. Furthermore, the specimen is no longer homogeneous once the deviatoric stress peak is reached, as a localization pattern (with a single or multiple shear bands) develops. It was recently shown that the development of a shear band pattern during the softening regime corresponds to a phase transition, marked by a break in symmetry in the kinematic field, consisting in a microstructural reorganization to optimize the inter-granular, frictional dissipation processes (optimal dissipative structure), and make the specimen able to resist against the external loading. The appearance of shear bands within granular materials, patterning of the material in several structured zones, has received much attention over the past decades. It is worth mentioning that such features can be observed irrespective of the scale considered: along faults during earthquakes, in large rocky cliffs, or in laboratory test specimens. Focusing on the lab scale, the occurrence of such events is well predicted by the Rice-Mandel criterion [9, 10]. This criterion was shown to be a subset of the second-order work criterion introduced in the middle of the past century by Hill [11]. More recently, the second-order work theory was formalized [12], by relating the occurrence of an outburst in kinetic energy, regarded as a quasi-unbounded fluctuation, and the vanishing of the second-order work. This vanishing corresponds to the loss of the positive definiteness of the symmetric part of the constitutive operator relating both incremental stress and strain operating on the material point scale. Along a given loading path, all mechanical states corresponding to the non-ellipticity of the symmetric part of the constitutive operator thus define the bifurcation domain. In such a domain, an effective failure of the material can occur according to the loading direction applied, and to the loading control adopted [13].

In Nicot and Darve [14], two material failure modes in granular material are presented: one associated with

localization and the other associated with a diffuse mode. Localization has been clarified through the contributions made by Rudnicki and Rice [10] and Vardulakis [15, 16] in which the vanishing of the determinant of the acoustic tensor occurs when a shear band develops. They use an elasto-plastic model with a non-associative plastic flow. More recently, La Ragione et al. [17], based upon a micro-mechanical approach, predict shear band in stressed granular material with a macroscopic stiffness tensor that does not have the major symmetry, consequently with a similar structure of that introduced in [10, 18]. The other failure is associated with a diffusive mode that can anticipate the localization band when the second order work vanishes (e.g. [19]).

Here we focus on the second order work and the possibility that at some stressed states, during the loading path, it becomes negative under proper perturbations [20, 21]. Our analysis involves both a theoretical model and numerical simulations. The merit of the simulations, based upon the distinct element method (DEM), is to reproduce the proper kinematics and statics of every single grain of the sample while both force and moment equilibrium are satisfied at each applied incremental macroscopic strain [22]. In so doing, it is possible to follow the behavior at contact level of the aggregate, including particle deletion, sliding and rolling.

The novelty of our work is the introduction of fluctuations in particles deformation that allows to predict a stressed state in which the second order work vanishes. Following the DEM analysis, it is possible to employ theoretical models that reproduce the behavior of the aggregate as seen in the simulations. The simplest hypothesis is that contacting particles move according with an average deformation. The resulting macroscopic response of the aggregate seems reasonable from the qualitative point of view. Particle sliding and deletion are predicted as well as the deviatoric stress and the volume strain along a monotonic triaxial test (e.g. [23, 24]). On the other hand, the comparison with numerical simulation shows a quantitative difference that is solved by relaxing the kinematics with the introduction of a fluctuation field along with equilibrium [25]. The kinematic of contacting particles with fluctuations also leads to a qualitative difference with respect to the average strain, under particular loading condition: the macroscopic stiffness of the aggregate, $A_{ijkm}$, that relates the increments in stress with the increment in strains, $\dot{\sigma}_{ij} = A_{ijkm} \dot{E}_{km}$, can be characterized by the loss of the major symmetry. This represents a necessary condition to predict localization.

Following this idea, we focus our investigation on a triaxial test conducted at constant pressure conditions. As
emphasized in [17], under specific conditions, the incremental response of the aggregate tends to be frictionless on average, resulting in a non-symmetric elasto-plastic stiffness tensor. Within this framework, we show that the second-order work linked with stressed states can vanish.

## 2 Triaxial loading

We focus on a triaxial test of an aggregate of N identical, frictional, elastic particles with diameter D which deform according to a macroscopic average strain, E. After an initial isotropic compression, $p_0$, the aggregate is uniaxial compressed, $E_{33} < 0$, keeping the pressure constant. We define the shear strain $\gamma = -1/2(E_{33} - E_{11})$, with $E_{11} = E_{22}$, and the volume strain, $\Theta = -(E_{33} + 2E_{11})$, positive in compression, while $\Theta_0$ is the volume strain associated with the pressure at the isotropic state. Particles interaction is given by a non central contact force with a normal component that follows the Hertz law and a tangential component that has an elasto-plastic response: a bilinear relation with an elastic displacement followed by a frictional sliding. The incremental force between two contacting particles is given by:

$$ \dot{F}_i = \dot{F}_i^N + \dot{F}_i^T \quad (1) $$

in which $\dot{F}_i^N$ and $\dot{F}_i^T$ are, respectively, the incremental normal and tangential component of the contact force. Given $\dot{F}$, the incremental average stress $\dot{\sigma}$ may be written as the average over all particles in a region of homogeneous deformation that is identified with the continuum point [25]:

$$ \dot{\sigma}_{ij} = \frac{1}{2} \frac{N}{V} \sum_{p=1}^{M^p} \frac{M^p}{N} \langle \dot{F}_i \rangle \dot{d}_j^p \quad (2) $$

where $M^p$ is the number of pairs of contacting particles with contact vectors (the vector that joins the particle centers) within the p-th element of R equal elements $\Delta\Omega^p$ of solid angle in which the surface of the unit sphere has been partitioned, V is the total volume of the aggregate and $\langle \dot{F}_i \rangle \dot{d}_j^p$ is the average of all the contact forces between pairs whose contact vectors, $d_i = D\hat{d}_i$, are within a given solid angle centered at $\dot{d}_j^p$. We have $\Delta\Omega^p = 4\pi/R$ and $\sum_{p=1}^R M^p = 2N_c$ in which $N_c$ is the total number of contacts in the aggregate.

In the monotonic triaxial test, the deviatoric stress is $q = -1/2(\sigma_{33} - \sigma_{11})$ and its components are $q_N$ and $q_T$; the former is the part of the deviatoric stress associated with the normal component of the contact force, the latter with the tangential component of the contact force. Therefore, the response of the aggregate is characterized by the $q-\gamma$ and the $\Theta-\gamma$ curves.

We investigate the incremental response of the aggregate, at given stressed state along the curve $q-\gamma$; in particular, we look at the second order work, $W_2 = \dot{\sigma}_{ij}\dot{E}_{ij}$, and the possibility that it may become negative. We do this through numerical simulations and a theoretical model.

## 2.1 Numerical simulation

We employ a numerical simulation based upon DEM on a sample made of 10.000 elastic particles with the shear modulus $G = 2.9 \times 10^9$ Pa, the Poisson ratio $\nu = 0.2$, and radii 9.5 and $10.5 \times 10^{-2}$ mm, initially isotropically compressed and then sheared. More detailed about the protocol adopted can be found in [26]. The initial compression is carried out without friction to generate a very dense state, volume fraction $\phi = 0.64$; next, the aggregate is uniaxial strained at constant pressure, $p_0 = 200$ KPa, with friction coefficient $\mu = 0.5$ and increments of deformation, at every cycle, of the order of $10^{-6}$. Along the loading path, there are intermediate states in which no external strains are applied and the aggregate relaxes towards the equilibrium. That is, equilibrium is guaranteed by fluctuation fields that allow particles deformation to deviate from the applied average strain.

The response of the aggregate is represented by stress-strain and volume strain-strain curves in which the deviatoric stress is normalized by the pressure $p_0$ while the strains by the volume strain $\Theta_0$.

In Fig. 1 we show the response of the aggregate: the evolution of the normalized deviatoric stress and its partition in $q^N$ and $q^T$, top; the evolution of the normalized volume strain in which dilatancy is present, bottom. As in Jenkins and Strack [23] and Thorthon and Antony [27], the contribution to the deviatoric stress associated with $q^N$ and $q^T$ are illustrative. We can distinguish an initial almost elastic response of the aggregate, for $0 \le \gamma/\Theta_0 \lesssim 0.2$, followed by an elasto-plastic response in which $q^T$ does not vary, essentially. In this latter regime, a typical hardening behavior occurs in which the incremental deviatoric stress depends only on the normal contact force while particles slide/roll and irreversibility occurs. Here, the incremental tangential component is assumed to be approximately zero, so $\dot{q}_T = 0$.

We also consider the network of contact vectors and introduce the fabric tensor (e.g. [27, 28]):

$$ \mathbf{\Phi}_{ij} = \frac{1}{N_c} \sum_{c=1}^{N_c} \hat{d}_i^c \hat{d}_j^c \quad (3) $$

We measure the evolution of the deviatoric fabric, $\hat{\mathbf{\Phi}} = \mathbf{\Phi}_{33} - (\mathbf{\Phi}_{11} + \mathbf{\Phi}_{22})/2$, normalized by the initial fabric
<figure><img src="image_3.png" /><figcaption>Fig. 1 Uniaxial deformation at constant pressure, numerical simulations: normalized deviatoric stress (top), normalized volume strain (bottom)</figcaption></figure>

$\Phi_{iso} = \Phi_{11} = \Phi_{22} = \Phi_{33}$, during the loading. In Fig. 2, top, we show the evolution of the fabric, $\hat{\Phi}/\Phi_{iso}$.

The fabric tensor, $\Phi_{ij}$, is also identified as the second moment of the contact distribution $f(\hat{\mathbf{d}})$ whose approximated expression is (e.g. [29]):

$$f(\hat{\mathbf{d}}) = \frac{k}{4\pi} \left[ 1 + \frac{15}{2}(\hat{d}_i\hat{d}_j - \frac{1}{3}\delta_{ij})(\Phi_{ij} - \frac{1}{3}\delta_{ij}) \right] \quad (4)$$

with $\int_{\Omega} f(\hat{\mathbf{d}})d\Omega = k$ and $\hat{\mathbf{d}}(\sin\theta \cos\varphi, \sin\theta \sin\varphi, \cos\theta)$ in which $\theta$ is the polar angle with respect to the axis of compression $y_3$ and $k$ is the coordination number, the average number of contacts per particle. If we knew the contact

<figure><img src="image_4.png" /><figcaption>Fig. 2 Uniaxial deformation at constant pressure, numerical simulations: fabric evolution (top), contact distribution (bottom)</figcaption></figure>

distribution we could have determined the second moment through the following equation (e.g. [30]):

$$\Phi_{ij} = \int_{\Omega} \hat{f}(\hat{\mathbf{d}}) \hat{d}_i \hat{d}_j d\Omega \quad (5)$$

in which $d\Omega = \sin\theta d\theta d\varphi$. In Fig. 2, bottom, we show the normalized contact distribution, $\hat{f}(\hat{\mathbf{d}}) = f(\hat{\mathbf{d}})/k$, at different strain states (for symmetry reason we look at $0 \le \theta \le \pi/2$). Both Fig. 2, top and bottom, show a development of the structural anisotropy during the loading with a region of increasing contacts around the pole of compression, about $\theta = 0$, and decreasing about $\theta = \pi/2$ associated with deletion.
<figure><img src="image_5.png" /><figcaption>**Fig. 3** Uniaxial deformation at constant pressure, numerical simulation: evolution of the coordination number with the normalized shear strain (solid line); polynomial approximation (dashed line)</figcaption></figure>

In Fig. 3 we show the evolution of the coordination number $k$, the average number of contacts per particle, whose initial value is $k_0 \approx 6$. We also evaluate a best fit approximation for the evolution of the coordination number, with a polynomial function, which offers the following relation: $k = ay^2/\Theta_0 + by/\Theta_0 + c$ with $a = 0.5, b = -1.6$ and $c = 6$.

At this point, incremental probes are applied at given stressed states. Because of the uniaxial compression, these stressed states are anisotropic and the incremental response depends on the magnitude of the perturbation [31–33]: we have an elastic response that corresponds to an unloading condition with respect to the previous triaxial loading, and an elasto-plastic response which corresponds to a forward loading. We are interested to evaluate $W_2$ when an elasto-plastic increment is applied because we expect to characterize the macroscopic stiffness tensor without the major symmetry [17]. In the last section, we determine $W_2$ under these particular conditions and test the theory against the simulation data.

## 2.2 Theory

The model is based on the theory of fluctuations proposed by Jenkins et al. [34] and developed in case of localization by La Ragione et al. [17]. We assume that contacting particles deviate from the applied average strain because of equilibrium.

We write the kinematic that governs the problem, we assume the contact law for contacting particles and solve, approximately, equilibrium. With a proper statistical hypothesis, we move from the local pair interaction to the aggregate in order to determine the average stress from which it is possible to obtain an expression for the second order work. It is the nature

<figure><img src="image_6.png" /><figcaption>**Fig. 4** The pair and its neighboord: particles A and B fluctuate while particle n is constrained to move with the average deformation</figcaption></figure>
of the macroscopic stiffness that relates incremental in strain with incremental in stress that provides an explanation of negative second order work.

## 2.2.1 Kinematics

We phrase the problem incrementally and we restrict our attention to the hardening regime in which, as seen in Fig. 1 (top), the macroscopic response can be approximated as incrementally frictionless. That is, we refer to incremental forward loading. In this regime particles slide, so plasticity occurs but we only have increments in normal components of the contact forces, approximately. Consequently, we focus only on the kinematic of contacting particles related to the motion of their centers. The incremental relative displacement between contacting particles A and B is given by:

$$ \dot{u}_i^{(BA)} = \dot{c}_i^{(B)} - \dot{c}_i^{(A)} \qquad (6) $$

in which $\dot{c}_i$ is the incremental particle center displacement. We allow particles to deviate from the average strain and introduce a fluctuation field. So we write

$$ \dot{c}_i^{(A)} = \bar{\dot{c}}_i + \tilde{c}_i^{(A)} \qquad (7) $$

in which we make explicit the average contribution, $\bar{\dot{c}}_i$, and the fluctuation, $\tilde{c}_i^{(A)}$. A simple way to enrich the kinematics is to focus on a pair A-B and its neighborhoods (see Fig. 4) and assume that the relative displacement of the contact point of a typical pair A-B is given by

$$ \dot{u}_i^{(BA)} = \dot{E}_{ij} d_j^{(BA)} + \dot{\Delta}_i^{(BA)} \qquad (8) $$

while for particle n in contact with A is

$$ \dot{u}_i^{(nA)} = \dot{c}_i^{(n)} - \dot{c}_i^{(A)} \qquad (9) $$

in which only particle A fluctuates. So

$$ \dot{u}_i^{(nA)} = \dot{E}_{ij} d_j^{(nA)} + \frac{1}{2} \dot{\Delta}_i^{(BA)} - \frac{1}{2} \dot{\Sigma}_i^{(BA)}. \qquad (10) $$

The difference in fluctuation is

$$ \dot{\Delta}_i^{(BA)} = \tilde{c}_i^{(B)} - \tilde{c}_i^{(A)} \qquad (11) $$

while the sum is

$$ \dot{\Sigma}_i^{(BA)} = \tilde{c}_i^{(A)} + \tilde{c}_i^{(B)}. \qquad (12) $$

Similar expressions can be employed for particle B and its neighborhood. The contact vector is indicated by $d_j^{(BA)} = D\hat{d}_j^{(BA)}$.

## 2.2.2 Statics

We write the increment $\dot{\mathbf{F}}^{(BA)}$ in the contact force exerted by particle B on particle A in terms of the increment $\dot{\mathbf{u}}^{(BA)}$ in the relative displacement of the points of contact:

$$ \dot{\mathbf{F}}_i^{(BA)} = K_{ij}^{(BA)} \dot{u}_j^{(BA)} \qquad (13) $$

in which, in case of incremental frictionless behavior, the contact stiffness is

$$ K_{ij}^{(BA)} = K_N \hat{d}_i^{(BA)} \hat{d}_j^{(BA)}. \qquad (14) $$

The normal contact stiffness depends on the normal component $\delta$ of the compressive displacement of the centers of the particles, the diameter D of the spheres, and their material properties:

$$ K_N^{(BA)} = \frac{GD^{1/2}}{(1 - \nu)} \delta^{1/2}. \qquad (15) $$

While the incremental displacement of the centers of contacting particles is given in terms of the averages of the increments, the total normal component, $\delta = u_i^{(BA)} \hat{d}_i^{(BA)}$ in Eq. (15), is given in terms of the average strain by

$$ \delta = -\hat{d}_i^{(BA)} E_{ij} d_j^{(BA)}. \qquad (16) $$

This is how the stiffness depend upon the existing average strain and it represents a measure of the strain level of the aggregate along the triaxial loading. Equation (16) can be written as:

$$ \delta = \frac{D}{3} (\Theta - 2\gamma + 6\gamma \cos^2 \theta). \qquad (17) $$

## 2.2.3 Equilibrium

We determine the fluctuations that ensure the equilibrium of the pair A-B through the following equations. For particle A equilibrium is

$$ \sum_{n \neq B}^{N^{(A)}} \dot{\mathbf{F}}_i^{(nA)} + \dot{\mathbf{F}}_i^{(BA)} = 0 \qquad (18) $$

while for particle B

$$ \sum_{n \neq A}^{N^{(B)}} \dot{\mathbf{F}}_{i}^{(nB)} + \dot{\mathbf{F}}_{i}^{(AB)} = 0, \qquad (19) $$

in which, for example, $N^A$ is the number of particles in contact with particle A. Difference of the equilibrium equations leads to
$$ \Delta_i^{(BA)} = -D \left[ A_{ik}^{(BA)-1} J_{klm}^{(BA)} - A_{ik}^{(AB)-1} J_{klm}^{(AB)} \right] \dot{E}_{lm} \quad (20) $$

in which, for example,

$$ A_{ik}^{(BA)} = \sum_{n=1}^{N^{(A)}} K_{N}^{(nA)} \hat{d}_{i}^{(nA)} \hat{d}_{k}^{(nA)} \quad (21) $$

and

$$ J_{klm}^{(BA)} = \sum_{n=1}^{N^{(A)}} K_{N}^{(nA)} \hat{d}_{k}^{(nA)} \hat{d}_{l}^{(nA)} \hat{d}_{m}^{(nA)}. \quad (22) $$

We employ the conditional average defined over all pairs whose contact vector is within an increment of solid angle:

$$ \begin{aligned} \langle A_{ik} \rangle_{\hat{d}^{(BA)}} &= \left\langle \sum_{n=1}^{N^{(A)}} K_N \hat{d}_i^{(nA)} \hat{d}_k^{(nA)} \right\rangle_{\hat{d}^{(BA)}} \\ &= \left\langle \sum_{n=1}^{N^{(B)}} K_N \hat{d}_i^{(nB)} \hat{d}_k^{(nB)} \right\rangle_{\hat{d}^{(AB)}} \end{aligned} \quad (23) $$

and

$$ \begin{aligned} \langle J_{klm} \rangle_{\hat{d}^{(BA)}} &= \left\langle \sum_{n=1}^{N^{(A)}} K_N^{(nA)} \hat{d}_k^{(nA)} \hat{d}_l^{(nA)} \hat{d}_m^{(nA)} \right\rangle_{\hat{d}^{(BA)}} \\ &= - \left\langle \sum_{n=1}^{N^{(B)}} K_N^{(nB)} \hat{d}_k^{(nB)} \hat{d}_l^{(nB)} \hat{d}_m^{(nB)} \right\rangle_{\hat{d}^{(AB)}} \end{aligned} \quad (24) $$

because $\hat{d}^{(BA)} = -\hat{d}^{(AB)}$. Eq. (20) becomes

$$ \langle \dot{\Delta}_i \rangle_{\hat{d}} = -2D \langle A_{ik}^{-1} J_{klm} \rangle_{\hat{d}} \dot{E}_{lm} \quad (25) $$

with the following representation for the third rank tensor

$$ \langle A_{ik}^{-1} J_{klm} \rangle_{\hat{d}} = \beta_1 \hat{d}_i^{(BA)} \hat{d}_l^{(BA)} \hat{d}_m^{(BA)} + \beta_2 (\delta_{il} \hat{d}_m^{(BA)} + \delta_{lm} \hat{d}_i^{(BA)} + \delta_{im} \hat{d}_l^{(BA)}). \quad (26) $$

The coefficients $\beta_1$ and $\beta_2$ depend on the coordination number $k$. The summation in the conditional averages employed in Eqs. (23) and (24) are replaced by an integration over all the directions in the neighbors of all pairs characterized by a contact vector $\hat{d}$ within an incremental solid angle, assuming a uniform distribution. With a simple statistical model, we obtain

$$ \beta_1 = \frac{3}{k} \left( \frac{166 - 11k}{128} \right) \quad (27) $$

and

$$ \beta_2 = -\frac{3}{k} \left( \frac{k + 14}{128} \right). \quad (28) $$

We are neglecting the contribution of anisotropy through the contact stiffness $K_N$ so fluctuations depend only on the contacts network. Both $\beta_1$ and $\beta_2$ decreases as $k$ increases, or, equivalently, the fluctuations contribution diminishes as the confining pressure increases [35].

Given Eq. (8), we derive:

$$ \begin{aligned} \langle \dot{u}_i \rangle_{\hat{d}^{(BA)}} = & \dot{E}_{ij} \hat{d}_j^{(BA)} - 2D \beta_1 \hat{d}_i^{(BA)} \hat{d}_l^{(BA)} \hat{d}_m^{(BA)} \dot{E}_{ml} \\ & - 2D \beta_2 (\delta_{il} \hat{d}_m^{(BA)} + \delta_{lm} \hat{d}_i^{(BA)} + \delta_{im} \hat{d}_l^{(BA)}) \dot{E}_{ml}. \end{aligned} \quad (29) $$

Therefore, for a given pair orientated along $\hat{d}$, we write the incremental contact force as

$$ \langle \dot{F}_i \rangle_{\hat{d}} = K_N \hat{d}_i \hat{d}_j \langle \dot{u}_j \rangle_{\hat{d}} \quad (30) $$

or

$$ \langle \dot{F}_i \rangle_{\hat{d}} = DK_N \hat{d}_i [\hat{d}_q \dot{E}_{qj} \hat{d}_j (1 - 2\beta_1 - 4\beta_2) - 2\beta_2 \delta_{lm} \dot{E}_{ml}]. \quad (31) $$

With respect to the average strain theory, it is the last term in the bracket, proportional to the volume strain and $\beta_2$, that characterizes the new ingredient in the model from the qualitative point of view, once fluctuations are included.

### 2.2.4 Incremental stress

An analytical expression for the stress increment may be obtained by employing the continuous analog of Eq. (2). This is phrased in terms of the number of particles per unit volume $n$ and a contact distribution function $f(\hat{d})$. In the pre-peak region, range of deformation of our interest, the distribution of contacts is assumed to be isotropic, $f(\hat{d}) = k/4\pi$, where $k$ is the coordination number that evolves during the loading. Therefore, we account for the induced anisotropy through Eqs. (15) and (17) while we ignore the structural anisotropy which seems reasonable according to Fig. 2 in which, at most, we have, in the region of compression around $\theta = 0$, a slight increments of number of contacts with respect to the isotropic distribution.

The average stress is:

$$ \dot{\sigma}_{ij} = \frac{n}{2} \int_{\Omega} f(\hat{d}) \langle \dot{F}_i \rangle_{\hat{d}} d_j d\Omega \quad (32) $$

or, given Eqs. (15), (17) and (31),

$$ \begin{aligned} \dot{\sigma}_{ij} = & \omega(1 - 2\beta_1 - 4\beta_2) \\ & \times \int_0^{2\pi} \int_0^{\theta_1} R^{1/2} \hat{d}_k \hat{d}_m \hat{d}_i \hat{d}_j \dot{E}_{km} \sin \theta d\theta d\varphi \\ & - 2\omega\beta_2 \int_0^{2\pi} \int_0^{\theta_1} R^{1/2} \hat{d}_i \hat{d}_j \delta_{km} \dot{E}_{km} \sin \theta d\theta d\varphi \end{aligned} \quad (33) $$
in which $R = \Theta - 2\gamma + 6\gamma \cos^2 \theta$, $\omega \equiv nkGD^3 / [4\pi\sqrt{3}(1-\nu)]$ and $\theta_1$ is the angle of deletion, when $\delta = 0$ in Eq. (17):

$$ \theta_1 = \pi/2 \qquad (34) $$

if $\Theta - 2\gamma \ge 0$, otherwise

$$ \theta_1 = \arccos \sqrt{(2\gamma - \Theta)/(6\gamma)}. \qquad (35) $$

The upper limit of the integration, $\theta_1$, account for deletion when the average contact displacement, $\delta$, goes to zero. In compact form,

$$ \dot{\sigma}_{ij} = A_{ijkm} \dot{E}_{km} \qquad (36) $$

where the stiffness tensor, $A_{ijkm}$, can be represented through the axial vector $\mathbf{h} = \mathbf{y}_3$ and Kronecker delta, $\delta_{ij}$, with

$$ \begin{aligned} A_{ijkm} = & \eta_1 h_m h_k h_j h_i + \eta_2 \delta_{ij} \delta_{km} + \eta_3 (\delta_{ik} \delta_{jm} + \delta_{im} \delta_{jk}) \\ & + \eta_4 \delta_{km} h_j h_i + \eta_5 h_m h_k \delta_{ij} \\ & + \eta_6 (\delta_{ik} h_j h_m + \delta_{jm} h_k h_i + \delta_{im} h_k h_j + \delta_{jk} h_i h_m). \end{aligned} \qquad (37) $$

The six coefficients are

$$ \eta_1 = 2\pi\omega(1 - 2\beta_1 - 4\beta_2)a_1, $$

$$ \eta_2 = 2\pi\omega(1 - 2\beta_1 - 4\beta_2)a_2 - 4\pi\omega\beta_2 a_4, $$

$$ \eta_3 = 2\pi\omega(1 - 2\beta_1 - 4\beta_2)a_2, $$

$$ \eta_4 = 2\pi\omega(1 - 2\beta_1 - 4\beta_2)a_3 - 4\pi\omega\beta_2(a_5 - a_4), $$

$$ \eta_5 = \eta_6 = 2\pi\omega(1 - 2\beta_1 - 4\beta_2)a_3, $$

with

$$ a_1 \equiv \int_0^{\theta_1} R^{1/2} \left( \cos^4 \theta - 3 \sin^2 \theta \cos^2 \theta + \frac{3}{8} \sin^4 \theta \right) \sin \theta d\theta, $$

$$ a_2 \equiv \int_0^{\theta_1} R^{1/2} \frac{1}{8} \sin^5 \theta d\theta, $$

$$ a_3 \equiv \int_0^{\theta_1} R^{1/2} \left( \frac{1}{2} \sin^2 \theta \cos^2 \theta - \frac{1}{8} \sin^4 \theta \right) \sin \theta d\theta, $$

$$ a_4 \equiv \int_0^{\theta_1} R^{1/2} \frac{1}{2} \sin^3 \theta d\theta, $$

and

$$ a_5 \equiv \int_0^{\theta_1} R^{1/2} \cos^2 \theta \sin \theta d\theta. \qquad (38) $$

As suggested by La Ragione et al. [17, 36], it is the presence of fluctuations in the kinematics of contacting particles that leads to $\eta_4 \ne \eta_5$ which implies the loss of major symmetry for the anisotropic stiffness tensor $A_{ijkm}$. It is the term proportional to $\beta_2$ and the volume strain in Eq. (31) that makes clear the difference between the simple average strain and

the presence of fluctuations. Without fluctuations, $A_{ijkm}$ is symmetric. The asymmetry in $A_{ijkm}$ is a necessary condition for the second order work to vanish at some stressed state, under a particular perturbation.

### 3 Second-order work

The second order work (e.g. [14]) is

$$ W_2 = \dot{\sigma}_{ij} \dot{E}_{ij}. \qquad (39) $$

or, with Eq. (36)

$$ W_2 = A_{ijkm} \dot{E}_{km} \dot{E}_{ij}. \qquad (40) $$

Given Eq. (37), we obtain

$$ W_2 = \eta_1 h_m h_k h_j h_i \dot{E}_{km} \dot{E}_{ij} + \eta_2 \dot{E}_{qq} \dot{E}_{qq} + 2\eta_3 \dot{E}_{ij} \dot{E}_{ij} \\ + (\eta_4 + \eta_5) \dot{E}_{qq} h_i h_j \dot{E}_{ij} + 4\eta_6 h_j h_m \dot{E}_{sm} \dot{E}_{sj}. \qquad (41) $$

The coefficients $\eta$'s depend on the relation between the volume strain $\Theta$ and the shear strain $\gamma$ which are present in the contact stiffness $K_N$. Here, following Jenkins and Strack [23], we use the simple average strain and determine the initial pressure

$$ p = \frac{nGD^2k_0}{27\sqrt{3}(1-\nu)} \Theta_0^{3/2} \qquad (42) $$

and then equate it to the pressure when an uniaxial deformation is applied:

$$ \int_0^{\theta_1} R^{3/2} \sin \theta d\theta - \frac{k_0}{k} \Theta_0^{3/2} = 0. \qquad (43) $$

Equation (43) provides a relation between $\gamma$ and $\Theta$ and with the evolution of the coordination number, $k = 0.5\gamma^2/\Theta_0 - 1.6\gamma/\Theta_0 + 6$, previously derived (see Fig. 3), we can predict $W_2$. We focus on the increments that may provide negative values of $W_2$, according with Hadda et al. [37]. We consider $\dot{E}_{33} < 0$ and $\dot{E}_{11} = \dot{E}_{22} > 0$, so

$$ W_2 = \eta_1 \dot{E}_{33}^2 + \eta_2 (2\dot{E}_{11} + \dot{E}_{33})^2 + 2\eta_3 (2\dot{E}_{11}^2 + \dot{E}_{33}^2) \\ + (\eta_4 + \eta_5) (2\dot{E}_{11} + \dot{E}_{33}) \dot{E}_{33} + 4\eta_6 \dot{E}_{33}^2. \qquad (44) $$

In particular, in the Rendulic strain plane, see Fig. 5 (strain positive in compression), we take $\dot{E}_{33} = \rho \sin(\alpha_e)$ and $\dot{E}_{11} = \dot{E}_{22} = -\rho \cos(\alpha_e)/\sqrt{2}$, in which $\rho$ is a positive constant and $\alpha_e$ indicates the incremental strain direction in the range $110 \le \alpha_e \le 160$. These are probes that we identify as forward incremental loadings for which the theory has been developed.
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Circular diagrams of the normalized second order work. Focus on the second quadrant. First negative $W_2$ at $\alpha_{\epsilon} = 130^{\circ}$ when $\gamma/\Theta_0 = 0.55$ according to numerical simulations</figcaption></figure>

<figure><img src="image_8.png" /></figure>

We predict a negative value for $W_2$ when $\gamma/\Theta_0 = 0.93$, with $\alpha = 150^{\circ}$ (see Matlab code in Appendix). This angle is close to what also seen in numerical simulations [37].

We also evaluate $W_2$ through our numerical simulations. At each stressed states, along the triaxial test, we apply probes to calculate the second order work through Eq. (39). The value of $\rho$ is crucial as we consider probes associated with an elasto-plastic response [38]. So, we

apply $\rho \approx 10^{-6}$ as in the triaxial test. In Fig. 6 we plot the normalized second order work, $(\dot{\sigma}_{ij}\dot{E}_{ij})/(|\dot{\sigma}|\rho)$. A first negative value for $W_2$ occurs at $\gamma/\Theta_0 = 0.55$ and $\alpha_{\epsilon} = 130^{\circ}$. Both the strain and the angle are smaller than those predicted by the theory. From simulations, we note an indication of the preferential angle of failure since the first stressed states.
<figure><img src="image_9.png" /><figcaption>**Fig. 7** Evolution of the term responsible of the asymmetry of the stiffnnes tensor $A_{ijkm}$ with the normalized shear strain, see the coefficient $\eta_4$ in Eq. (37)</figcaption></figure>

The evolution of the asymmetry of the stiffnnes tensor $A_{ijkm}$, in the hardening regime, is visualized in Fig. 7.

A possible way to improve the theory is to employ fluctuations in the stiffness, Eq. (16), that will lead to a more complicated relation between $\gamma$ and $\Theta$, now given by Eq. (43), and other terms in the coefficients of the macroscopic stiffness tensor $A_{ijkm}$, without changing the structure of the tensor. It is also possible to include fabric and structural anisotropy in both fluctuations, through Eq. (25) and stress, through Eq. (32). However, we think that the present model is more transparent, and it makes clear and illustrative the role played by the fluctuations in relation to the vanishing of the second order work.

Finally, we take Eq. (44) and dived it by $\dot{E}_{11}^2$ and, with $\chi = \dot{E}_{33}/\dot{E}_{11}$, we obtain

$$ \frac{W_2}{\dot{E}_{11}^2} = \eta_1 \chi^2 + \eta_2 (2 + \chi^2)^2 + 2\eta_3 (2 + \chi^2)^2 + (\eta_4 + \eta_5) (2\chi + \chi^2)^2 + 4\eta_6 \chi^2 \quad (45) $$

or

$$ \frac{W_2}{\dot{E}_{11}^2} = A_1 \chi^2 + A_2 \chi + A_3 \quad (46) $$

with

$$ A_1 = \eta_1 + \eta_2 + 2\eta_3 + \eta_4 + \eta_5 + 4\eta_6, \quad (47) $$

$$ A_2 = 4\eta_2 + 2\eta_4 + 2\eta_5, \quad (48) $$

$$ A_3 = 4(\eta_2 + \eta_3). \quad (49) $$

The discriminant is $A_2^2 - 4A_1A_3$ and it is zero at $\gamma/\Theta_0 = 0.93$.

## 4 Conclusion

We have employed a model that characterizes the second work with respect to failure. We have introduced a kinematics of contacting particles that deviates from the average strain and written an explicit but approximate equations of equilibrium for pairs of contacting particles in the aggregate. Because of fluctuations, it is possible, in the hardening regime, to predict a macroscopic stiffness tensor without the major symmetry which is essential to have a stressed, anisotropic state in which $W_2$ goes to zero. That is, in the context of a simplified model, we have linked the extend kinematics of contacting particles with fluctuations to the asymmetry of the macroscopic stiffness tensor, and to the possibility that the second work $W_2$ vanishes. The present model overpredicts the angle of failure and the value of the strain when the second order work vanishes. On the other hand, the theory exhibits the important ingredients from the microscopic point of view that may explain failure of the aggregate.

## Appendix

clear; clc;

% k is the coordination number

% c is the normalized shear strain

% a is the normalized volume strain

% strain negative in compression

E3 = -sin(30*pi/36); %it means compression along the vertical direction

E1 = -cos(30*pi/36)/sqrt(2); % it means traction along the horizontal direction

% input for the volume strain and shear strain from Eq. 33

load('input_average_10k.mat')

syms x

for j=1:100

% coordination number evolution

kn(j)=0.5+gp(j).+gp(j)-1.6+gp(j)+6;

% Coefficient for the fluctuation Delta, beta1 and beta2, evaluated from Jenkins et al. (JMPS

2005)

beta1(j) = 3/kn(j)=(166-11+kn(j))/128;

beta2(j) = -3/kn(j)=(kn(j)+14)/128;

%% Lp and gp are the volume strain and shear strain given by eq. 33

a(j)=Lp(j);

c(j)=gp(j);

thn(j)=thp(j);

% Numerical integration of the integrals for a1,a2,a3,a4,a5

TR1 =@(x) sqrt(a1(j)-2+c(j))+6+c(j).*.(cos(x)).^2).*.(cos(x)).^4-3+sin(x).^2.*.cos(x).^3+8+sin(x).^4).*.(sin(x));

a1(j) = integral(TR1,0,thn(j));

TR2 =@(x) sqrt(a1(j)-2+c(j))+6+c(j).*.(cos(x)).^2).*1/8.*sin(x).^5;

a2(j) = integral(TR2,0,thn(j));

TR3 =@(x) sqrt(a1(j)-2+c(j))+6+c(j).*.(cos(x)).^2).*1/(2+sin(x).^2).*cos(x).^2-1/8+sin(x).^4).*.(sin(x));

a3(j) = integral(TR3,0,thn(j));

TR4 =@(x) sqrt(a1(j)-2+c(j))+6+c(j).*.(cos(x)).^2).*1/(2+sin(x).^3);

a4(j) = integral(TR4,0,thn(j));

TR5 =@(x) sqrt(a1(j)-2+c(j))+6+c(j).*.(cos(x)).^2).*.(cos(x).^2).*.(sin(x));

a5(j) = integral(TR5,0,thn(j));

eta1(j) = real(a1(j))+(1-2+beta1(j))-4+beta2(j));

eta2(j) = (1-2+beta1(j)-4+beta2(j))*real(a2(j))-2+beta2(j)+real(a4(j));

eta3(j) = (1-2+beta1(j)-4+beta2(j))*real(a2(j));

eta4(j) = (1-2+beta1(j)-4+beta2(j))*real(a3(j))-2+beta2(j)+real(a5(j))-real(a4(j));

eta5(j) = (1-2+beta1(j)-4+beta2(j))*real(a3(j));

W2(j)=eta1(j)+E3^2+eta2(j)+(2+E1+E3)^2+2+eta3(j)*(2+E1^2+E3^2)+
(eta4(j)+eta5(j))+(2+E1+E3)+E3+4+eta5(j)+E3^2;

end

**Acknowledgements** Luigi La Ragione is grateful to the Gruppo Nazionale della Fisica Matematica (G.N.F.M.), Italy and CN-00000013 - National Centre for HPC. Big Data and Quantum Computing - Spoke 5 "Environment and Natural Disasters".
**Funding** Open access funding provided by Politecnico di Bari within the CRUI-CARE Agreement.

**Declarations**

**Compliance with Ethical Standards** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

**Competing interest** The authors have no competing interests to declare that are relevant to the content of this article.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

1. Nicot, F., Darve, F.: Basic features of plastic strains: from micro-mechanics to incrementally nonlinear models. Int. J. Plast. **23**, 1555–1588 (2007)

2. Nicot, F., Darve, F.: Micro-mechanical bases of some salient constitutive features of granular materials. Int. J. Solids Struct. **44**, 7420–7443 (2007)

3. Tordesillas, A.: Force chain buckling, unjamming transitions and shear banding in dense granular assemblies. Phil. Mag. **87**(32), 4987–5016 (2007)

4. Tordesillas, A., Muthuswamy, M.: On the modeling of confined buckling of force chains. J. Mech. Phys. Solids **57**(4), 706–727 (2009)

5. Walker, D.M., Tordesillas, A.: Topological evolution in dense granular materials: a complex networks perspective. Int. J. Solids Struct. **47**, 624–639 (2010)

6. Radjai, F., Roux, J.N., Daouadji, A.: Modeling granular materials: century-long research across scales. J. Eng. Mech. **143**(4), 04017002–20 (2017)

7. Glansdorff, P., Prigogine, I.: Thermodynamic Theory of Structure, Stability and Fluctuations. Wiley-Interscience, New York (1971)

8. Nicolis, G., Prigogine, I.: Self-Organization in Nonequilibrium Systems: From Dissipative Structures to Order through Fluctuations. John Wiley & Sons, New York (1977)

9. Rice, J.R.: Continuum mechanics and thermodynamics of plasticity in relation to microscale deformation mechanisms. In: Argon, A.S. (ed.) Constitutive Equations in Plasticity, pp. 23–79. MIT Press, Cambridge (1975)

10. Rudnicki, J.W., Rice, J.R.: Conditions for the localization of deformations in pressure-sensitive dilatant materials. J. Mech. Phys. Solids **23**, 371–394 (1975)

11. Hill, R.: A general theory of uniqueness and stability in elastic-plastic solids. J. Mech. Phys. Solids **6**(3), 236–249 (1958)

12. Nicot, F., Sibille, L., Darve, F.: Failure in rate-independent granular materials as a bifurcation toward a dynamic regime. Int. J. Plast. **29**, 136–154 (2012)

13. Wan, R., Nicot, F., Darve, F.: Failure in Geomaterials, a Contemporary Treatise. Wiley (2017)

14. Nicot, F., Darve, F.: Diffuse and localized failure modes: two competing mechanisms. Int. J. Num. Anal. Method. Geomech. **35**, 586–601 (2011)

15. Vardoulakis, I.: Equilibrium theory of shear bands in plastic bodies. Mech. Res. Commun. **3**, 209–214 (1976)

16. Vardoulakis, I.: Shear band inclination and shear modulus of sand in biaxial tests. Int. J. Numer. Anal. Methods Geomech. **4**, 103–119 (1980)

17. La Ragione, L., Prantil, V.C., Jenkins, J.T.: A micromechanical prediction of localization in a granular material. J. Mech. Phys. Solids **83**, 146–159 (2015)

18. Wren, J.R., Borja, R.I.: Micromechanics of granular media Part II: overall tangential moduli and localization model for periodic assemblies of circular disks. Comput. Methods Appl. Mech. Eng. **141**, 221–246 (1997)

19. Darve, F., Flavigny, E., Meghachou, M.: Constitutive modeling and instabilities of soil behaviour. Comp. Geotech. **17**, 203–224 (1995)

20. Nicot, F., Darve, F., Koah, H.-D.-V.: Bifurcation and second-order work in geomaterials. Int. J. Num. Anal. Meth. Geomech. **31**, 1007–1032 (2007)

21. Nicot, F., Sibille, L., Donze, F., Darve, F.: From microscopic to macroscopic second-order work in granular assemblies. Mech. Mat. **39**, 664–684 (2007)

22. Cundall, P.A., Strack, O.D.L.: A discrete numerical model for granular assemblies. GEOT **29**, 47–65 (1979)

23. Jenkins, J.T., Strack, O.D.L.: Mean-field inelastic behavior of random arrays of identical spheres. Mech. Mat. **16**, 25–33 (1993)

24. Nicot, F., Darve, F.: A multiscale approach to granular materials. Mech. Mat. **37**(9), 980–1006 (2005)

25. La Ragione, L., Jenkins, J.T.: The initial response of an idealized granular material. Proc. R. Soc. A. **463**, 735–758 (2007)

26. Magnanimo, V., La Ragione, L., Jenkins, J.T., Wang, P., Makse, H.A.: Characterizing the shear and bulk moduli of an idealized granular material. Europhys. Lett. **81**, 34006–27 (2008)

27. Thornton, C., Antony, S.J.: Quasi-static deformation of particulate media. Philos. Trans. R. Soc. Lond. Ser. A **356**, 2763–2782 (1998)

28. Luding, S.: Micro-macro transition for anisotropic, frictional granular packings. Int. J. Solids Struct. **41**, 5821–5836 (2004)

29. Adani, S.G., Tucker, C.L.: The use of tensors to describe and predict fiber orientation in short fiber composites. J. Rheol. **31**(8), 751–784 (1987)

30. Jenkins, J.T., Seto, R., La Ragione, L.: Predictions of microstructure and stress in planar extensional flows of a dense viscous suspension. J. Fluid Mech. **927**(A27), 1–20 (2021)

31. Kuhn, M., Daouadji, A.: Quasi-static incremental behavior of granular materials: elastic-plastic coupling and micro-scale dissipation. J. Mech. Phys. Solids **114**, 219–237 (2018)

32. Recchia, G., Magnanimo, V., Cheng, H., La Ragione, L.: DEM simulation of anisotropic granular materials: elastic and inelastic behavior. Gran. Matter **22**(85), 1–13 (2020)

33. Cheng, H., Luding, S., Saitoh, K., Magnanimo, V.: Elastic wave propagation in dry granular media: effects of probing characteristics and stress history. Int. J. Solids Struct. **187**, 85–99 (2020)

34. Jenkins, J.T., Johnson, D.L., La Ragione, L., Makse, H.: Fluctuations and the effective moduli of an isotropic, random aggregate of identical, frictionless spheres. J. Mech. Phys. Solids **53**, 197–225 (2005)

35. Makse, H.A., Gland, N., Johnson, D.L., Schwartz, L.: Granular packings: nonlinear elasticity, sound propagation, and collective relaxation dynamics. Phys. Rev. E **70**, 061302 (2004)

36. La Ragione, L., Oger, L., Recchia G., Sollazzo, A.: Anisotropy and lack of symmetry for a frictionless, random, aggregate of elastic particles: theory and numerical simulations. Proc. R. Soc. A, **471** (2015)
37. Hadda, N., Nicot, F., Bourrier, F., Sibille, L., Radjai, F., Darve, F.: Micromechanical analysis of second order work in granular media. Gran. Matter **15**(2), 221–235 (2013). (**Springer Verlag**)

38. La Ragione, L.: The incremental response of a stressed, anisotropic granular material: loading and unloading. J. Mech. Phys. Solids **95**, 147–168 (2016)

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_10.png" /></figure>