# Shear banding as a dissipative structure from a thermodynamic viewpoint

F. Nicot, X. Wang, A. Wautier, R. Wan, F. Darve

► To cite this version:

F. Nicot, X. Wang, A. Wautier, R. Wan, F. Darve. Shear banding as a dissipative structure from a thermodynamic viewpoint. Journal of the Mechanics and Physics of Solids, 2023, 179, pp.105394. 10.1016/j.jmps.2023.105394. hal-04184141

HAL Id: hal-04184141

https://hal.science/hal-04184141v1

Submitted on 21 Aug 2023

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 CC BY 4.0 - Attribution - International License
# Shear banding as a dissipative structure from a thermodynamic viewpoint

F. Nicot¹, X. Wang¹,², A. Wautier³, R. Wan⁴, and F. Darve⁵

(1) Université Savoie Mont Blanc, ISTerre, Chambéry (France)

(2) University of Science and Technology Beijing, Dept. of Civil Engineering, Beijing (China)

(3) Aix-Marseille Université, INRAE, Unité de Recherche RECOVER, Aix-en-Provence (France)

(4) University of Calgary, Calgary (Canada)

(5) Université Grenoble-Alpes, Laboratoire Sols Solides Structures, Grenoble (France)

Corresponding author: francois.nicot@univ-smb.fr

## ABSTRACT:

Granular materials are now known to be an illustration of complex materials as they display emergent macroscopic properties when loaded. An initially homogenous response can bifurcate into a heterogeneous one with the appearance of a rich variety of structured kinematical patterns. The shear banding that ensues illustrates a symmetry-breaking transition with multiple choices of macroscopic behaviours, a common feature of dynamical complex systems. Even though the phenomenon has been studied for decades, this regime transition remains mostly mysterious in geomaterials, with no convincing arguments that could link it to the underlying microscopic mechanisms. The paper investigates this issue by invoking the fundamental minimum entropy production theorem established by Prigogine in the past century to seek any connection with the second-order work theory in the mechanics of failure. A general equation linking the derivatives of the entropy of a mechanical system to the second-order work is thus inferred, which leads to a thermodynamic interpretation of bifurcations in the failure behaviour of granular materials under a given loading. This is verified through discrete element simulations that highlight the fundamental role played by the elastic energy stored within a granular material before a bifurcation occurs, which also corresponds to a minimization of the entropy production. The analysis suggests a new interpretation of the intriguing shear banding phenomenon as a bifurcation with the emergence of ordered dissipative structures germane to nonequilibrium thermodynamics of open systems.

## KEYWORDS:

Granular materials, Shear band, Instability, Entropy production, Bifurcation, Break in symmetry, Elastic energy, Dissipative structure, Microstructure.
# 1. INTRODUCTION

The term "dissipative structure" was coined in 1969 by Ilya Prigogine to emphasize that far from thermodynamic equilibrium, when systems are traversed by flow of matter and energy, structuring processes leading to spontaneous organization can occur within these systems (Glansdorff and Prigogine, 1971; Nicolis and Prigogine, 1977). The association between structure and dissipation is apparently paradoxical since the former evokes order while the latter calls for waste, disorder, and degradation. This is because the second principle of thermodynamics, which deals with dissipative processes, has traditionally been related to the irreversible evolution of a system toward a state of equilibrium identified as the state of maximum disorder where all the usable energy of the system has been degraded. However, the discovery of dissipative structures reveals that those irreversible processes, far from equilibrium, can also play a constructive role and become a source of order (Prigogine and Lefever, 1975; Nicolis and Prigogine, 1977; Lefever, 1978).

It is worth recalling that a given system is said to be in thermodynamic equilibrium under constant temperature when any fluctuation that takes place stays bounded and eventually vanishes. When a system is in thermodynamic equilibrium<fnref n="1" />, it is automatically in mechanical equilibrium which is a less restrictive equilibrium type. 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. Conversely, when a system is not in mechanical equilibrium, it is out of thermodynamic equilibrium. In that case, any fluctuation developing within the system (that can be spontaneous or imposed) is likely to grow, making it evolve in an irreversible way 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. In this regard, the notion of time and space scales is fundamental, as discussed in a review paper by Veveakis and Regenauer-Lieb (2015). In accordance with the Fluctuations Theorems (Evans and Searle, 2002), for systems at equilibrium, fluctuations vanish through dissipative mechanisms occurring at elementary (microscopic) scales over small time scales. Dissipative processes initially stay uncorrelated with each other within the system. When dissipative mechanisms extend and grow, they become correlated with a marked incidence on the behaviour of the system. In such a case, the

<footnote n="1">¹ In the following, when no confusion is possible, the word “equilibrium” used alone will refer to “thermodynamic equilibrium”</footnote>
system enters progressively an out-of-equilibrium regime where internal fluctuations grow to
eventually lead to an irreversible process from a macroscopic point of view both in space and
over time.

Classical thermodynamics has always included the antagonistic concepts of order and disorder.
Following Boltzmann's illuminating contribution, the ordered macroscopic state is known to
be a rare state, i.e., when huge numbers of the elementary units of a complex system are found
in a limited number of microscopic configurations such as molecules restricted to move in the
same direction. By contrast, a disordered state is achieved when a huge number of
microscopic configurations is necessary to describe the state of the system as, for example,
when molecules move incoherently in all possible directions (Parisi and Sourlas, 2002). The
classical laws of thermodynamics refer to these "disordered" states through entropy, and the
system at equilibrium continuously fluctuates around average values determined by these laws
in the context of the mean field theory. It is of paramount importance to note that the second
principle of thermodynamics states that, close to the equilibrium, any fluctuation which brings
this system out of the state of equilibrium should vanish. In short, the disorder is 'stable'.

However, when a macroscale system is brought far from equilibrium (out-of-equilibrium
state), disorder is no longer 'stable'. Fluctuations at the micro scale, instead of vanishing, can
be amplified. The system then bifurcates towards a new response mode, where a gain in order
(as a restriction of the number of reachable microscopic configurations) can be observed
together with a new topological organization of the energy dissipation. A well-known
example is the Rayleigh-Bénard cells, that form within a liquid layer whose lower surface is
heated (Bénard, 1901; Rayleigh, 1916; Chandrasekhar, 1961). Such a liquid contains billions
and billions of molecules which, instead of following a disordered motion, organize
themselves into structured macroscopic cells. In systems where chemical reactions take place,
the deviation from equilibrium is often not large enough for dissipative structures to appear.
But on the other hand, when so-called non-linear reactions (e.g., autocatalysis and mutual
catalysis) are present, fluctuation amplifications can occur and lead to very diverse chemical
patterns (Turing, 1952), including among others the famous chemical clock pattern
(Bélousov-Zhabotinsky chemical reaction; Hudson and Mankin, 1981). In brief, a variety of
bifurcation modes can appear, making the system no longer homogeneous, with a spatial or a
time structuration that gives the system intrinsic dimensions.
It is worth noting that the discovery of dissipative structures has a significance that goes beyond physics, by addressing fundamental issues such as self-organization in biology, life and the becoming of Nature (Prigogine and Stengers, 1979). In fact, dissipative structures seem to emerge in complex systems, that are known to give rise to a wide spectrum of emerging properties (Damper, 2000; Aziz-Alaoui and Bertelle, 2009). Among complex systems, granular materials have attracted a marked interest in the past few years in materials sciences. Granular materials are involved in a variety of engineering purposes, such as pharmaceutical engineering, food particle storage, geological mass-driven hazards, and civil engineering. It was highlighted that granular materials exhibit a wide spectrum of emergent features that have been intensively investigated (see for example Nicot and Darve, 2007a and 2007b; Tordesillas, 2007; Tordesillas and Muthuswamy, 2009; Walker and Tordesillas, 2010). As a manifestation of complexity, several properties are observed at the specimen scale (also referred to as the macro scale), whereas those properties are absent at smaller (micro) scales, namely the particle scale or even larger scales involving a few grains. For example, the existence of a plastic flow rule that can be regular or irregular (Nicot and Darve, 2007a and 2007b), or the existence of a bifurcation domain in which a variety of failure modes can be encountered (Bigoni, and Hueckel, 1991; Petryk, 1993; Nicot and Darve, 2011) are some of the salient constitutive properties that can be observed on the granular assembly scale.

As for the existence of a bifurcation domain, the sudden appearance of a shear band within granular materials when sheared (as for example during the conventional drained triaxial compression test – axial compression under constant lateral stress; Desrues and Chambon, 2002), patterning of the material in several structured zones, has received much attention over the past decades. 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 fundamental reasons why a specimen along a continuous loading can spontaneously lose its (kinematic and microstructural) homogeneity to give way to a structured pattern remain largely mysterious. Taking advantage of the background developed by Hill (1962) and Mandel (1966) in relation to the occurrence of acceleration wave propagation discontinuities (see also Regenauer-Lieb et al., 2021 for a very interesting extension), the occurrence of such events can also be modelled by the Rice-Mandel criterion (Rice, 1975; Rudnicki and Rice, 1975). This criterion was recognized later as a subset of the broader second-order work criterion introduced in the middle of the past century by Hill (1958). More recently, the second-order work theory was formalized further (Nicot et al.,
2012), by establishing a direct relation between the sudden occurrence of an outburst in kinetic energy, and the vanishing of the (internal) second-order work corresponding to the non-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 operative (effective) failure of the material can only occur according to the loading direction applied, and to the loading control adopted (Wan et al., 2017).

Granular plasticity, through mesostructural rearrangements, was recognized to be a key ingredient for the emergence of potential material instabilities (Wautier et al., 2018). However, the link between material instabilities and the framework of the out-of-equilibrium thermodynamics remains to be investigated in the light of the pioneering contributions of Bazant (1988 and 1989) who was the first to relate in a proper formalism the second-order work to entropy-based thermodynamic considerations. Thus, the purpose of this paper is to address this issue by deriving an equation relating a potential increase in kinetic energy within the system (marking a transition from a quasi-static toward a dynamic response regime) to the second-order time derivative of the total entropy of the system. This result augurs well with the fundamental minimum entropy production theorem established in the past century by Prigogine (Glansdorff and Prigogine; 1954, 1963 and 1964). The findings are then discussed within the framework of granular materials subjected to shear loading conditions and offer a proper conceptual foundation to interpret the development of shear bands as the emergence of dissipative structures.

Throughout this paper, only rate-independent materials are considered. Soil mechanics convention will be used with strains taken as positive in compression. Vectors, matrices, and tensors are represented in either bold face or index notation. Einstein notational convention is used such that repeated indices mean summation. Moreover, time and spatial derivatives of any term $\psi$ will be distinguished by denoting $\delta\psi$ the incremental variation with time of $\psi$ (defined as the product of the particulate derivative $\dot{\psi}$ by the infinitesimal time increment $\delta t$) with respect to a given frame, and by denoting $d\psi$ the spatial differential of $\psi$ with respect to the spatial coordinates $x_i$, with $d\psi = \frac{\partial\psi}{\partial x_i} dx_i$. Here, it must be cautioned that the common
notation $\dot{\psi}$ should apply only to state functions that are not path dependent. For example, if $\dot{\psi}$ corresponds to the rate of heat and work supplied to the system, it cannot be interpreted as a state function and hence deals with inexact differential (Teixeira-Dias, 2017). In this paper, as the uncertainty related to the path dependent integration of inexact differentials is bounded by extrema of entropy production, it is not necessary to introduce a special notation for inexact differentials. Thus, the same notations will be used for path-dependent terms such as heat in heat-related terms, or internal and external entropy. Moreover, whenever $\psi$ is a force or a stress term, the rate of $\psi$ should be understood as a material, objective derivative such as the Jaumann rate whose definition is frame independent. Here again, in order not to obfuscate the notations, the same overdot symbol will be used throughout the paper.

## 2. THERMODYNAMICS BACKGROUND

The incremental variation of the total energy $E^{tot}$ of a given closed system can be expressed as in Eq. (1), where $E_c$ denotes the kinetic energy of the system, $E^{int}$ is the internal energy, $P^{ext}$ is the external power of forces applied to the system, and $\dot{Q}$ is the heat power provided to the system:

$$ \delta E^{tot} = \delta(E^{int} + E_c) = (P^{ext} + \dot{Q}) \delta t \quad (1) $$

The external power comprises a volume term involving body forces $\gamma = \rho g$ acting within the volume $V$ ($\rho$ being the mass per volume of the material, and $g$ the acceleration due to gravity), and a surface term with tractions $f$ applied to the system boundary $\partial V$, as follows:

$$ P^{ext} = \int_V \gamma_i \dot{u}_i \, dv + \int_{\partial V} f_i \dot{u}_i \, ds \quad (2) $$

Likewise, the heat power reads:

$$ \dot{Q} = \int_V \rho \dot{r} \, dv - \int_{\partial V} \dot{q}_i n_i \, ds \quad (3) $$

where the first integral corresponds to the radiated heat within the system body, and the second integral accounts for the boundary conduction heat transfer with a heat flux $\dot{q}$.
Invoking the theorem of kinetic energy that reads:

$$ \delta E_c = \int_{\partial V} f_i \delta u_i ds - \int_V \sigma_{ij} \frac{\partial(\delta u_i)}{\partial x_j} dv \quad (4) $$

and combining with Eq. (1) finally give:

$$ \dot{E}^{\text{int}} = \int_V \rho \dot{r} \, dv + \int_V \sigma_{ij} \frac{\partial \dot{u}_i}{\partial x_j} \, dv - \int_{\partial V} \dot{q}_i n_i \, ds \quad (5) $$

Equation (5) expresses the internal energy $E^{\text{int}}$ within the system in an integral form. Noting $E^{\text{int}} = \int_V \rho e \, dV$, where e denotes the specific internal energy, the equation can be written with a local formulation as follows:

$$ \rho \dot{e} = \rho \dot{r} + \sigma_{ij} \frac{\partial \dot{u}_i}{\partial x_j} - \frac{\partial \dot{q}_i}{\partial x_i} \quad (6) $$

The resulting Eq. (6) corresponds to the first law of thermodynamics.

Under the external thermal and mechanical loading, the entropy S of the system should evolve. This evolution can be split into two parts, as follows:

$$ \dot{S} = \dot{S}^{\text{ext}} + \dot{S}^{\text{int}} \quad (7) $$

where $\dot{S}^{\text{ext}}$ denotes the rate of entropy exchanged with the external medium, and $\dot{S}^{\text{int}}$ denotes the rate of internal entropy. It then follows that:

$$ \dot{S}^{\text{ext}} = \frac{\dot{Q}}{\theta} = \int_V \rho \frac{\dot{r}}{\theta} dv - \int_{\partial V} \frac{\dot{q}_i}{\theta} n_i ds \quad (8) $$

where $\theta$ is the thermodynamic temperature.

Furthermore, the second law of thermodynamics imposes that $\dot{S}^{\text{int}} \ge 0$ (zero in case of reversible processes). This term can be expressed as $\dot{S}^{\text{int}} = \int_V \frac{\Phi}{\theta} dv$, where $\Phi$ denotes the total dissipation rate taking place within the system. Noting $S = \int_V \rho \eta \, dv$, where $\eta$ is the specific entropy, it can be shown that:
$$ \Phi = \rho \theta \dot{\eta} - \rho \dot{r} + \frac{\partial \dot{q}_i}{\partial x_i} - \frac{\dot{q}_i}{\theta} \frac{\partial \theta}{\partial x_i} \quad (9) $$

where $\Phi^{\text{int}} = \rho \theta \dot{\eta} - \rho \dot{r} + \frac{\partial \dot{q}_i}{\partial x_i}$ stands for the intrinsic dissipation rate, whereas $\Phi^{\text{th}} = -\frac{\dot{q}_i}{\theta} \frac{\partial \theta}{\partial x_i}$ represents the conduction heat dissipation rate. Introducing the Helmholtz specific free energy $\Psi = e - \theta \eta$, and taking advantage of the first law of thermodynamics as expressed in Eq. (6), yield:

$$ \Phi^{\text{int}} = -\rho (\dot{\Psi} + \dot{\theta} \eta) + \sigma_{ij} \frac{\partial \dot{u}_i}{\partial x_j} \quad (10) $$

The second law of thermodynamics states that:

$$ -\rho (\dot{\Psi} + \dot{\theta} \eta) + \sigma_{ij} \frac{\partial \dot{u}_i}{\partial x_j} - \frac{\dot{q}_i}{\theta} \frac{\partial \theta}{\partial x_i} \geq 0 \quad (11) $$

This relation is known as the Clausius-Duhem inequality.

The intrinsic dissipation, as expressed in Eq. (10), corresponds to any process for which dissipation is not solely thermal. This will be explored further in the next section, by addressing the case where a material system is subjected to a mechanical loading applied on its boundary under isothermal conditions.

# 3. EVOLUTION OF A MATERIAL SYSTEM UNDER EXTERNAL FORCES

## 3.1 Preamble

Throughout this section, a system consisting of a volume $V_o$ (with $\Gamma_o = \partial V_o$) of a given material, initially in a configuration $C_o$ is considered. After a loading history, the system is in a strained configuration $C$ and occupies a volume $V$ (with $\Gamma = \partial V$), in mechanical equilibrium under a prescribed external loading. This loading is controlled by specific boundary static or kinematic parameters, referred to as the control parameters. It will be assumed hereafter that no chemical reactions take place within the system, and that no matter is being added to or removed from the system, i.e., a closed system. For the sake of simplicity,
the thermal capacity of the material will be assumed to be very large, so that the temperature of the system is virtually constant. This assumption does not weaken the scope of the approach, in the sense that we mainly investigate the contribution of mechanical processes responsible for the intrinsic dissipation (Collins and Houlsby, 1997).

We introduce the transformation $\chi$ relating each material point $\mathbf{x}$ of the current configuration $C$ to the corresponding material point $\mathbf{X}$ of the initial configuration $C_o$. The continuity of the matter ensures that $\chi$ is bijective, i.e., the transformation is invertible with a one-to-one correspondence. One consequence is that the Jacobian $J$ of the tangent linear transformation (deformation gradient) $\mathbf{F}$, with $F_{ij} = \partial x_i / \partial X_j$, and its determinant being strictly positive. The displacement field $\mathbf{u}$ of material points between both initial and current configurations is defined by the relation $\mathbf{x} = \mathbf{X} + \mathbf{u}$, with $\mathbf{u}$ being a function of $\mathbf{X}$.

In what follows, we investigate the stability of the system at a given thermodynamic state, under a given externally prescribed loading.

## 3.2 The second-order work equation

Let us consider the material system at equilibrium at time $t$. If a change in the loading parameters is applied at time $t$ over a small duration $\Delta t$, the evolution of the system is shown to follow in Lagrangian formulation (See Appendix 1):

$$E_c(t + \Delta t) = E_c(t) + \frac{(\Delta t)^2}{4} \left( \int_{\Gamma_o} \dot{f}_i \dot{u}_i ds_o - \int_{V_o} \dot{\Pi}_{ij} \frac{\partial \dot{u}_i}{\partial X_j} dv_o \right) \quad (12)$$

where $\Pi$ corresponds to the Piola-Kirchhoff stress tensor. Following Murnaghan's first attempts and Hill's pioneering contribution (Murnaghan, 1944; Hill, 1958), the two terms inside brackets in Eq. (12) involve second order works herein called $W_2^{ext}$ and $W_2^{int}$ due to the loading variables controlled by the external world on the boundary of the system, and the intrinsic constitutive properties of the material at hand, respectively. Further details can be found in Nicot and Darve, 2007; Nicot et al., 2012; Wan et al., 2017.
It thus transpires from Eq. (12) that the incremental evolution of the kinetic energy of the system within a finite time interval $[t; t + \Delta t]$ is governed by the difference between the external and internal second-order works. If we further assume that the kinetic energy of the system is zero at time $t$, namely $E_c(t) = 0$,

$$E_c(t + \Delta t) = \frac{(\Delta t)^2}{4} (W_2^{\text{ext}}(t + \Delta t) - W_2^{\text{int}}(t + \Delta t)) \quad (13)$$

This equation indicates that the kinetic energy of the system may abruptly increase from an equilibrium configuration, according to the competition between the external loading applied to the system (term $W_2^{\text{ext}}$), and the capacity of the system to adapt against this external loading through its constitutive behaviour responsible for internal stress and strain fields (term $W_2^{\text{int}}$). This is also reminiscent of the zero-acceleration wave condition first proposed by Hill (1962). In such eventuality, the system bifurcates from a quasi-static regime towards a dynamical regime. Conversely, when the system is able to sustain the external loading, both terms $W_2^{\text{ext}}$ and $W_2^{\text{int}}$ are equal, and Eq. (13) indicates that no increase in kinetic energy occurs. The system then stays in a quasi-static regime over the time increment $\Delta t$ such that $E_c(t + \Delta t) = E_c(t) = 0$.

To illustrate this theoretical result, one can refer to the classical drained triaxial laboratory test in geomechanics to explore the constitutive properties of geomaterials. During this test, a material sample is first isotropically compressed under a prescribed confining pressure $\Pi_o$. Then, the lateral pressure exerted onto the lateral walls is kept constant ($\Pi_o$), while the upper platen is kinematically controlled with a downward motion at a constant displacement rate. In short, the axial (direction 1) strain rate $\dot{F}_1$ is positive and constant, while the lateral (directions 2 and 3) pressures are kept constant: $\Pi_2 = \Pi_3 = \Pi_o$. For an initially dense specimen, a typical stress-strain response curve, as depicted in Fig. 1, can be obtained.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Typical stress-strain response curve of a dense soil specimen along a drained triaxial loading path.</figcaption></figure>

The axial stress $\Pi_1$ monotonously increases until reaching a peak ($\hat{\Pi}_1$), then decreases until approaching a steady-like, asymptotic regime. This latter regime is known as the Critical State regime, with no more volume change nor stress evolution. Close to the peak stress during the hardening regime, it is often observed the development of a localized kinematic pattern consisting of a single (or multiple) band(s) crossing the specimen from bottom to top. This pattern is referred to as a shear band and can have various topologies with sometimes multiple reflections on the sides according to various factors. The latter include: specimen geometry, presence of defects or heterogeneities, loading conditions, as well as boundary conditions imposed by the loading platens as illustrated in Fig. 2.

<figure><img src="image_3.png" /><figcaption>Fig. 2. Example of shear band formation within a geomaterial specimen subjected to a drained triaxial loading; after Desrues and Ando, 2015.</figcaption></figure>
If at the peak stress, the loading control was switched from a kinematical control to a statical
one (for example, a small axial stress increase could be imposed instead of prescribing a
constant axial strain rate), the specimen would abruptly collapse, marked by violent outbursts
in kinetic energy and large deformations. This is in perfect accord with the theoretical
predictions of Eq. (13) (Nicot et al., 2012). At the peak stress, Fig. 1 indicates that the internal
stress should decrease as the material has exhausted all mobilization of mechanical strength
against applied loads. In fact, if external loading were to be next increased, an imbalance
would occur between what is imposed on the boundary of the system and what the
constitutive properties of the material allow in terms of internal stress. As a result, a (time)
bifurcation occurs from a quasi-static response to a dynamic one (Nicot et al., 2012; Wan et
al., 2017).

On the other hand, if the loading control were kept unchanged, no outburst in kinetic energy
would be observed, but a (spatial) kinematic bifurcation would occur. The system dissipates
energy in terms of the emergence of a patterned deformation. The incremental strain field is
no longer homogeneous (as roughly observed before the peak) — the incremental strain field
is heterogenous, with one or more shear bands crossing the specimen.

It is worth noting that both loose and dense specimens will behave in totally different ways
under a drained triaxial loading. The dense case gives rise to a peak stress followed by a
softening regime corresponding to the formation of a shear band, whereas the loose one will
experience a hardening regime until reaching an asymptotic plateau in the absence of any
shear localization. Also, heterogeneities in local density or void distribution that develop
during loading history are a determinant for shear band branching in the form of various mode
switchings and the associated entropy rate productions.

The above phenomenon can also be verified from numerical discrete simulations (based on a
discrete element method; see for example Cundall and Strack, 1979), as shown in Fig. 3. For
comparison, two specimens of different initial density are considered: a dense specimen on
the left, and a loose specimen on the right in Fig. 3. The distribution of the incremental
deviatoric strain $\varepsilon_d$ (see footnote 2) reveals a clear heterogeneous pattern in the case of a dense specimen, with a shear band crossing the specimen roughly along a diagonal. The shear band appears as a narrow domain in which $\varepsilon_d$ is much more intense than in the rest of the specimen. Interestingly, this is no longer the case when the initial porosity is increased.

<figure><img src="image_4.png" /><figcaption>Fig. 3. Distribution of incremental deviatoric strain within a granular assembly subjected to a drained biaxial loading. A shear band has clearly developed within the dense specimen after the stress peak is reached (a), whereas the distribution remains much more random and roughly homogeneous within the loose specimen (b).</figcaption></figure>

The distribution has been computed for an axial strain increment $\delta\varepsilon_1 = 10^{-3}$; (Liu et al., 2018).

The purpose of what follows is to discuss this loss of homogeneity as a kinematic bifurcation in relation to thermodynamics.

## 3.3 The second-order work approach in the light of thermodynamics

### 3.3.1 General setting

<footnote n="2">² In 2D conditions, $\varepsilon_d = \sqrt{\frac{1}{2}\left(\delta\varepsilon_{ij} - \frac{1}{2}(\delta\varepsilon_{11} + \delta\varepsilon_{22})I_{ij}\right)\left(\delta\varepsilon_{ij} - \frac{1}{2}(\delta\varepsilon_{11} + \delta\varepsilon_{22})I_{ij}\right)}$, where $I_{ij}$ is the identity matrix</footnote>
In the example discussed above (Fig. 3), the occurrence of a shear band at the peak stress could be regarded as a structuring event within the material sample. Focusing on granular materials, as will be done in section 4, the shear band coincides with a gain in order in the sense that the kinematic pattern is no longer random as it might have been before the peak stress for a dense specimen (Fig. 3a) or in a loose specimen (Fig. 3b), but organized within specific spatial zones. The objective of this section is to analyze this transition by invoking the thermodynamic equations coupled with the second-order work equation (Eq. 13).

As the thermodynamic temperature within the material system is assumed to be constant (isothermal conditions) and homogeneous, the dissipation potentials written in Lagrangian formulation are:

$$ \Phi_o^{th} = 0 \qquad (14) $$

$$ \Phi_o^{int} = -\rho_o \dot{\Psi} + \Pi_{ij} \frac{\partial \dot{u}_i}{\partial X_j} = -\rho_o \dot{\Psi} + \Pi_{ij} \dot{F}_{ij} \qquad (15) $$

During loading, the specimen stores a part of the external work into elastic energy at contacts which corresponds to the free energy in the sense that it is usable energy available in the system that can be later retrieved and transformed. At the contact scale, this elastic energy is reversible; however, this is no longer true at the specimen scale where geometrical effects come into play. The latter are also responsible for the strong coupling between both elastic and plastic mechanisms. The reader can refer to Nicot and Darve (2006) for a discussion on this issue. In essence, the free energy can also be defined as the maximum stored energy that can be used as elastic work exchanged with the outside or dissipated within the system into plastic work. Even if the local elastic strain energy, on the contact scale, is recoverable, this is no longer true on the macroscopic scale: the macroscopic elastic strain energy is not totally recoverable, which led to the notion of frozen energy (Collins and Muhunthan, 2003; Collins, 2005; Nicot and Darve, 2006).

The internal power $\dot{W} = \Pi_{ij} \dot{F}_{ij}$ can therefore be split in an elastic part $\dot{W}^e$ and a plastic counterpart $\dot{W}^p$:

$$ \dot{W} = \dot{W}^e + \dot{W}^p \qquad (16) $$
The classical decomposition $\dot{F}_{ij} = \dot{F}_{ij}^e + \dot{F}_{ij}^p$ ensues from Eq. (16), with $\dot{W}^e = \Pi_{ij} \dot{F}_{ij}^e$ and $\dot{W}^p = \Pi_{ij} \dot{F}_{ij}^p$. Thus, on the boundary of the system, the displacements can also be decomposed as: $u_i = u_i^e + u_i^p$, where $u_i^e$ corresponds to the elastic part of the displacement taking place on the boundary, and $u_i^p$ to its plastic counterpart.

Noting that the rate of the Helmholtz free energy is $\rho_o \dot{\Psi} = \Pi_{ij} \dot{F}_{ij}^e$, the internal dissipation potential in Eq. (15) reduces to $\Phi_o^{int} = \Pi_{ij} \dot{F}_{ij}^p$, which corresponds to the intrinsic dissipation rate solely emanating from all contacts, which is always positive, i.e.

$$ \Phi_o^{int} = \Pi_{ij} \dot{F}_{ij}^p \geq 0 \qquad (17) $$

which follows the second law of thermodynamics.

As $\dot{S}^{int} = \int_{V_o} \frac{\Phi_o^{int}}{\theta} dv_o = \int_{V_o} \rho_o \dot{\eta}^{int} dv_o$, it can also be written that:

$$ \rho_o \theta \dot{\eta}^{int} = \Pi_{ij} \dot{F}_{ij}^p \qquad (18) $$

where $\dot{\eta}^{int}$ corresponds to the specific internal entropy within the system.

Furthermore, the external entropy flow can be related to the dissipation taking place on the boundary of the system, namely:

$$ \theta^{ext} \dot{S}^{ext} = -\int_{\Gamma_o} f_i \dot{u}_i^p ds_o \qquad (19) $$

where $\theta^{ext} = \theta$, as isothermal conditions are assumed.

Finally, it can be written that:

$$ \dot{S} = \int_{V_o} \frac{\Pi_{ij} \dot{F}_{ij}^p}{\theta} dv_o - \int_{\Gamma_o} \frac{f_i \dot{u}_i^p}{\theta} ds_o \qquad (20) $$

### 3.3.2 A second-order entropy equation

Let us now assume that at time $t$, the system is in mechanical equilibrium under a prescribed external loading. An incremental loading distribution ($\Delta\vec{f}$ or $\Delta\vec{u}$) is smoothly imposed on the
boundary of the system from time $t$ to time $t+\Delta t$, making the system evolve. The first order time derivatives can therefore be supposed equal to zero at time $t$. Let us examine the different terms at time $t+\Delta t$. Equation (18) reads:

$$ \rho_o \theta \dot{\eta}^{\text{int}}(t+\Delta t) = \Pi_{ij}(t+\Delta t) \dot{F}_{ij}^p(t+\Delta t) \quad (21) $$

Adopting a first-order time series expansion, the following relations hold:

$$ \dot{\eta}^{\text{int}}(t+\Delta t) = \dot{\eta}^{\text{int}}(t) + \frac{\Delta t}{2} (\ddot{\eta}^{\text{int}}(t) + \dot{\eta}^{\text{int}}(t+\Delta t)) = \frac{\Delta t}{2} \dot{\eta}^{\text{int}}(t+\Delta t) \quad (22) $$

$$ \Pi_{ij}(t+\Delta t) = \Pi_{ij}(t) + \frac{\Delta t}{2} (\dot{\Pi}_{ij}(t) + \dot{\Pi}_{ij}(t+\Delta t)) = \Pi_{ij}(t) + \frac{\Delta t}{2} \dot{\Pi}_{ij}(t+\Delta t) \quad (23) $$

$$ \dot{F}_{ij}^p(t+\Delta t) = \dot{F}_{ij}^p(t) + \frac{\Delta t}{2} (\ddot{F}_{ij}^p(t) + \dot{F}_{ij}^p(t+\Delta t)) = \frac{\Delta t}{2} \ddot{F}_{ij}^p(t+\Delta t) \quad (24) $$

Thus, by virtue of Eqs. (22)-(24), and after rearranging the terms, Eq. (21) gives:

$$ \rho_o \theta \dot{\eta}^{\text{int}}(t+\Delta t) = \Pi_{ij}(t) \ddot{F}_{ij}^p(t+\Delta t) + \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}^p(t+\Delta t) \quad (25) $$

which also can be written as:

$$ \rho_o \theta \dot{\eta}^{\text{int}}(t+\Delta t) = \Pi_{ij}(t) \ddot{F}_{ij}^p(t+\Delta t) + \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}(t+\Delta t) - \dot{\Pi}_{ij}(t+\Delta t) \ddot{F}_{ij}^e(t+\Delta t) \quad (26) $$

Thus, combining with Eq. (12), Eq. (26) yields:

$$ E_c(t+\Delta t) = \frac{(\Delta t)^2}{4} \left( \int_{\Gamma_o} \dot{f}_i(t+\Delta t) \dot{u}_i(t+\Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}^e(t+\Delta t) dv_o + \int_{V_o} \Pi_{ij}(t) \ddot{F}_{ij}^p(t+\Delta t) dv_o - \int_{V_o} \rho_o \theta \dot{\eta}^{\text{int}}(t+\Delta t) dv_o \right) \quad (27) $$

Furthermore, using Green's formula one gets:

$$ \int_{V_o} \Pi_{ij}(t) \ddot{F}_{ij}^p(t+\Delta t) dv_o = \int_{\Gamma_o} \Pi_{ij}(t) N_j \ddot{u}_i^p(t+\Delta t) ds_o - \int_{V_o} \frac{\partial \Pi_{ij}(t)}{\partial X_j} \ddot{u}_{ij}^p(t+\Delta t) dv_o \quad (28) $$

As the system is supposed to be at equilibrium at time $t$, it follows that $\frac{\partial \Pi_{ij}(t)}{\partial X_j} = 0$ and $\Pi_{ij}(t) N_j = f_i(t)$. Thus, Eq. (27) can be rewritten as:
$$E_c(t+\Delta t) = \frac{(\Delta t)^2}{4} \left( \begin{array}{l} \displaystyle \int_{\Gamma_o} \dot{f}_i(t+\Delta t) \dot{u}_i(t+\Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}^e(t+\Delta t) dv_o + \\ \displaystyle \int_{\Gamma_o} f_i(t) \ddot{u}_i^p(t+\Delta t) ds_o - \int_{V_o} \rho_o \theta \ddot{\eta}^{\text{int}}(t+\Delta t) dv_o \end{array} \right) \quad (29)$$

Furthermore, Eq. (19) gives:

$$\theta \dot{S}^{\text{ext}}(t+\Delta t) = -\int_{\Gamma_o} f_i(t+\Delta t) \dot{u}_i^p(t+\Delta t) ds_o \quad (30)$$

which can be rewritten as:

$$\theta \ddot{S}^{\text{ext}}(t+\Delta t) = -\int_{\Gamma_o} f_i(t) \ddot{u}_i^p(t+\Delta t) ds_o - \int_{\Gamma_o} \dot{f}_i(t+\Delta t) \dot{u}_i^p(t+\Delta t) ds_o \quad (31)$$

Combining Eqs. (29) and (31) finally yields:

$$E_c(t+\Delta t) = \frac{(\Delta t)^2}{4} \left( \begin{array}{l} \displaystyle \int_{\Gamma_o} \dot{f}_i(t+\Delta t) \dot{u}^e(t+\Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}^e(t+\Delta t) dv_o \\ \displaystyle -\theta \ddot{S}^{\text{ext}}(t+\Delta t) - \int_{V_o} \rho_o \theta \ddot{\eta}^{\text{int}}(t+\Delta t) dv_o \end{array} \right) \quad (32)$$

By virtue of Eq. (20), Eq. (35) also reads:

$$E_c(t+\Delta t) + \frac{(\Delta t)^2}{4} \theta \ddot{S}(t+\Delta t) = \frac{(\Delta t)^2}{4} \left( \begin{array}{l} \displaystyle \int_{\Gamma_o} \dot{f}_i(t+\Delta t) \dot{u}_i^e(t+\Delta t) ds_o \\ \displaystyle - \int_{V_o} \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}^e(t+\Delta t) dv_o \end{array} \right) \quad (33)$$

It is herein emphasized that in absence of dissipative processes, the existence of an elastic potential precludes any outburst of kinetic energy except for the case of geometrical instabilities. Thus, according to Eq. (13), both internal and external second-order works are equal, which reads:

$$\int_{\Gamma_o} \dot{f}_i(t+\Delta t) \dot{u}_i^e(t+\Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t+\Delta t) \dot{F}_{ij}^e(t+\Delta t) dv_o = 0 \quad (34)$$

Combining with Eq. (33) finally gives:

$$-\theta \ddot{S}(t+\Delta t) = W_2^{\text{ext}}(t+\Delta t) - W_2^{\text{int}}(t+\Delta t) = \frac{4}{(\Delta t)^2} E_c(t+\Delta t) \quad (35)$$

which highlights the key role of the imbalance between the internal mechanisms within the system related to material constitutive properties ($W_2^{\text{int}}$ term) and external loading ($W_2^{\text{ext}}$ term).
term). This imbalance triggers the emergence of an effective failure characterized by an outburst in kinetic energy and a decrease in the rate of the total entropy. The system is no longer able to store, nor dissipate all the input energy provided by the external loading.

The occurrence of an effective failure is marked by an increase in kinetic energy from zero value ($E_c(t) = 0$) to a non-zero ($E_c(t + \Delta t) > 0$), turning the mechanical regime from quasi-static to dynamic. According to Eq. (35), this bifurcation in the mechanical regime is equivalently signaled by the condition:

$$ \ddot{S}(t + \Delta t) \le 0 \qquad (36) $$

In combining with Eq. (13), it can also be derived from Eq. (33), that:

$$ \theta \ddot{S}(t + \Delta t) + \int_{\Gamma_o} \dot{f}_i(t + \Delta t) \dot{u}_i^p(t + \Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t + \Delta t) \dot{F}_{ij}^p(t + \Delta t) dv_o = 0 \qquad (37) $$

Recalling that $\dot{S} = \dot{S}^{\text{int}} + \dot{S}^{\text{ext}} = \int_{V_o} \frac{\Pi_{ij} \dot{F}_{ij}^p}{\theta} dv_o - \int_{\Gamma_o} \frac{f_i \dot{u}_i^p}{\theta} ds_o$, it can be thus inferred from Eq.

(37) that:

$$ \theta \dot{P}^{\text{int}}(t + \Delta t) = \int_{V_o} \dot{\Pi}_{ij}(t + \Delta t) \dot{F}_{ij}^p(t + \Delta t) dv_o \qquad (38) $$

where $P^{\text{int}} = \dot{S}^{\text{int}}$ refers to as the entropy production.

This can also be written locally as:

$$ \theta \Delta P^{\text{int}} \Delta t = \Delta \Pi_{ij} \Delta F_{ij}^p \qquad (39) $$

Noting that $\Delta \Pi_{ij} \Delta F_{ij} = \Delta \Pi_{ij} \Delta F_{ij}^e + \Delta \Pi_{ij} \Delta F_{ij}^p$, yields:

$$ \theta \Delta P^{\text{int}} \Delta t + q_2(\Delta F_{ij}^e) = \Delta \Pi_{ij} \Delta F_{ij} \qquad (40) $$

where $q_2^e(\Delta F_{ij}^e) = \Delta \Pi_{ij} \Delta F_{ij}^e$ is a necessarily a positive definite quadratic form.

Also, $q_2(\Delta F_{ij}) = \Delta \Pi_{ij} \Delta F_{ij}$ takes a quadratic form, but not necessarily positive definite, the sign of which can be positive or negative depending on whether: (a) the mechanical state lies within to the bifurcation domain, and (b) the incremental loading $\Delta F_{ij}$ (or $\Delta \Pi_{ij}$) belongs to the instability cone (see Wan et al., 2017, for a thorough review). Thus, as
$\theta \Delta P^{\text{int}} \Delta t \le \Delta \Pi_{ij} \Delta F_{ij}$, if $q_2(\Delta F_{ij})$ is negative, then $\Delta P^{\text{int}}$ is negative too (negative entropy production).

It can therefore be inferred from Eq. (40) that the entropy production $P^{\text{int}}$ decreases and tends to be minimum along an incremental loading direction belonging to the instability cone, i.e., such that $q_2(\Delta F_{ij}) < 0$

Finally, inequality (36) should be compared to the main result inferred by Prigogine in the context of out-of-equilibrium thermodynamics (Glansdorff and Prigogine, 1964 and 1971) which states that the negativity of the second-order derivative of the entropy is a necessary condition of instability around an equilibrium state. Fluctuations can increase without turning back to zero, transporting the system toward a new state that is different from the initial stationary state. When the system stays close to the equilibrium state, Eqs. (39) and (40) are reminiscent of the minimal entropy production theorem, which was key in the subsequent expansion toward the out-of-equilibrium thermodynamics (Prigogine, 1945 and 1947; Ziegler and Wehrli, 1987). In fact, once the system is brought out-of-equilibrium, it reorganizes itself in a way that extremizes the entropy production to reach a new, more stable equilibrium state (Ziegler, 1983; Dewar, 2005; Veveakis and Regenauer-Lieb, 2015).

In the context of the thermodynamics of irreversible processes, and particularly when dealing with diffusion or thermal mechanisms, instability occurs once a certain level of loading is reached. Indeed, Lefever (1978) writes that a certain critical amount of dissipation should have occurred for instability to develop. For example, Rayleigh-Bénard cells develop when the heating at the bottom of the fluid reaches a critical value of the Rayleigh number. Then, the system evolves abruptly toward a new regime giving rise to ordered cells. Such structures were denoted as *dissipative structures* by Prigogine, highlighting the fact that the system enters a more ordered regime to dissipate all the external energy provided to the system.

Extending this reasoning to the mechanical system considered throughout this section, we now turn to analytical argumentations. Under the effect of the external loading, the system stores a part of the external work in elastic energy, the remaining part being dissipated through plastic mechanisms. Before the plastic limit surface is reached, the system can become potentially unstable, as soon as the mechanical state belongs to the bifurcation
domain<fnref n="3" /> (Darve et al., 2004; Darve and Vardoulakis, 2005; Nicot et al., 2012; Wan et al., 2017). In such a situation, it was proved that the symmetric part of the constitutive operator is no longer elliptic, with at least one eigenvalue being negative. Thus, incremental loading directions exist, along which the second-order work is negative. For a suitable choice of control parameters (Nicot et al., 2012; Wan et al., 2017), the system can bifurcate from a quasi-static regime towards a dynamical one usually marked by a sudden outburst in kinetic energy. The essential notion of out-of-equilibrium state far from the equilibrium state, suggested by Prigogine, could therefore be compared to the concept of potential instability that occurs as soon as the mechanical state of the system belongs to the bifurcation domain. Such potential instabilities require that the constitutive behaviour is no longer conservative (plastic dissipation should take place), as the constitutive operator should admit at least one negative eigenvalue. Basically, as will be shown in the next section, the system approaches its limit capacity for storing elastic energy. At the stress peak, no more elastic energy is stored.

The onset of inertial mechanisms, detected by a local increase in kinetic energy, and
repeatedly observed along the softening regime during a drained triaxial test (section 3.2),
coincides with the formation of a fully developed shear band, crossing the material specimen
(as observed in Fig. 2 from lab experiments, or in Fig. 3 from discrete numerical simulations).
This shear band corresponds to a bifurcation from a homogeneous stress-strain pattern, toward
a structured one. The system is more ordered, with a narrow band where most of the
dissipative processes concentrate. The shear band can therefore be regarded as an example of
dissipative structure, in the sense suggested by Prigogine. As emphasized in several seminal
works, dissipative structures develop far from the equilibrium, marking the break in symmetry
of the initial homogenous configuration. The purpose of the next section is to explore with
more details this bifurcation, by considering the mechanical response of a granular assembly
simulated with a discrete numerical approach.

<footnote n="3">³ The bifurcation domain, for a given material, is a part of the stress space where the constitutive operator relating both incremental strain and stress is no longer positive-definite, and admits at least one nil or negative eigenvalue. When the stress state of a given material belongs to the bifurcation domain, it can be shown that a set of incremental stress (resp. strain) loading directions exist, along which the second-order work takes nil or negative values. This set of incremental loading directions is referred to as the instability cone. This cone belongs to the incremental stress (resp. strain) space.</footnote>
# 4. MICROMECHANICAL INVESTIGATION USING A DISCRETE APPROACH

## 4.1 General setting

A granular assembly composed of $N$ grains 'p', with $1 \le p \le N$, is now considered (Fig. 4). Each grain 'p' is in contact with $n_p$ other grains 'q'. Each grain 'p' belonging to the boundary $\partial V$ of the volume is subjected to an external force $\mathbf{F}^{\text{ext},p}$. It will be assumed that no torque is applied to the particles of $\partial V$:

$$ \forall p \in \partial V, \delta \mathbf{M}^{\text{ext},p} = \mathbf{0} \qquad (41) $$

In addition, the external loading will be considered equilibrating at any time:

$$ \sum_{p \in \partial V} \mathbf{F}^{\text{ext},p} = \mathbf{0} \qquad (42) $$

Under the effect of the external loading ($\mathbf{F}^{\text{ext},p}$, $p \in \partial V$), the granular assembly rearranges mainly through dissipative mechanisms (such as intergranular sliding) and marginally through grain deformation. In most discrete approaches, such as the Discrete Element Method (Cundall and Strack, 1979), the grain deformation is described by grain overlapping. Although this process contributes marginally to the global deformation, it is the main source of elastic energy storage. While the external loading is increased (the external work provided to the system is increasing), the elastic energy stored within the system through *grain overlapping* is increasing too.

<figure><img src="image_5.png" /><figcaption>Fig. 4. Definition of the granular assembly: Boundary of the REV and external forces.</figcaption></figure>
## 4.2 Analysis of the post-bifurcated regime from a DEM approach

The drained triaxial loading test as described in section 3.2 is herein used as a prototype example to numerically explore the thermodynamics aspects of shear banding. As such, the mechanical response of a dense specimen subjected to a drained biaxial loading was simulated using YADE DEM code (see Appendix 2 for more details) in 2D conditions with model parameters summarized in Table 1.

<table><caption>Table 1. Parameters and numerical values used during the simulations.</caption><thead><tr><td>Parameters</td><td>Numerical values</td></tr></thead><tbody><tr><td>Specimen initial height (m) × width (m)</td><td>1.5 × 1.0</td></tr><tr><td>Range of particles radii (m)</td><td>0.003 – 0.006</td></tr><tr><td>Initial porosity</td><td>0.153</td></tr><tr><td>Interparticle friction coefficient</td><td>0.5</td></tr><tr><td>Damping coefficient</td><td>0.01</td></tr><tr><td>Contact normal stiffness k<sub>n</sub> (N/m)</td><td>5 10<sup>8</sup></td></tr><tr><td>Contact tangential stiffness k<sub>t</sub> (N/m)</td><td>5 10<sup>8</sup></td></tr><tr><td>Confining pressure σ<sub>o</sub> (kPa)</td><td>400</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Evolution of the deviatoric stress (q) and of the volumetric strain (ε<sub>v</sub>) with axial strain. The dashed line A corresponds to the characteristic point (maximal contractancy), the dashed line B corresponds to the deviatoric stress peak, whereas the dashed line C corresponds to the beginning of the critical state regime.</figcaption></figure>
Fig. 5 depicts the computed deviatoric stress-strain-volumetric strain response which is typical for a dense sand. Characteristic points of the response are identified by lines A, B and C marking: maximum contractancy, peak deviatoric stress and beginning of the critical state regime, respectively. As the specimen is incrementally loaded, the granular assembly rearranges and adapts itself to external loading until the peak stress is reached.

<figure><img src="image_7.png" /><figcaption>Fig. 6. Evolution of both: (a) the external work, and (b) the elastic energy stored within the system along an axial compression path under constant lateral stress.</figcaption></figure>

In line with the previous sections, the energetics of the triaxial test is next explored in Fig. 6. As loads are being incrementally applied so that external work continuously increases in Fig. 6a, the system stores elastic energy $E^{el}$ at contacts that can be easily computed (see Appendix 2). It is found that the accumulation of stored elastic energy is limited as signaled by a peak value followed by a decrease, hinting that part of the elastic energy is being released and dissipated; see Fig. 6b. Interestingly, this limit in elastic energy storage is reached exactly at the peak stress when the shear band simultaneously develops (Fig. 7). Thus, the bifurcation of the system from a homogenous regime toward a heterogeneous one (corresponding therefore to a proper break in symmetry) is strongly governed by the maximal capacity of the system to store elastic energy (Sun et al., 2015).
<figure><img src="image_8.png" /><figcaption>Fig. 7. Distribution of incremental deviatoric strain ratio ($\delta\varepsilon_d / \delta\varepsilon_1$) within the granular assembly at the peak stress (point B, left) and at the beginning of the critical state regime (point C, right).</figcaption></figure>

The elastic energy stored within the system wanes after the peak stress to reach a non-zero plateau. This suggests an important feature characterizing the post-bifurcation critical state regime: the system is not able to store more elastic energy. This is perfectly in line with the recent findings of Liu et al. (2018, 2020 and 2022) who showed that shear band development during the critical state regime for dense granular specimens is marked by continuous microstructural rearrangements. Elastic energy continuously released from opening contacts is dissipated into plastic energy. Short and fast grain motions occur, associated with recurrent, short-range outbursts in kinetic energy. The nature of the critical state regime that emerges is therefore stationary only on an average sense at the macroscopic scale. From a microstructural point of view, the critical state consists of intense and continuous grain reorganization (Deng et al., 2022).

As far as dense specimens are concerned, the post-bifurcated regime beyond peak stress corresponds to a macroscopic increase in order. The specimen becomes structured with the occurrence of shear banding which divides it into several zones. Thus, an increase in disorder on the microscopic scale is observed through repeated and intense microstructural rearrangements (Parisi et al., 2017). Above all, the free energy is minimal during the post-bifurcated critical state regime. Hence, the critical state can be regarded as a proper state of equilibrium, where all the usable energy newly stored in the system is immediately degraded through a fast-fading memory process; see Deng et al., 2022. This is likely the reason why
such a regime is known to be so robust, acting as an attractor, irrespective of the initial fabric of the material (Deng et al., 2021).

Interestingly, the entropy production along the loading path can be determined from the relation:

$$ \theta P^{\text{int}} = \int_{v_o} \prod_{ij} \dot{F}_{ij}^p dv_o \qquad (43) $$

where the right-hand side integral corresponds to the plastic power, and the thermodynamic temperature $\theta$ is assumed to be constant (isothermal conditions). Thus, the evolution of the entropy production $P^{\text{int}}$ corresponds to that of the plastic power within the system.

The plastic power can be readily obtained by subtracting the elastic power from the external power provided to the system. It should be noted that it corresponds to the plastic dissipation pertaining to the frictional sliding occurring at the contacts between grains. This quantity is accessible from the DEM computations, at each time step. As the evolution of the plastic dissipation energy follows a noisy signal, the derivative of such a signal would inevitably amplify the noise, even though the amplitude of the fluctuations in plastic dissipation energy remains small with respect to the mean signal. For this reason, to avoid amplification of noise, the raw fluctuating plastic dissipation curve has been fitted piecewise using polynomials (see Appendix 3), and the entropy production curve subsequently obtained by direct analytical differentiation to get a smooth curve (Fig. 8).

As observed in Fig. 8, the entropy production reaches a peak at the end of the softening regime, when the critical state is reached (point C). After this peak, the entropy production decreases until reaching a constant value while the critical state continues to develop. As a result, the bifurcation of the system after the peak stress from a homogeneous state toward a heterogeneous one corresponds to a structuring of the material that minimizes the entropy production. Through a structuring pattern, the shear band will make it possible to optimize the plastic dissipation within the system, with a high level of dissipation within the shear band, whereas the dissipation level remains very small outside the band. Indeed, it has repeatedly been shown that sliding contacts are mainly located within the shear band (see for example Liu et al, 2020).
The shear band, as a structural medium, can therefore be regarded as a proper dissipative structure in the sense given by Prigogine (Prigogine and Lefever, 1968 and 1975; Glansdorff and Prigogine, 1971; Lefever, 1978). The emergence of this structure coincides with a gain in macroscopic ordering to optimize the level of dissipation within the system. It is worth emphasizing that the minimal value of entropy production is not reached at the beginning of the critical state, but only later. This suggests that while the microstructural reorganization after peak stress is responsible for a transition in the energy dissipation modalities, it is after all a long-lasting process that stabilizes after the critical state regime has started. This could be due to fluctuations in the shear band topology that are still operative at the beginning of the critical state regime.

Further discussions are necessary to link the microstate continua entropies to the critical state regime. In fact, Prigogine's minimal entropy production principle applies to near-equilibrium states. More specifically, this principle establishes that the steady state reached by a system corresponds to a minimal entropy production only if: (i) the system is close enough to the equilibrium (Onsager's reciprocity relations hold), and (ii) the steady state corresponding to the minimal entropy production is unique. Otherwise, when (ii) is denied (several steady states exist), the maximal entropy production principle applies (MEP), and the steady state effectively reached by the system (local minimal entropy production) is the one corresponding to the largest of the different minimal entropy productions. The extremal entropy production, i.e., minimal entropy production vs. maximal entropy production, is still a debatable question, although about to be clarified (Ziegler, 1983; Dewar, 2005; Veveakis and Regenauer-Lieb, 2015). Interestingly, Veveakis and Regenauer-Lieb (2015) discuss the extremal (min vs. max) entropy production diffusion mechanism across both space and time scales. In line with Veveakis and Regenauer-Lieb's conjecture, we hypothesize that the delay in getting the minimal entropy production value well after state B stems from a diffusion mechanism from the elementary grain contact scale to the whole assembly scale. Investigating the entropy production over the scales at hand (from the grain scale to the specimen scale, and the related time scales) is outside the scope of this work.
<figure><img src="image_9.png" /><figcaption>Fig. 8. Evolution of the entropy production (times the thermodynamic temperature considered constant) along an axial compression path under constant lateral stress. The dashed line A corresponds to the maximal contractancy state, the dashed line B corresponds to the deviatoric stress peak, and the dashed line C corresponds to the beginning of the critical state regime. To avoid amplification of noise, the raw fluctuating plastic dissipation curve has been fitted piecewise using polynomials, and the entropy production curve subsequently obtained by direct analytical differentiation.</figcaption></figure>

5. CONCLUSION

The paper intends to shed a new light on the bifurcation mechanism operating within granular materials by advocating the field of out-of-equilibrium thermodynamics. Juxtaposing the fundamental minimum entropy production theorem with the second-order work theory, a general equation linking the derivatives of the entropy of the system to the second-order work is inferred. It is therefore proven than the negativity of the second-order derivative of the entropy is a necessary condition of instability around a steady state. Fluctuations within the granular assembly can increase without returning to zero, transporting the system toward a new state that is different from the initial stationary state. This is reminiscent of the minimal entropy production theorem, which was key to the development of out-of-equilibrium thermodynamics. More recent contributions belonging to the post Prigogine’s era have established variational principles applied to the time derivative of the second law of
thermodynamics and have contributed significantly to a better understanding of dynamic
systems. Thus, Prigogine's minimal entropy production principle should be thought as a
subset of the more complex case with extremal entropy production involving various time and
length scales.

For verification purposes, the above framework of out-of-equilibrium thermodynamics is
presented through the discrete element method numerical simulation of a biaxial loading test
under constant lateral pressure. These simulations have pointed out the fundamental role
played by the elastic energy stored within the material before the onset of bifurcation. In
particular, the transition of the system from a homogenous regime toward a heterogeneous
one is strongly governed by the maximal capacity of the system to store elastic energy and
which reaches a non-zero plateau after the peak stress. The system is no longer able to store
more elastic energy during the post-bifurcation critical state regime. Simultaneously, the
entropy production reaches an extremum upon entering the critical state regime. After peak
stress, the entropy production decreases to a constant value as the critical state is sustained.
Thus, the bifurcation of the system from a homogeneous state toward a heterogeneous one
into a shear band corresponds to a structuring of the material that ultimately minimizes the
entropy production. In essence, for the system to dissipate energy in an optimal manner, it
must be made through a shear band where plastic dissipation is the most intense. This result
suggests a new interpretation of the intriguing shear banding phenomenon as the emergence
of spontaneous dissipative structures by exhaustion of available energy.
ACKNOWLEDGMENTS

We gratefully acknowledge the CNRS International Research Network GeoMech for having
offered the opportunity to make this project possible through long-standing collaboration of
all the authors (http://gdrmege.univ-lr.fr/).

APPENDIX 1

Derivation of equation (12)

The material system at hand is assumed to be at equilibrium at time t. An infinitesimal change
in the loading parameters is applied between time t, and time t + Δt. The evolution of the
system should satisfy the kinetic energy theorem, that expresses in Lagrangian form as:

$$
\dot{E}_c(t + \Delta t) = \int_{\Gamma_o} f_i(t + \Delta t) \dot{u}_i(t + \Delta t) ds_o - \int_{V_o} \Pi_{ij}(t + \Delta t) \frac{\partial \dot{u}_i(t + \Delta t)}{\partial X_i} dv_o \quad (A1)
$$

At time *t*, the first-orders time derivatives in Eq. (A1) can be assumed to equal zero.

Expanding $\Pi_{ij}(t + \Delta t)$ and $f_i(t + \Delta t)$ to the first time-order gives:

$$
\Pi_{ij}(t + \Delta t) = \Pi_{ij}(t) + \frac{\Delta t}{2} (\dot{\Pi}_{ij}(t + \Delta t) + \dot{\Pi}_{ij}(t)) = \Pi_{ij}(t) + \frac{\Delta t}{2} \dot{\Pi}_{ij}(t + \Delta t) \quad (A2)
$$

$$
f_i(t + \Delta t) = f_i(t) + \frac{\Delta t}{2} (\dot{f}_i(t + \Delta t) + \dot{f}_i(t)) = f_i(t) + \frac{\Delta t}{2} \dot{f}_i(t + \Delta t) \quad (A3)
$$

Combining Eqs. (A1) and (A2) yields:

$$
\begin{equation}
\begin{aligned}
\dot{E}_c(t + \Delta t) = \frac{\Delta t}{2} & \left( \int_{\Gamma_o} \dot{f}_i(t + \Delta t) \dot{u}_i(t + \Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t + \Delta t) \frac{\partial \dot{u}_i(t + \Delta t)}{\partial X_i} dv_o \right) + \dots \right. \\
& \left. \qquad + \int_{\Gamma_o} f_i(t) \dot{u}_i(t + \Delta t) ds_o - \int_{V_o} \Pi_{ij}(t) \frac{\partial \dot{u}_i(t + \Delta t)}{\partial X_i} dv_o \right)
\end{aligned}
\tag{A4}
\end{equation}
$$

By virtue of Gauss theorem and with help of Eq. (A3), Eq. (A4) reads:

$$
\begin{equation}
\begin{aligned}
\dot{E}_c(t + \Delta t) = \frac{\Delta t}{2} & \left( \int_{\Gamma_o} \dot{f}_i(t + \Delta t) \dot{u}_i(t + \Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t + \Delta t) \frac{\partial \dot{u}_i(t + \Delta t)}{\partial X_i} dv_o \right) + \dots \right. \\
& \left. + \int_{\Gamma_o} (f_i(t) - \Pi_{ij}(t) n_j) \dot{u}_i(t + \Delta t) ds_o + \Delta t \int_{V_o} \frac{\partial \Pi_{ij}(t)}{\partial X_i} \dot{u}_i(t + \Delta t) dv_o \right)
\end{aligned}
\tag{A5}
\end{equation}
$$
As the system is at equilibrium at time $t$, $\frac{\partial \Pi_{ij}(t)}{\partial X_i} = 0$ and $f_i(t) - \Pi_{ij}(t) n_j = 0$

Furthermore,

$$E_c(t + \Delta t) = E_c(t) + \frac{\Delta t}{2} (\dot{E}_c(t + \Delta t) + \dot{E}_e(t)) = E_c(t) + \frac{\Delta t}{2} \dot{E}_e(t + \Delta t) \quad (A6)$$

which finally gives:

$$E_c(t + \Delta t) = E_c(t) + \frac{(\Delta t)^2}{4} \left( \int_{\Gamma_o} \dot{f}_i(t + \Delta t) \dot{u}_i(t + \Delta t) ds_o - \int_{V_o} \dot{\Pi}_{ij}(t + \Delta t) \frac{\partial \dot{u}_i(t + \Delta t)}{\partial X_i} dv_o \right) \quad (A7)$$

# APPENDIX 2

The discrete Element Method (DEM) in brief

Down to the microscale, granular materials are modeled as poly-disperse assemblies of non-deformable spheres interacting through contact laws. In the present paper, we use the simple elasto-frictional contact law proposed by Cundall (1979) and illustrated in Fig. A1.

Two spherical particles are said to be in contact if they overlap. Then, the repulsive normal force reads $F_n = k_n u_n$ where $u_n$ is the overlapping distance and $k_n$ is the contact normal stiffness. In addition to the normal force, a tangential force $F_t$ develops depending on the grains' relative displacements. This tangential force is defined in an incremental form by computing the tangential component of the incremental relative displacement $\Delta u_t$ at the contact point between two contacting grains. $\Delta u_t$ depends on the relative translation and rotation speeds of the two grains in contact. $F_t$ is then updated according to the tangential stiffness $k_t = \alpha k_n$ expressed as a fraction $\alpha$ of its normal counterpart. The tangential elasticity is capped by a maximum value depending on the contact friction angle $\varphi$. Following Coulomb's friction law, the ratio $F_t / F_n$ must remain smaller than $\tan \varphi$.

After computing all inter-particle contact forces, the induced particles displacements are integrated based on Newton's second law of motion using an explicit integration scheme over a time step $dt$.
<figure><img src="image_10.png" /><figcaption>Fig. A1. Elasto-frictional contact law used in DEM simulations.</figcaption></figure>

In the simulations presented in this paper, all the boundary walls are considered rigid, although flexible membranes or periodic boundary conditions could be used as well. The choice of the type of boundary conditions influences both the occurrence and the pattern of the shear band. Hence, the results presented here should be understood in the context of the boundary condition chosen. However, it is believed that these results remain generic enough to illustrate the matter of bifurcation as addressed in this work.

From DEM simulations, it is straightforward to compute the elastic energy of the system as the sum of elastic energy stored at each contact ‘c’. The elastic energy thus reads:

$$E^{el} = \sum_c \left( \frac{F_n^2}{2k_n} + \frac{F_t^2}{2k_t} \right) \qquad (A8)$$

Computing the plastic dissipation directly at the contact level is trickier as it requires decomposing the incremental tangential displacement into elastic and plastic parts at every time step. For the sake of computation time, we wrote the energy balance directly at the sample scale. From the knowledge of the forces and velocities acting on the sample boundary, and assuming mechanical equilibrium between the inner granular specimen and the boundary loading, the external power $\dot{W}^{ext}$ provided to the system reads:

$$\dot{W}^{ext} = \int_{\partial V} f_i \dot{u}_i ds = \dot{E}^{el} + \dot{E}^{pl} \qquad (A9)$$

where $\dot{E}^{pl}$ is the plastic power. The advantage of using such a macroscopic approach for computing $\dot{E}^{pl}$ is that we can use a time step $\Delta t$ much larger than the DEM time step $dt$, which is very convenient from a computation point of view. We have therefore:
$$ \dot{E}^{pl}(t+\Delta t) = \dot{W}^{ext}(t+\Delta t) - \frac{E^{el}(t+\Delta t) - E^{el}(t)}{\Delta t} \quad (A10) $$

# APPENDIX 3

## Piecewise fitting method using polynomials

Based on the energy computation described in Appendix 2, the evolution of the plastic energy $E^{pl}$ can be obtained directly along the biaxial loading path. However, due to the numerical fluctuations, the rate of energy dissipation $\theta P^{int} = \frac{dE^{pl}}{dt}$ displays a very noisy evolution, which makes it difficult to be analyzed. Thus, the raw evolution of plastic energy $E^{pl}$ with respect to axial strain $\epsilon_1$ is fitted piecewise with a sixth-order polynomial function over every 1.0% of axial strain range. By differentiating the resulting fitted function $f(\epsilon_1)$, the rate of energy dissipation is directly obtained as follows:

$$ \theta P^{int} = \frac{df}{d\epsilon_1} \dot{\epsilon}_1 \quad (A11) $$
REFERENCE

Aziz-Alaoui, M., and Bertelle, C. (Eds.) (2009): From system complexity to emergent properties. Springer Science & Business Media.

Bažant, Z.P (1988): Stable states and paths of structures with plasticity or damage. Journal of Engineering Mechanics, ASCE, Vol. 114(12), pp. 2013-2033.

Bažant, Z.P, and Feng-bao, Lin (1989): Stability against localization of softening into ellipsoids and bands: Parameter study. International Journal of Solids and Structures, Vol. 25(12), pp. 1483-1498.

Bénard, H. (1901): Les tourbillons cellulaires dans une nappe liquide transportant de la chaleur par convection en régime permanent. Annales de chimie et de physique, Vol. 7(23), pp. 62-144.

Bigoni, D., and Hueckel, T. (1991): Uniqueness and localization, I. Associative and non-associative elastoplasticity. International Journal of Solids and Structures. Vol. 28(2), pp. 197-213.

Chandrasekhar, S. (1961): Hydrodynamic and hydromagnetic stability. Oxford, Clarendon Press.

Collins, I. F., and Houlsby, G. T. (1997): Application of thermomechanical principles to the modelling of geotechnical materials. Proceedings - Royal Society of London, A 453, pp. 1975-2001.

Collins, I. F., and Muhunthan, B. (2003): On the relationship between stress-dilatancy, anisotropy, and plastic dissipation for granular materials. Geotechnique, 53(7): 611-618.

Collins, I.F. (2005): The concept of stored plastic work or frozen elastic energy in soil mechanics. Géotechnique, 55(5), pp. 373-382.

Cundall, P.A, and Strack, O.D.L. (1979): A discrete numerical model for granular assemblies. Geotechnique, Vol. 29(1), pp. 47-65.

Damper, R.I. (2000): Emergence and levels of abstraction. Editorial for the Special Issue on 'Emergent Properties of Complex Systems', International Journal of Systems Science, Vol. 31(7), pp. 811-818.

Darve, F., and Vardoulakis, I. (2005): Instabilities and degradations in geomaterials. Darve and Vardoulakis Eds., Springer publ.

Deng, N., Wautier, A., Thiery, Y., Yin, Z.Y., Hicher, P.Y., and Nicot, F. (2021): On the attraction power of critical state in granular materials. J. Mech. and Physics of Solids, Vol. 149, 104300.
Deng, N., Wautier, A., Tordesillas, A., Thiery, Y., Yin, Z.Y., Hicher, P.Y., and Nicot, F. (2022): Lifespan dynamics of cluster conformations in stationary regimes in granular materials. Physical Review/E, in press.

Desrues, J., and Andò, E. (2015): Strain localisation in granular media. Comptes Rendus. Physique, Académie des sciences (Paris), Vol. 16 (1), pp.26-36.

Desrues, J., and Chambon, R. (2002): Shear band analysis and shear moduli calibration. International Journal of Solids and Structures. Vol. 39, pp. 3757-3776.

Dewar, R. C. (2005). Maximum entropy production and the fluctuation theorem. Journal of Physics A: Mathematical and General, 38(21): L371-L381.

Evans, D., and Searle, D. (2002): The fluctuation theorem. Advances in Physics, 51(7): 1529-1585.

Glansdorff, P., and Prigogine, I. (1954): Sur les propriétés différentielles de la production d'entropie. Physica, Vol. 20(7-12), pp. 773-780.

Glansdorff, P., and Prigogine, I. (1963): Generalised entropy production and hydrodynamic stability. Physics Letters, Vol. 7(4), pp. 243-244.

Glansdorff, P., and Prigogine, I. (1964): On a general evolution criterion in macroscopic physics. Physica, Vol. 30(2), pp. 351-374.

Glansdorff, P., and Prigogine, I. (1971): Thermodynamic theory of structure, stability and fluctuations. Wiley-Interscience, N.Y.

Hill, R. (1958): A general theory of uniqueness and stability in elastic-plastic solids. Journal of the Mechanics and Physics of Solids, Vol. 6 (3), pp. 236-249.

Hill, R. (1962): Acceleration waves in solids. Journal of the Mechanics and Physics of Solids, Vol. 10(1), pp. 1-16.

Hudson, J.L., and Mankin, J.C. (1981). Chaos in the Belousov-Zhabotinskii reaction. J. Chem. Phys., Vol. 74 (11), pp. 6171-6177.

Lefever., R. (1978): Stabilité globale et mécanismes de fluctuations des structures dissipatives. Journal de Physique, 39 (C5), pp.C5-83-C5-92.

Liu, J., Nicot, F., and Zhou, W. (2018): Sustainability of internal structures during shear band forming in 2D granular materials. Powder Technology, Vol. 338, pp. 458-470.

Liu, J., Wautier, A., Bonelli, S., Nicot, F., and Darve, F. (2020): Macroscopic softening in granular materials from a mesoscale perspective. Int. J. Solids Structures, 93-194, pp. 222-238. DOI 10.1016/j.ijsolstr.2020.02.022.
Liu, J., Wautier, A., Nicot, F., Darve, F., and Zhou, W. (2022): How meso shear chains bridge multiscale shear behaviors in granular materials: a preliminary study. Int. J. Solids Structures, in Press.

Mandel, J. (1966) : Conditions de stabilite et Postulat de Drucker. In: Kravtechenko, Sirieys (Eds.), Rheology and Soil Mechanics. Springer-Verlag Publisher, pp. 58–68.

Murnaghan, F. D. (1944). The Compressibility of Media under Extreme Pressures. Proc Natl Acad Sci USA, 30(9): 244-247.

Nicolis, G. (1970): Thermodynamic theory of stability, structure and fluctuations. Pure and Applied Chemistry, Vol. 22, pp. 379 – 392.

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

Nicot, F., and Darve, F. (2006): On the elastic and plastic strain decomposition in granular materials. Granular Matter, Vol. 8(3-4), pp. 221-237.

Nicot, F., and Darve, F. (2007a): Basic features of plastic strains: from micro-mechanics to incrementally nonlinear models. Int. Journal of Plasticity, Vol. 23, pp. 1555-1588.

Nicot, F., and Darve, F. (2007b): Micro-mechanical bases of some salient constitutive features of granular materials. Int. J. of Solids and Structures, Vol. 44, pp. 7420–7443.

Nicot, F., and Darve, F. (2011): Diffuse and localized failure modes: two competing mechanisms. International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 35, Issue 5, pp. 586–601.

Nicot, F., Sibille, L., and Darve, F. (2012): Failure in rate-independent granular materials as a bifurcation toward a dynamic regime. Int. Journal of Plasticity, Vol. 29, pp. 136–154.

Parisì, G., and Sourlas, N. (2002): Scale invariance in disordered systems: the example of the random-field ising model. Phys. Rev. Lett., Vol. 89, 257204.

Parisì, G., Procaccia, I., Rainone, C., Singh, M. (2017): Shear bands as manifestation of a criticality in yielding amorphous solids. Proc. Natl. Acad. Sci. U. S. A., Vol. 114, pp. 5577–5582. https://doi.org/10.1073/pnas.1700075114.

Petryk, H. (1993): Theory of bifurcation and instability in time-independent plasticity. In Bifurcation and Stability of Dissipative Systems, Nguyen QS (ed.). CISM Courses and Lecturers, vol. 327. Springer: Berlin, pp. 95–152.

Prigogine I., and Stengers, I. (1979): The New Alliance. Gallimard Ed., Paris.

Prigogine, I. (1945). Modération et transformations irréversibles des systèmes ouverts. Bulletin de la Classe des Sciences., Académie Royale de Belgique. Vol. 31, pp. 600–606.
Prigogine, I. (1947): *Étude thermodynamique des phénomènes irréversibles*. Paris et Liège, Dunod et Desoer Eds.

Prigogine, I., and Lefever, R. (1968): Symmetry Breaking Instabilities in Dissipative Systems. The Journal of Chemical Physics, vol. 48, no 4, pp. 1695–1700.

Prigogine, I., and Lefever, R. (1975): Stability and Self-Organization in Open Systems. Advances in Chemical Physics, Vol. 29, pp. 1–28.

Prigogine, I., Lefever, R., et al. (1969): Symmetry Breaking Instabilities in Biological Systems. Nature, Vol. 223(5209), pp. 913–916.

Rayleigh, J. W. (1916): On convection currents in a horizontal layer of fluid, when the higher temperature is on the underside. The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, Sixth series, Vol. 32(192), pp. 529-546.

Regenauer-Lieb, K., Hu, M., Schrank, C., Chen, X. et al. (2021): Cross-diffusion waves resulting from multiscale, multi-physics instabilities: theory. Solid Earth, 12(4): 869-883.

Rice, J.R. (1975): Continuum mechanics and thermodynamics of plasticity in relation to microscale deformation mechanisms. Constitutive equations in plasticity, A.S. Argon Ed., MIT Press, Cambridge, pp. 23-79.

Rudnicki, J.W., and Rice, J. (1975): Conditions for the localization of deformation in pressure sensitive dilatant materials. International Journal of Solids and Structures. Vol. 23, pp. 371–394.

Sun, Q., Jin, F., Wang, G., Song, S., and Zhang, G. (2015): On granular elasticity. Sci. Rep. 5, 9652; DOI:10.1038/srep09652.

Teixeira-Dias, J.C. (2017): Molecular physical chemistry: a computer-based approach using Mathematica and Gaussian. Springer International Publishing, pp. 1-82.

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

Tordesillas, A., and Muthuswamy, M. (2009): On the modeling of confined buckling of force chains. Journal of the Mechanics and Physics of Solids, Volume 57(4), pp. 706-727.

Turing, A.M. (1952): The Chemical Basis of Morphogenesis », Philosophical Transactions of the Royal Society of London B: Biological Sciences, Vol. 237(641), pp. 37–72.

Veveakis, E., and Regenauer-Lieb, K. (2015): Review of extremum postulates. Current Opinion in Chemical Engineering, 7(C), pp. 40-46.

Walker, D.M., and Tordesillas, A. (2010): Topological evolution in dense granular materials: a complex networks perspective. Int. J. of Solids and Structures, Vol. 47, pp. 624-639.
Wan, R., Nicot, F., and Darve, F. (2017): Failure in geomaterials, a contemporary treatise.
Wiley.

Wautier, A., Bonelli, S., and Nicot, F. (2018): Flow impact on granular force chains and
induced instability. Physical Review E, Vol. 98(4): 042909.

Ziegler, H. (1983): An Introduction to Thermomechanics. Amsterdam, North Holland.

Ziegler, H., and Wehrli, C. (1987). On a principle of maximal rate of entropy production. J.
Non-Equilib. Thermodyn., Vol. 12 (3), pp. 229–243.