# A Cosserat Breakage Mechanics model for brittle granular media

Nicholas Anton Collins-Craft, Ioannis Stefanou, Jean Sulem, Itai Einav

► To cite this version:

Nicholas Anton Collins-Craft, Ioannis Stefanou, Jean Sulem, Itai Einav. A Cosserat Breakage Mechanics model for brittle granular media. Journal of the Mechanics and Physics of Solids, 2020, 141, pp.103975. 10.1016/j.jmps.2020.103975. hal-03120686

HAL Id: hal-03120686

https://hal.science/hal-03120686

Submitted on 25 Jan 2021

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.
# A Cosserat Breakage Mechanics model for brittle granular media

Nicholas Anton Collins-Craft¹,², Ioannis Stefanou²,³, Jean Sulem², and Itai Einav¹

¹Particles and Grains Laboratory, The School of Civil Engineering, The University of Sydney,
Sydney, Australia

²Laboratoire Navier, École des Ponts ParisTech, Université Gustave Eiffel, CNRS,
Champs-sur-Marne, France

³GeM (Institute de Recherche en Genie Civile et Mécanique), École Centrale de Nantes, Université
de Nantes, CNRS, Nantes, France

## Abstract

This paper introduces a Breakage Mechanics model formulated in the framework of Cosserat continuum theory. This formulation allows us to account for the interaction of grain size evolution with shear band formation in a physically motivated and thermodynamically consistent way. We provide an upscaling procedure that introduces material parameters accounting for the contribution of the Cosserat rotations, as well as the micro-inertia of the entire grain size distribution. A particular highlight of the new formulation is the capacity to describe an evolving internal Cosserat length that takes account of the complete grain size distribution and its evolution. In addition, we specify a particular model that requires no additional calibration relative to the same model in the classical continuum. We apply this model to the study of a layer sheared under constant volume which simulates the fast undrained deformation observed in seismogenic faults. Through linear stability analysis we use the model to examine the effect of grain size polydispersity on the thickness of shear bands. By implementing the model using the finite element method we provide an explanation for the geological formation of double cataclastic shear bands.

**Keywords:** Granular materials, Breakage, Upscaling, Shear bands, Cosserat continuum

## 1 Introduction

The accurate prediction of the failure modes of geomaterials, such as rocks and brittle granular soils, is of wide interest in both the fields of geology and geotechnical engineering. In particular, failure due to various modes of localisation, including compaction, dilation and shear bands is of particular importance, as it has been claimed that these are the failure mechanisms associated with landslides (e.g Veveakis, Vardoulakis and Toro, 2007) and seismic slip (e.g Rice, 2006). Localised failure is typically described with the Thomas-Hill-Mandel shear band model, and can be studied using bifurcation analysis to describe the material behaviour both before and after the localisation occurs (Rudnicki and Rice, 1975; Vardoulakis and Sulem, 1995). Particularly in a geological context, multi-physical couplings have a strong influence on the formation of localisation zones with temperature and pore fluid pressure interactions being well studied (e.g. Sulem, 2010; Sulem, Stefanou and Veveakis, 2011) and chemical changes such as decomposition and dehydration (Brantut and Sulem, 2012; Brantut, Sulem and Schubnel, 2011; Stefanou and Sulem, 2013, 2014; Sulem and Famin, 2009; Veveakis, Alevizos and Vardoulakis, 2010; Veveakis, Stefanou and Sulem, 2013, among others) have also been shown to contribute towards forming but also potentially limiting the size of localisation zones. Recently, fully coupled systems that integrate the effects of temperature, pressurisation and chemical reactions have been explored (Platt, Rudnicki and Rice, 2014; Rattez, Stefanou and Sulem, 2018; Rattez et al., 2018a,b; Rice, Rudnicki and Platt, 2014; Sulem and Stefanou, 2016).

However, notably for the present study, the evolution of the size of the microstructure of the considered material also plays an important role in shear band development and apparent softening behaviour. The latter is connected to earthquake stability and control (Stefanou, 2019). Fault cores reveal that there are distinct zones within a more extensive damaged zone. Extremely thin zones of localisation, known as principal slipping zones may be observed, and it is suggested that these zones accommodate the majority of slip in a fault (Chester and Chester, 1998; Sibson, 2003; Sulem, 2007; Wibberley and Shimamoto, 2002). Within these zones, we may observe that the particle size distribution is shifted considerably towards finer sizes than those in the surrounding fault core, which implies a further evolution of the cataclastic grain crushing mechanism from damaged rock to granular material.
The evolution of the mechanical properties in fault cores is tightly coupled to the other physical phenomena, as changes in grain size will correlate to changed surface area, having an effect on chemical reactions (Buscarnera and Das, 2016; Stefanou and Sulem, 2014; Viswanath and Das, 2019; Zhang and Buscarnera, 2018), will feature heat energy being released due to friction and grain breakage, and are associated with changes in porosity, which will affect the pore water pressure (Rattez, 2018; Rattez, Stefanou and Sulem, 2018; Rattez et al., 2018a,b). As such, being accurately able to characterise this evolution in grain size is a pressing need in order to be able to develop a comprehensive model of fault mechanics.

A successful continuum mechanics theory for the prediction of grain size evolution is the Breakage Mechanics theory, a thermodynamically-consistent family of models that introduce an internal variable B, capturing the change of the grain size distribution (Einav, 2007a,b). By making the claim that under indefinite loading, the grain size distribution tends towards a definable ultimate distribution, B can be defined so that it tracks the progress of the grain size distribution from an initial state towards the assumed ultimate state. A power-law distribution, which is well supported experimentally and in field observations (Sammis, King and Biegel, 1987; Turcotte, 1986) is typically chosen. Using the thermodynamic principles of constitutive modelling, Breakage Mechanics models can be defined which range from simple models drawing analogies with fracture mechanics (Einav, 2007c,d) to models that describe granular mixtures (Einav and Valdes, 2008), large flow problems (Zhang, 2012; Zhang, Einav and Nguyen, 2012), the energetics of cataclasis (Einav and Nguyen, 2008; Nguyen and Einav, 2009), the formation of compaction bands (Das, 2013; Das, Nguyen and Einav, 2011, 2013), partial saturation effects (Buscarnera and Einav, 2012), inter-particle cementation effects (Das et al., 2014; Tengattini, 2015; Tengattini, Das and Einav, 2014; Tengattini et al., 2014), coupling with fracture mechanics (Zhang and Buscarnera, 2014, 2015; Zhang, Buscarnera and Einav, 2016), shock loading conditions (Herbold, Homel and Rubin, 2020), and the non-unique dependence of critical state on porosity (Tengattini, Das and Einav, 2016). Of particular note is a formulation of Breakage mechanics in a non-local continuum (Nguyen and Einav, 2010).

When we consider how exactly to model these phenomena, it can be shown that without some form of regularisation, the predicted width of any localisation band will be zero (Sulem, Stefanou and Veveakis, 2011). There are a number of different ways that this can be resolved, but they all involve equipping the continuum with a length scale. The most commonly used technique is viscous regularisation, which adds a rate-dependent effect. However, at slow rates of loading where inertial effects are practically negligible, these techniques result in predicted localisations widths that depend on the size of the perturbations that may arise due to imperfections or boundary conditions. Non-local models, where the behaviour at a material point depends on the behaviour of other material points within a certain distance of it has already been integrated with Breakage Mechanics (Nguyen and Einav, 2010). However, in this implementation, the non-local length scale is dependent on the choice of the weighting function (e.g. Gaussian, Lucy etc), and does not evolve with the underlying physical structure of the grains. While second gradient theories have also demonstrated success in modelling localisations in geomaterials (Zervos, Papanastasiou and Vardoulakis, 2001), and there has been some recent advances to include an evolving length scale in second gradient (Scherer et al., 2019) and non-local damage mechanics (Pijaudier-Cabot, Haidar and Dubé, 2004; Rastiello et al., 2018), or to account for system evolution by changing the non-local weighting or interaction functions (Nguyen, 2011; Poh and Sun, 2017; Wang and Poh, 2018), the length scale has not yet been attached to any particular micro-mechanical quantity that could give it physical meaning. The mathematical structures of non-local models do not imply any particular expression for the length scale evolution as a consequence of changes in the physical media, and this has yet to be attempted in the context of Breakage Mechanics. However, as we will show in §3, the Cosserat continuum does suggest a particular expression that captures the evolving length scale as a consequence of its inherent structure when combined with Breakage Mechanics. This continuum, where additional rotational degrees of freedom are given to the material point (Cosserat and Cosserat, 1909; Germain, 1973), has also previously demonstrated success when modelling physical phenomena where the micro-structure plays an important role (Mühlhaus and Vardoulakis, 1987; Papanastasiou and Vardoulakis, 1989, 1992; Stefanou, Sulem and Vardoulakis, 2008, 2010; Vardoulakis and Papanastasiou, 1988; Vardoulakis and Sulem, 1995). We also develop our model with an eye to integrating it with previously developed multi-physics models that have been regularised in the Cosserat framework (Rattez, Stefanou and Sulem, 2018; Rattez et al., 2018a,b). These Cosserat models associate the internal lengths of the Cosserat continuum with that of the mean grain size. Some recent works (Bauer, 2016; Bauer, Li and Huang, 2015; Bauer, Li and Khosravi, 2017) have included an empirically based evolution equation for the mean grain size using a hypoplastic Cosserat formulation that is indifferent to thermodynamics, and thus cannot be used to correctly capture thermal processes in fault physics.

Hence, it is clear that combining the Cosserat continuum with Breakage mechanics in a consistent thermodynamic way can potentially offer substantial benefits to our capacity to model physical phenomena by accounting for both the entire grain size distribution, and its evolution. In particular, being able to predict the formation of localisation bands in dry cohesion-less granular media can form the foundation of richer constitutive models that can eventually
be applied to problems arising from fault mechanics.

# 2 Definitions

## 2.1 Strain and curvature rates

In a three-dimensional Cosserat continuum, the material points have six degrees of freedom, those of translation that exist also in the classical continuum, and three degrees of rotational freedom.

The infinitesimal (non-symmetric) strain rate tensor is defined by

$$ \dot{\gamma}_{ij} = \dot{u}_{i,j} + \epsilon_{ijk} \dot{\omega}_k^c, \quad (1) $$

and the infinitesimal (non-symmetric) curvature rate tensor is defined by

$$ \dot{\kappa}_{ij} = \dot{\omega}_{i,j}^c, \quad (2) $$

where $\dot{u}_i$ and $\dot{\omega}_i^c$ are the translation and rotation rates along and about the $x_i$ axes, $\epsilon_{ijk}$ is the Levi-Civita symbol, while we follow the Einstein summation convention.

The strain rate tensor can be split into symmetric and anti-symmetric parts, denoted $\dot{\gamma}_{(ij)}$ and $\dot{\gamma}_{[ij]}$ respectively (similarly with the curvature rate tensor). The symmetric part of the strain rate tensor is identical to the infinitesimal strain rate tensor in the classical continuum, $\dot{\epsilon}_{ij}$.

Both strain and curvature rate tensors can also be split into trace and deviatoric parts:

$$ \dot{\gamma}_{ij} = \frac{1}{3} \dot{\epsilon}_{kk} \delta_{ij} + \dot{e}_{ij}, \quad (3) $$

$$ \dot{\kappa}_{ij} = \frac{1}{3} \dot{\kappa}_{kk} \delta_{ij} + \dot{k}_{ij}, \quad (4) $$

where $\delta_{ij}$ is the Kronecker delta. The deviatoric parts can be further decomposed into symmetric and anti-symmetric components if so desired.

## 2.2 Plasticity framework

Later in this paper, we develop a constitutive model within the context of rate-independent processes. The strain and curvature rates are decomposed into elastic (recoverable) and plastic (non-recoverable) parts:

$$ \dot{\gamma}_{ij} = \dot{\gamma}_{ij}^e + \dot{\gamma}_{ij}^p, \quad (5) $$

$$ \dot{\kappa}_{ij} = \dot{\kappa}_{ij}^e + \dot{\kappa}_{ij}^p. \quad (6) $$

The decompositions of the elastic and plastic parts into their respective volumetric and deviatoric, symmetric and anti-symmetric parts follow as for the total rate tensors above.

We hold to the convention that constitutive models must not rely on arbitrary reference-dependent variables (Rubi, 2001). As such, the elastic strain and curvature tensors, $\gamma_{ij}^e$ and $\kappa_{ij}^e$, are used as state variables in the model, but not the plastic and total strain and curvature tensors.

## 2.3 Stresses and couple-stresses

The stress tensor $\tau_{ij}$ in a Cosserat continuum is, in general, non-symmetric like the strain rate tensor. The tensor can be decomposed into symmetric and anti-symmetric parts, and the symmetric part is identified with the Cauchy stress tensor of the classical continuum, $\sigma_{ij}$. The couple-stress tensor $\mu_{ij}$ can be decomposed similarly.

Both of these tensors can also be decomposed into volumetric and deviatoric parts:

$$ \tau_{ij} = \frac{1}{3} \tau_{kk} \delta_{ij} + s_{ij}, \quad (7) $$

$$ \mu_{ij} = \frac{1}{3} \mu_{kk} \delta_{ij} + m_{ij}, \quad (8) $$

where $\frac{1}{3}\tau_{kk}$ is equivalent to the mean stress $p$. As with the strain, the deviatoric stresses and couple stresses can be further split into symmetric and anti-symmetric parts.
## 2.4 Breakage state variable

The breakage state variable B is defined by the ratio of the area between the current and initial cumulative grain size distributions (GSDs) to the area between the ultimate and initial cumulative grain size distributions, as shown in Figure 1.

<figure><img src="image_1.png" /><figcaption>Figure 1: Definition of the breakage variable B as the ratio of two grading areas Buscarnera and Einav, 2012; Einav, 2007a. $F_0(x)$, $F_u(x)$ and $F(x)$ are the initial, ultimate and current cumulative GSDs respectively. $d_{min}$ and $d_{Max}$ are the minimum and maximum grain sizes present in the GSDs.</figcaption></figure>

The maximum grain size $d_{Max}$ is constant between the initial and ultimate distributions, and $d_{min}$ is the smallest grain size possible to obtain due to comminution, typically taken to be one micron (Kendall, 1978).

This definition of breakage is consistent with that given in Einav (2007a) and Buscarnera and Einav (2012). In particular, we make use of the universal initial distribution from Buscarnera and Einav (2012), using a Heaviside step function at $d_{Max}$ for the cumulative initial distribution. We express the probability density function for the current GSD as being a linear function of B, ranging between those corresponding to the initial and ultimate distributions:

$$p(x, B) = p_0(x)(1 - B) + p_u(x)B. \qquad (9)$$

Using the Heaviside step function for the cumulative initial distribution means that the probability density function for the initial distribution is given by

$$p_0(x) = \delta(x - d_{Max}), \qquad (10)$$

where $\delta$ represents the Dirac delta function. We note that the introduction of B also requires the introduction of its thermodynamic conjugate, the breakage energy $E_B$, per Einav (2007a,b)

# 3 Upscaling procedures

## 3.1 Elastic upscaling

The previous models of Breakage Mechanics that we discussed in §1 have been formulated in either the classical continuum, or a non-local continuum. To formulate Breakage Mechanics in a Cosserat continuum, a statistical homogenisation process (upscaling) is undertaken that takes account of the contribution of the stored energy of the elastic curvatures when upscaled from the granular medium to the continuum formulation.

This process starts by considering the Helmholtz free energy density per unit volume. The function is given by

$$\Psi = \Psi(\gamma_{ij}^e, \kappa_{ij}^e, B) = \int_{d_{min}}^{d_{Max}} \hat{\psi}(\gamma_{ij}^e, \kappa_{ij}^e, x)p(x, B) dx. \qquad (11)$$

That is to say, the total Helmholtz free energy density is the Helmholtz free energy density at a given grain size x, weighted by the grain size distribution $p(x, B)$, integrated over the entire distribution.

The energy density is assumed to be further decomposed into a summation of elastic strain contributions and elastic curvature contributions:

$$\hat{\psi}(\gamma_{ij}^e, \kappa_{ij}^e, x) = \hat{\psi}_{\gamma}(\gamma_{ij}^e, x) + \hat{\psi}_{\kappa}(\kappa_{ij}^e, x). \qquad (12)$$
We make this assumption for the sake of simplicity, and also because a multiplicative form would cause homogeneous elastic strains to induce couple effects, and we are not aware of any evidence for this occurring. This form also has the benefit of having been shown to result in unique solutions (Eringen, 1967).

Now we use the assumption of fractional independence from Einav (2007a) to make a separation of variables. The Helmholtz free energy contribution for a grain size $x$ due to the strains is

$$ \hat{\psi}_{\gamma}(\gamma_{ij}^{e}, x) = \left(\frac{x}{x_r}\right)^2 \psi_{r\gamma}(\gamma_{ij}^{e}), \qquad (13) $$

where the subscript $r$ represents a reference grain size fraction.

We can apply the same assumption of fractional independence to the contribution of the curvatures. However, in order for the energy density to have the correct dimensionality to contribute correctly to the Helmholtz free energy the elastic curvature terms must be further multiplied by some length scale $\ell^2$. Micro-mechanically, smaller particles have fewer particle contacts and so participate in fewer force chains. They are also typically less angular, meaning there is less particle interlocking (Andò et al., 2012). Hence, we expect that the rotation rate of the grains and thus ultimately the elastic curvatures will scale with the grain size $x$. Hence, we set $\ell^2 = x^2$ and state that the Helmholtz free energy contribution for a grain size $x$ due to the elastic curvatures is

$$ \hat{\psi}_{\kappa}(\kappa_{ij}^{e}, x) = \left(\frac{x}{x_r}\right)^2 \hat{\psi}_{r\kappa}(\kappa_{ij}^{e}) = \left(\frac{x}{x_r}\right)^4 x_r^2 \psi_{r\kappa}(\kappa_{ij}^{e}). \qquad (14) $$

Making successive back substitutions of the relevant expressions, the Helmholtz free energy expression becomes

$$ \Psi = \int_{d_{min}}^{d_{Max}} \left( \left(\frac{x}{x_r}\right)^2 \hat{\psi}_{r\gamma}(\gamma_{ij}^e) + \left(\frac{x}{x_r}\right)^4 x_r^2 \psi_{r\kappa}(\kappa_{ij}^e) \right) p(x, B) dx. \qquad (15) $$

We can split (15) into two integrals:

$$ \Psi = \Psi_{\gamma} + \Psi_{\kappa} = \frac{\psi_{r\gamma}(\gamma_{ij}^{e})}{x_r^2} \int_{d_{min}}^{d_{Max}} x^2 p(x, B) dx + x_r^2 \frac{\psi_{r\kappa}(\kappa_{ij}^{e})}{x_r^4} \int_{d_{min}}^{d_{Max}} x^4 p(x, B) dx. \qquad (16) $$

Now, we note that the $n$th statistical moment of the distribution (about the origin) is defined by

$$ J^{(n)}(B) = \int_{d_{min}}^{d_{Max}} x^n p(x, B) dx. \qquad (17) $$

Simple inspection of (16) shows that the two integrals give the second and fourth moments of the current grain size distribution. We also note that due to our choice in (10), evaluating (17) for $p_0(x, B)$ gives the $n$th moment of that distribution:

$$ J_{0}^{(n)} = \int_{d_{min}}^{d_{Max}} x^n p_0(x, B) dx = d_{Max}^{n}. \qquad (18) $$

As the maximum grain size is the only one present in all the GSDs, from the initial to the ultimate, it is the only physically meaningful choice for the reference grain size $x_r$. Therefore we use

$$ x_r = d_{Max}. \qquad (19) $$

We now turn our attention to the strain integral $\Psi_{\gamma}$. Substituting in (9) and evaluating the integral returns the moments of the initial and ultimate grain size distributions. Conveniently, because of the equality in (19), we may introduce the grading constant $\theta_{\gamma}$ by

$$ \theta_{\gamma} = 1 - \frac{J_{u}^{(2)}}{d_{Max}^{2}}, \qquad (20) $$

which allows us to re-arrange $\Psi_{\gamma}$ in (16) to obtain

$$ \Psi_{\gamma} = (1 - \theta_{\gamma} B) \psi_{r\gamma}(\gamma_{ij}^{e}), \qquad (21) $$

which is the classical Breakage Mechanics expression for the Helmholtz free energy per Einav (2007a), noting however that the strain tensor is not symmetric here.

This expression can be further separated into symmetric and anti-symmetric contributions:

$$ \Psi_{\gamma} = \Psi_{(\gamma)} + \Psi_{[\gamma]} = (1 - \theta_{\gamma} B)(\psi_{r(\gamma)}(\gamma_{ij}^{e}) + \psi_{r[\gamma]}(\gamma_{[ij]}^{e})). \qquad (22) $$
This process can be repeated for the curvature integral $\Psi_{\kappa}$ in (16), once again substituting in (9) and (19) and evaluating the integral to obtain the fourth moments of the initial and ultimate distributions. Using these relationships, we can introduce a new granulometry constant $\theta_{\kappa}$ defined by

$$ \theta_{\kappa} = 1 - \frac{J_u^{(4)}}{d_{Max}^4}, \qquad (23) $$

which allows us to rewrite $\Psi_{\kappa}$ in the form

$$ \Psi_{\kappa} = (1 - \theta_{\kappa} B) d_{Max}^2 \psi_{r\kappa}^*(\kappa_{ij}^e). \qquad (24) $$

This formulates the stored energy due to the elastic curvatures in a manner consistent with classical Breakage Mechanics.

Hence, the total Helmholtz free energy can be written

$$ \Psi = (1 - \theta_{\gamma} B)(\psi_{r(\gamma)}(\gamma_{(ij)}^e) + \psi_{r[\gamma]}(\gamma_{[ij]}^e)) + (1 - \theta_{\kappa} B) d_{Max}^2 \psi_{r\kappa}(\kappa_{ij}^e). \qquad (25) $$

We now compare this expression for the total elastic energy with that given in the extensive treatment of Cosserat continuum mechanics by Vardoulakis (2019) for isotropic linear Cosserat elasticity:

$$ \Psi = \frac{2\nu G}{1 - 2\nu} \epsilon_{ii}^{e} \epsilon_{jj}^{e} + G\epsilon_{ij}^{e} \epsilon_{ij}^{e} + \eta_{1} G\gamma_{[ij]}^{e} \gamma_{[ij]}^{e} + 6\eta_{1} G\gamma_{k}^{*} \gamma_{k}^{*} + \frac{1}{2} G\ell^{2}(\kappa_{(ij)}^{e} \kappa_{(ij)}^{e} + \eta_{2} \kappa_{(ii)}^{e} \kappa_{(jj)}^{e} + \eta_{3} \kappa_{[ij]}^{e} \kappa_{[ij]}^{e}), \qquad (26) $$

where $\gamma_k^*$ is an axial vector encompassing the relative rotation (the difference between the macro-rotations and the Cosserat rotations). The elastic length terms multiplying the terms associated with the elastic curvatures can be directly identified in the form given by Vardoulakis (2019). Identifying the comparable terms in our energetic expression (25) gives an expression for the elastic length when Breakage Mechanics is upscaled in the Cosserat continuum:

$$ \ell^e = d_{Max} \sqrt{1 - \theta_{\kappa} B}. \qquad (27) $$

We highlight here that the elastic length scale we have derived depends on the entire grain size distribution through the new grading constant $\theta_{\kappa}$, and evolves as the Breakage variable evolves. It should also be noted that as a consequence of the length scale being multiplied out of the expression for the energy due to curvature, the elastic parameters will have units solely of stress, and not couple-stress unlike other authors (for example Rattez, Stefanou and Sulem (2018)).

## 3.2 Inertial upscaling

In the Cosserat continuum the rotational momentum balance must be maintained in addition to the linear momentum balance, and hence the moment of rotational inertia, $I$, of the medium must be defined in order to fully specify the continuum. Previously, this has been taken as the moment of inertia associated only with the average particle size (Rattez, Stefanou and Sulem, 2018; Rattez et al., 2018a,b; Sulem, Stefanou and Veveakis, 2011). However, Breakage Mechanics gives us access to information about the entire grain size distribution, therefore, we take the average of the moments of inertia of the entire grain size distribution as representative. We model the evolution of this average value as the grain size distribution evolves due to grain crushing.

An upscaling on the three dimensional case is performed, as it necessitates the definition of a new material parameter, whereas the two dimensional case does not. The rotational moment of inertia of a solid sphere of diameter $x$ about its own axis is given by

$$ I = \frac{\pi}{60} \rho x^5, \qquad (28) $$

where $\rho$ is the density of the particle.

As particle breakage progresses, more (and smaller) particles appear inside the representative volume element (RVE). It is considered that each of these particles rotates around its own proper axis, but not around any other point. To find the average of the moments of inertia, the integration is performed over the continuous particle size distribution, using the substitution in (19):

$$ I(B) = \frac{\pi}{60} \rho \int_{d_{min}}^{d_{Max}} x^5 p(x, B) dx = \frac{\pi}{60} \rho d_{Max}^5 \int_{d_{min}}^{d_{Max}} \left( \frac{x}{d_{Max}} \right)^5 (p_0(x)(1 - B) + p_u(x)B) dx. \qquad (29) $$

Evaluating the integral and using the definitions of the fifth moments $J_0^{(5)}$ and $J_u^{(5)}$ results in

$$ I(B) = \frac{\pi}{60} \rho d_{Max}^5 \frac{(1 - B)J_0^{(5)} + BJ_u^{(5)}}{d_{Max}^5}. \qquad (30) $$
Once again exploiting that $J_0^{(5)} = d_{Max}^5$ and introducing another granulometry constant:

$$ \theta_I = 1 - \frac{J_u^{(5)}}{d_{Max}^5}, \qquad (31) $$

the expression for the rotational inertia can be re-arranged to

$$ I(B) = \frac{\pi}{60} (1 - \theta_I B) \rho d_{Max}^5, \qquad (32) $$

in three dimensions.

It should be noted that the rotational moment of inertia for a disc of diameter $x$ in two dimensions is

$$ I = \frac{\pi}{32} \rho x^4. \qquad (33) $$

If the above analysis were repeated in two dimensions, the definitions of $\theta_I$ and $\theta_\kappa$ would be identical. Hence, the rotational inertia is

$$ I(B) = \begin{cases} \frac{\pi}{32} (1 - \theta_\kappa B) \rho d_{Max}^4 & \text{in two dimensions,} \\ \frac{\pi}{60} (1 - \theta_I B) \rho d_{Max}^5 & \text{in three dimensions.} \end{cases} \qquad (34) $$

## 3.3 Limit values

The expressions for $\theta_\gamma$, $\theta_\kappa$ and $\theta_I$ can be written by evaluating the integrals in their explicit forms. Using the universal initial distribution and the ultimate fractal distribution consistent with Buscarnera and Einav (2012) and Einav (2007a), and the choice of reference length scale (19), the expressions are

$$ \theta_\gamma = 1 - \frac{3 - \alpha}{5 - \alpha} \left( \frac{1 - (d_{min}/d_{Max})^{5-\alpha}}{1 - (d_{min}/d_{Max})^{3-\alpha}} \right), \qquad (35) $$

$$ \theta_\kappa = 1 - \frac{3 - \alpha}{7 - \alpha} \left( \frac{1 - (d_{min}/d_{Max})^{7-\alpha}}{1 - (d_{min}/d_{Max})^{3-\alpha}} \right), \qquad (36) $$

$$ \theta_I = 1 - \frac{3 - \alpha}{8 - \alpha} \left( \frac{1 - (d_{min}/d_{Max})^{8-\alpha}}{1 - (d_{min}/d_{Max})^{3-\alpha}} \right). \qquad (37) $$

First, considering the limit $d_{min} \to 0$

$$ \theta_\gamma = \frac{2}{5 - \alpha}, \qquad (38) $$

$$ \theta_\kappa = \frac{4}{7 - \alpha}, \qquad (39) $$

$$ \theta_I = \frac{5}{8 - \alpha}. \qquad (40) $$

For typical values of $\alpha$ (2.5-2.7 (Sammis, King and Biegel, 1987)), $\theta_\gamma = 0.8-0.83$, $\theta_\kappa = 0.89-0.91$ and $\theta_I = 0.91-0.93$.

Next, considering the limit $\alpha \to 3$

$$ \theta_\gamma = 1 + \frac{1 - (d_{min}/d_{Max})^2}{2 \ln(d_{min}/d_{Max})}, \qquad (41) $$

$$ \theta_\kappa = 1 + \frac{1 - (d_{min}/d_{Max})^4}{4 \ln(d_{min}/d_{Max})}, \qquad (42) $$

$$ \theta_I = 1 + \frac{1 - (d_{min}/d_{Max})^5}{5 \ln(d_{min}/d_{Max})}. \qquad (43) $$

If a value such as $d_{min} = 0.00001d_{Max}$ is assigned, $\theta_\gamma = 0.957$, $\theta_\kappa = 0.978$ and $\theta_I = 0.98$.

From this analysis it can be concluded that $\theta_\gamma$, $\theta_\kappa$ and $\theta_I$ approach 1 only in the extreme limit as both $\alpha \to 3$ and $d_{min} \to 0$ at the same time.
<figure><img src="image_2.png" /><figcaption>Figure 2: (a) $\theta_{\gamma}$, (b) $\theta_{\kappa}$ and (c) $\theta_{I}$ as a function of $d_{min}/d_{Max}$ and $\alpha$. The insets are the respective plots holding $d_{min}/d_{Max} = 10^{-3}$ and varying $\alpha$ between 2 and 3.</figcaption></figure>

In Figure 2, we observe that $\theta_{\gamma}$, $\theta_{\kappa}$ and $\theta_{I}$ exhibit qualitatively the same behaviour, with stable values over a wide range of $d_{min}/d_{Max}$ ratios.

# 4 Constitutive model

## 4.1 Thermodynamic considerations

We introduce a constitutive model to describe and predict the behaviour of the crushable granular medium. The model is drawn from Nguyen and Einav (2009) and Einav (2007a,b), modified for consistency with the Cosserat continuum.

The model is constructed in the framework of hyperplasticity (Houlsby and Puzrin, 2007), so the first step is to make a statement encapsulating the First and Second laws of Thermodynamics:

$$W = \dot{\Psi} + \Phi, \quad \Phi \ge 0. \qquad (44)$$

Only the Helmholtz free energy is a proper state function, so its time rate is denoted with a superposed dot. The rate of work and rate of dissipation are definable only as rates, and so are left without a dot to denote this conceptual distinction.

For the time being we take the general form of $\Psi$ given in (25), and leave the explicit specification until later.

## 4.2 Dissipation

The total dissipation rate can be written as the sum of the products of the dissipative variables and their respective conjugate stresses, i.e.

$$\Phi(\dot{B}, \dot{\epsilon}_v^p, \dot{\gamma}^p) = E_B \dot{B} + p \dot{\epsilon}_v^p + q \dot{\gamma}^p. \qquad (45)$$

We assume that this dissipation can be given by a specific pseudo-potential, which takes the following specific form:

$$\Phi(\dot{B}, \dot{\epsilon}_v^p, \dot{\gamma}^p) = \sqrt{[D_B(\dot{B})]^2 + [D_v(\dot{\epsilon}_v^p)]^2 + [D_s(\dot{\gamma}^p)]^2}, \qquad (46)$$

where each of $D_B$, $D_v$ and $D_s$ are homogeneous first order functions of the dissipative rates.

The Legendre transformation for first order homogeneous functions (in the rates of the internal variables) is used to pass from the dissipation pseudo-potential to the yield surface in dissipative stress space. After some re-arrangement, we obtain:

$$y_{\Phi}(E_B, p, q) = \left( \frac{E_B}{\partial D_B(\dot{B}) / \partial \dot{B}} \right)^2 + \left( \frac{p}{\partial D_v^p(\dot{\epsilon}_v^p) / \partial \dot{\epsilon}_v^p} \right)^2 + \left( \frac{q}{\partial D_s^p(\dot{\gamma}^p) / \partial \dot{\gamma}^p} \right)^2 - 1 \le 0. \qquad (47)$$

The mathematical form of the dissipation potentials can now be specified:

$$D_B(\dot{B}) = \frac{\sqrt{E_B E_c}}{(1 - B) \cos(\omega)} \dot{B}, \qquad (48)$$

$$D_v^p(\dot{\epsilon}_v^p) = \frac{p}{(1 - B) \sin(\omega)} \sqrt{\frac{E_c}{E_B}} \dot{\epsilon}_v^p, \qquad (49)$$

$$D_s^p(\dot{\gamma}^p) = Mp|\dot{\gamma}^p|, \qquad (50)$$
where $E_c$ is the critical breakage energy, $M$ is the critical state parameter satisfying $q = Mp$ in triaxial compression and $\omega$ is a coupling angle between dissipation due to breakage and dissipation due to plastic volumetric straining, as per Einav (2007b) and Nguyen and Einav (2009).

Following the Legendre transformation, the evolution of the breakage state variable $B$, and the rates of the plastic strains and plastic curvatures are given by the flow rules:

$$ \dot{B} = \lambda \frac{\partial y_{\Phi}}{\partial E_B} = \lambda \frac{2E_B}{\left[ \partial D_B(\dot{B}) / \partial \dot{B} \right]^2} = \lambda \frac{2(1-B)^2 \cos^2(\omega)}{E_c}, \quad (51) $$

$$ \begin{aligned} \dot{\gamma}_{ij}^p &= \lambda \frac{\partial y_{\Phi}}{\partial \tau_{ij}} = \lambda \left( \frac{\partial y_{\Phi}}{\partial p} \frac{\partial p}{\partial \tau_{ij}} + \frac{\partial y_{\Phi}}{\partial q} \frac{\partial q}{\partial \tau_{ij}} \right) = \lambda \left( \frac{2p}{\left[ \partial D_v^p (\dot{\varepsilon}_v^p) / \partial \dot{\varepsilon}_v^p \right]^2} \frac{1}{3} \delta_{ij} + \frac{2q}{\left[ \partial D_s^p (\dot{\gamma}_s^p) / \partial \dot{\gamma}_s^p \right]^2} \frac{1}{q} (h_1^* s_{ij} + h_2^* s_{ji}) \right) \\ &= \lambda \left( \frac{2E_B (1-B)^2 \sin^2(\omega)}{3pE_c} \delta_{ij} + \frac{2(h_1^* s_{ij} + h_2^* s_{ji})}{(Mp)^2} \right), \end{aligned} \quad (52) $$

$$ \dot{\kappa}_{ij}^p = \lambda \frac{\partial y_{\Phi}}{\partial \mu_{ij}} = \lambda \frac{\partial y_{\Phi}}{\partial q} \frac{\partial q}{\partial \mu_{ij}} = \lambda \frac{2q}{\left[ \partial D_s^p (\dot{\gamma}_s^p) / \partial \dot{\gamma}_s^p \right]^2} \frac{1}{\ell^{e2}q} (h_3^* m_{ij} + h_4^* m_{ji}) = \lambda \frac{2(h_3^* m_{ij} + h_4^* m_{ji})}{\ell^{e2}(Mp)^2}, \quad (53) $$

where $h_1, h_2, h_3$ and $h_4$ are weighting factors discussed and defined in §4.6. It should be noted that the yield surface depends only on the trace part of the couple-stress tensor, simplifying the full expansion of the flow rule in the curvatures (53) by allowing us to set the trace part of the derivatives to zero.

As a matter of convenience, a change from the yield surface in dissipative stress space to the yield surface in mixed stress-breakage-energy space can be achieved by directly substituting (48), (49) and (50) into (47) to obtain

$$ y_{mix} = \frac{E_B (1-B)^2}{E_c} + \left( \frac{q}{Mp} \right)^2 - 1 \le 0. \quad (54) $$

It is then possible to express $E_B$ in terms of $p$ and $q$, so that instead of $y = y_{mix}(E_B, p, q, B)$ one may write $y = y_{pq}(p, q, B)$. The flow rules are associated with $y_{\Phi}$ but not with $y$ that is often drawn experimentally. Hence, they will appear non-associated if viewed only in the $p-q$ stress space.

We note here that the definitions of $E_B$ and $q$ are expanded compared to their classical formulation due to the addition of Cosserat effects. $M$ is also calculated with respect to the Cosserat, rather than classical, invariants (Mühlhaus and Vardoulakis, 1987) that are further defined in §4.6. When these Cosserat effects are induced during loading, yield will occur at lower stress values due to the presence of couple-stresses. However, there is no effect during isotropic compaction, meaning the definition of $E_c$ remains unchanged.

## 4.3 Elasticity

The elastic behaviour can be recovered by specifying the Helmholtz free energy density in the reference grain fraction so that linear Cosserat elasticity is recovered from the generic potential (25). This can be specified by

$$ \Psi = \frac{1 - \theta_{\gamma} B}{2} C_{ijkl}^{e} \gamma_{ij}^{e} \gamma_{kl}^{e} + d_{Max} \frac{1 - \theta_{\kappa} B}{2} D_{ijkl}^{e} \kappa_{ij}^{e} \kappa_{kl}^{e}, \quad (55) $$

where

$$ C_{ijkl}^{e} = \left( K - \frac{2}{3} G \right) \delta_{ij} \delta_{kl} + (G + G_c) \delta_{ik} \delta_{jl} + (G - G_c) \delta_{il} \delta_{jk}, \quad (56) $$

$$ D_{ijkl}^{e} = \left( L - \frac{2}{3} H \right) \delta_{ij} \delta_{kl} + (H + H_c) \delta_{ik} \delta_{jl} + (H - H_c) \delta_{il} \delta_{jk}. \quad (57) $$

$K$ and $G$ are the bulk and shear elastic moduli as per the classical continuum, and $G_c, L, H$ and $H_c$ are additional elastic parameters introduced for the Cosserat continuum. As previously noted in §3.1, the parameters $L, H$, and $H_c$ have units of stress in this work.
Now the thermodynamic stresses can be derived from the Helmholtz free energy (55):

$$ \tau_{ij} = \frac{\partial \Psi}{\partial \gamma_{ij}^e} = (1 - \theta_{\gamma B}) C_{ijkl}^e \gamma_{kl}^e, \qquad (58) $$

$$ \mu_{ij} = \frac{\partial \Psi}{\partial \kappa_{ij}^e} = d_{Max}^2 (1 - \theta_{\kappa B}) D_{ijkl}^e \kappa_{kl}^e, \qquad (59) $$

$$ E_B = -\frac{\partial \Psi}{\partial B} = \frac{\theta_{\gamma}}{2} C_{ijkl}^{e} \gamma_{ij}^{e} \gamma_{kl}^{e} + \frac{d_{Max}^{2} \theta_{\kappa}}{2} D_{ijkl}^{e} \kappa_{ij}^{e} \kappa_{kl}^{e}. \qquad (60) $$

## 4.4 Plastic multiplier

The model is completely specified by determining an explicit expression for the plastic multiplier $\lambda$ by applying the consistency condition, i.e. when the stress state is on the yield surface $y = 0$, so for any further inelastic loading $\dot{y} = 0$. By standard processes this allows a series of substitutions and rearrangements that results in an expression for $\lambda$, and since the process is quite laborious and the expression long, we do not write it explicitly here.

For the purposes of performing numerical integrations, a smoothing can be achieved by implementation in the framework of $h^2$ plasticity (Einav, 2012). This framework allows the evolution equations to be smooth and continuous for numerical purposes. Some small modifications to the flow rules (51), (52) and (53) result in:

$$ \dot{B} = \langle \lambda \rangle (1 + y) \xi \frac{2(1 - B)^2 \cos^2(\omega)}{E_c}, \qquad (61) $$

$$ \dot{\gamma}_{ij}^p = \langle \lambda \rangle (1 + y) \xi \left( \frac{2E_B(1 - B)^2 \sin^2(\omega)}{3pE_c} \delta_{ij} + \frac{2(h_1^* s_{ij} + h_2^* s_{ji})}{(Mp)^2} \right), \qquad (62) $$

$$ \dot{\kappa}_{ij}^p = \langle \lambda \rangle (1 + y) \xi \frac{2(h_3^* m_{ij} + h_4^* m_{ji})}{\ell^{e2} (Mp)^2}, \qquad (63) $$

where $\langle \cdot \rangle$ are Macaulay brackets, $y$ is the yield function (any expression of $y$ will work as they all have the same scalar value), and $\xi$ is a parameter controlling the instantaneous stiffness of the response. As $\xi \to \infty$, the classical hyper-plastic model is recovered. Smaller values of $\xi$ represent “softer” behaviours with dissipation being more significantly activated at lower values of load.

## 4.5 Incremental constitutive response

In order to complete the description of the model, the expression for the incremental stiffness relationship is required. Even though the expressions for $\tau_{ij}$ and $\mu_{ij}$ are independent of $\kappa_{ij}^e$ and $\gamma_{ij}^e$ respectively, there is coupling present in the incremental relationship due to $\lambda$ being a function of both $\dot{\gamma}_{kl}$ and $\dot{\kappa}_{kl}$. A linear decomposition of $\lambda$ is performed to obtain

$$ \lambda = \lambda_{kl}^{\gamma} \dot{\gamma}_{kl} + \lambda_{kl}^{\kappa} \dot{\kappa}_{kl}. \qquad (64) $$

Using (64), the incremental constitutive relationship under increasing load in a compact differential form is given by:

$$ \begin{bmatrix} \dot{\tau}_{ij} \\ \dot{\mu}_{ij} \end{bmatrix} = \begin{bmatrix} \frac{\partial \tau_{ij}}{\partial \gamma_{kl}^e} - \lambda_{kl}^\gamma (1+y)^\xi \left( \frac{\partial y_\Phi}{\partial \tau_{ij}} - \frac{\partial \tau_{ij}}{\partial B} \frac{\partial y_\Phi}{\partial E_B} \right) & -\lambda_{kl}^\kappa (1+y)^\xi \left( \frac{\partial y_\Phi}{\partial \tau_{ij}} - \frac{\partial \tau_{ij}}{\partial B} \frac{\partial y_\Phi}{\partial E_B} \right) \\ -\lambda_{kl}^\gamma (1+y)^\xi \left( \frac{\partial y_\Phi}{\partial \mu_{ij}} - \frac{\partial \mu_{ij}}{\partial B} \frac{\partial y_\Phi}{\partial E_B} \right) & \frac{\partial \mu_{ij}}{\partial \kappa_{kl}^\kappa} - \lambda_{kl}^\kappa (1+y)^\xi \left( \frac{\partial y_\Phi}{\partial \mu_{ij}} - \frac{\partial \mu_{ij}}{\partial B} \frac{\partial y_\Phi}{\partial E_B} \right) \end{bmatrix} \begin{bmatrix} \dot{\gamma}_{kl} \\ \dot{\kappa}_{kl} \end{bmatrix}. \qquad (65) $$

It is convenient to write this relationship as

$$ \begin{bmatrix} \dot{\tau}_{ij} \\ \dot{\mu}_{ij} \end{bmatrix} = \begin{bmatrix} E_{ijkl}^{ep} & F_{ijkl}^{ep} \\ K_{ijkl}^{ep} & M_{ijkl}^{ep} \end{bmatrix} \begin{bmatrix} \dot{\gamma}_{kl} \\ \dot{\kappa}_{kl} \end{bmatrix}. \qquad (66) $$

## 4.6 Invariants

The formulation of the first stress and strain invariants in a Cosserat continuum is the same as for the classical continuum, as outlined above in §2.1 and §2.3. However, the formulations of the second deviatoric stress and strain invariants are not trivial in the Cosserat continuum. Micro-mechanical arguments have been advanced by Mühlhaus and Vardoulakis (1987), Rattez, Stefanou and Sulem (2018) and Vardoulakis and Sulem (1995), allowing the invariant for plastic shear strain rate to be defined by consideration of kinematics, or the invariant for the deviatoric stress to be defined by consideration of statics. Comparison with DEM simulations has shown that the static model is more
suitable for two dimensional materials, while the kinematic is more suitable for a true three dimensional granular material. The two invariants are respectively

$$ \dot{\gamma}^p = \sqrt{g_1^* \dot{e}_{ij}^p \dot{e}_{ij}^p + g_2^* \dot{e}_{ij}^p \dot{e}_{ji}^p + \ell^\gamma (g_3^* \dot{k}_{ij}^p \dot{k}_{ij}^p + g_4^* \dot{k}_{ij}^p \dot{k}_{ji}^p)}, \quad (67) $$

$$ q = \sqrt{h_1^* s_{ij} s_{ij} + h_2^* s_{ij} s_{ji} + \frac{1}{\ell^{q2}} (h_3^* m_{ij} m_{ij} + h_4^* m_{ij} m_{ji})}. \quad (68) $$

Certain weighting factors are given to the values of $g_i^*$ in the kinematic expression of $\dot{\gamma}^p$, and $h_i^*$ in the static expression of $q$, in order to recover the classical expressions when the Cosserat effects vanish. Two pairs of values are defined, $q - \dot{\gamma}^p$ kinematic and $q - \dot{\gamma}^p$ static, with different weighting factors $g_i^*$ and $h_i^*$. The length scale $\ell^\gamma$ is identified with the grain size in the micro-mechanical derivations. It is essential to note here that the definition of these invariants in the work of Vardoulakis and co-workers implicitly takes $\ell^\gamma = \ell^q = R$, that is to say that both lengths are equal to some constant multiple of the mean grain radius (Mühlhaus and Vardoulakis, 1987; Vardoulakis and Sulem, 1995). As a consequence, the derivation of the weighting factors relies on all particles having the same diameter. In a medium featuring grain breakage, this will certainly not be the case, and so the invariants must be at least partially re-visited to determine appropriate expressions for the medium.

By relaxing the assumption that the internal length is the same for the stress and strain invariants, but retaining the assumption of coaxiality (that is to say that the stress and plastic strain rates act in the same direction), and taking the weighting factors (either $h_i^*$ or $g_i^*$) on an ad-hoc basis, new expressions can be determined by appropriate substitution of the flow rules and simplification. We leave open the possibility that $\ell^\gamma \neq \ell^q$, and we identify $\ell^q$ with $\ell^e$, defined as per (27). However, there is no immediately obvious identification to make that will give an explicit expression for $\ell^\gamma$. As such, it is supposed that it may also evolve with breakage in a similar way to the elastic length, and the ratio between the elastic and plastic length is a fixed ratio $\zeta$, which acts as a new material property. It should also be noted that the definition of $q$ used in previous examples of Breakage Mechanics is consistent with $q = \sigma_a - \sigma_r$, that is to say with the difference between the axial and radial stresses in a triaxial test. This leads to a tensorial definition of

$$ q = \sqrt{\frac{3}{2} s_{ij} s_{ij}} \quad (69) $$

in the classical continuum, which the definition must collapse to in the absence of Cosserat effects. A table of weighting factors for the invariants (67) and (68) is arrived at:

<table><caption>Table 1: Values of the coefficients for stress and plastic strain rate invariants in a Cosserat continuum made consistent with Breakage Mechanics</caption><thead><tr><th></th><th>2D model</th><th>3D model</th></tr></thead><tbody><tr><td>Static model</td><td>$g_i^* = \{1/2, 1/6, \zeta^2/3, 0\}$<br>$h_i^* = \{9/4, -3/4, 3, 0\}$</td><td>$g_i^* = \{8/15, 2/15, 8\zeta^2/15, 2\zeta^2/15\}$<br>$h_i^* = \{2, -1/2, 2, -1/2\}$</td></tr><tr><td>Kinematic model</td><td>$g_i^* = \{1, -1/3, 4/3, 0\}$<br>$h_i^* = \{9/8, 3/8, 3\zeta^2/4, 0\}$</td><td>$g_i^* = \{8/9, -2/9, 8/9, -2/9\}$<br>$h_i^* = \{6/5, 3/10, 6\zeta^2/5, 3\zeta^2/10\}$</td></tr></tbody></table>

The values in Table 1 are compatible with previous works (Mühlhaus and Vardoulakis, 1987; Sulem and Vardoulakis, 1990), adjusted appropriately to account for the slightly different definition of the stress invariant.

# 5 Parameter values

## 5.1 Cosserat stiffness relationships

Before specifying the parameter values, it is convenient to derive some relationships between the additional Cosserat elastic and the classical elastic moduli. The essence of this process is to set $B = 0$ (as we are deriving the stiffness of unbroken material), perform a Legendre transformation on (55) to the corresponding Gibbs free energy, and then make the constitutive assumption that the Gibbs free energy can also be expressed in terms of the stress invariants:

$$ \Gamma = -\frac{1}{2} \left( \frac{p^2}{K} + \frac{q^2}{3G} \right). \quad (70) $$

The compliance relationships derived from this form are expressed in terms of $K, G, h_1^*$ and $h_2^*$, and may be compared with the relationships obtained from the Legendre transformation, which are expressed in terms of $K, G$ and $G_c$. This
allows the derivation of a formula for the value of $G_c$. We compare the results considering the 1212 and 2112 entries of the two compliance relationships to ultimately obtain:

$$G_c = \frac{3G}{2(h_1^* - h_2^*)} \qquad (71)$$

Expressions are derived in a similar way for the couple-stress and curvature compliance relationships, leading ultimately to

$$H = \frac{3G}{2(h_3^* + h_4^*)}, \qquad (72)$$

and

$$H_c = \frac{3G}{2(h_3^* - h_4^*)}. \qquad (73)$$

Finally, it should be noted that the Gibbs free energy expression excludes any contribution from the trace part of the elastic curvatures, similarly to how the flow rules for the plastic curvatures ignore the trace component. This is equivalent to stating that there are no twisting motions contributing to the stored energy, and hence the material parameter $L$ makes no contribution to the model.

## 5.2 Calibration of the model parameters

In order to apply the model to a real geomaterial, certain experiments must be carried out in order to determine the values of the various parameters in the model. These parameters are the bulk stiffness $K$, the shear stiffness $G$, the Cosserat elastic moduli $G_c$, $H$ and $H_c$, the maximum grain size $d_{Max}$, the grading indices $\theta_\gamma$, $\theta_\kappa$ and $\theta_I$, the critical state parameter $M$, the crushing energy $E_c$, and the coupling angle $\omega$.

By measuring the grain size distribution, $d_{Max}$ is immediately obtained as the largest grain size present. Either the particular material can be placed in some sort of apparatus that allows arbitrarily large shear, in order to determine an ultimate grading, or one may select values of the power law coefficient from the literature. Hence, the full grain size distribution allows calculation of $\theta_\gamma$, $\theta_\kappa$ and $\theta_I$ by using the universal initial distribution and the ratio of the areas under the initial and ultimate grain size distributions calculated to give an initial value of $B \ge 0$.

The density of the material can be easily measured using standard techniques (i.e. measuring the mass of a known volume of material), or if the composition of the grains is known, their solid density divided by the initial solid volume fraction will give the bulk density.

The bulk and shear moduli $K$ and $G$ can be measured using standard techniques such as shear wave measurements. While in reality, these moduli are pressure dependent, we assume isotropic linear elasticity for this model, meaning that the calibration should be performed at a representative pressure. We may calculate $G_c$, $H$, $H_c$ from these values after selecting the kinematic or static model and the choice of dimension.

The critical state parameter $M$ is the slope of the critical state line in $p-q$ space, typically obtained in a triaxial test.

The critical breakage energy can be found by subjecting a sample to isotropic compression until grain crushing starts. By rearranging (54) with $q = 0$, an expression for $E_c$ is straightforwardly obtained. It should be noted that this value can also be inferred from scaling laws based on the maximum grain size $d_{Max}$ (Zhang, Buscarnera and Einav, 2016).

The only parameter of the model which we adjust to fit the plastic response is the coupling angle $\omega$. Calibration may be performed by seeing which value of $\omega$ best matches the stress-strain response after plasticity. In the model by Tengattini, Das and Einav (2016), this freedom was resolved, but this is not implemented in the current work for simplicity.

Noting the above general principles for model calibration, for this paper we take the values of the Breakage model that we explore from those used in Das, Nguyen and Einav (2011) for modelling Bentheim sandstones, which in turn were calculated from Baud, Klein and Wong (2004), Baud, Vajdova and Wong (2006), Tembe, Baud and Wong (2008) and Wong, Baud and Klein (2001). Following these papers, we also set the ultimate fractal parameter $\alpha = 2.6$ and take 0.105 mm as the maximum grain size $d_{Max}$. Consistent with Buscarnera and Einav (2012) and Kendall (1978) and the arguments therein, we set the minimum grain size to 0.001 mm. These choices set the values for $\theta_\gamma$, $\theta_\kappa$ and $\theta_I$. We have additional sets of model parameters associated with our re-framing of the model in the Cosserat continuum, which are calculated as per §5.1. The kinematic model is selected for the reasons discussed in §4.6, and the model parameter relating the elastic and plastic length scales is set $\zeta = 1$. We treat the $h^2$ parameter $\xi$ (see §4.4) merely
as a numerical smoothing parameter that allows the system to be numerically integrated more easily. Hence, we take $\xi = 100$ so that we are close to the response of the underlying hyper-plastic model. Finally, the total set of physical parameters in units of grams/millimetres/milliseconds is summarised in Table 2:

<table><caption>Table 2: List of material constants.</caption><thead><tr><th colspan="2">Mechanical Parameters</th><th colspan="2">Granulometry Parameters</th></tr></thead><tbody><tr><td>K</td><td>13833</td><td>MPa</td><td>θ<sub>γ</sub></td><td>0.80</td></tr><tr><td>G</td><td>7588</td><td>MPa</td><td>θ<sub>κ</sub></td><td>0.89</td></tr><tr><td>G<sub>c</sub></td><td>12646.7</td><td>MPa</td><td>θ<sub>I</sub></td><td>0.91</td></tr><tr><td>H</td><td>7588</td><td>MPa</td><td>d<sub>Max</sub></td><td>0.105 mm</td></tr><tr><td>H<sub>c</sub></td><td>12646.7</td><td>MPa</td><td>ρ</td><td>2.5 × 10<sup>-3</sup> g/mm<sup>3</sup></td></tr><tr><td>E<sub>c</sub></td><td>4.65</td><td>MPa</td><td></td><td></td></tr><tr><td>M</td><td>1.7</td><td></td><td></td><td></td></tr><tr><td>ω</td><td>70</td><td>°</td><td></td><td></td></tr></tbody></table>

# 6 Linear stability analysis

## 6.1 Method

We assume that there is a homogeneous solution obeying the governing equations, and perform a perturbation about said solution. We may re-write the kinematic fields as

$$u(x_i, t) = u^0(x_i, t) + \tilde{u}(x_i, t), \qquad (74)$$

$$\omega^c(x_i, t) = \omega^{c0}(x_i, t) + \tilde{\omega}^c(x_i, t), \qquad (75)$$

where the 0 superscript denotes the homogeneous solution and the tilde the perturbation term. We assume that the resultant perturbations in stresses, strains, couple-stresses and curvatures follow the incremental constitutive relationships (66), linearised around the homogeneous solution.

The governing equations are satisfied by the homogeneous field, and hence by subtraction they must also be satisfied by the perturbation terms. The momentum balance equations for the Cosserat continuum are given by:

$$\tau_{ij,j} - \rho \frac{\partial^2 u_i}{\partial t^2} = 0, \qquad (76)$$

$$\mu_{ij,j} - \epsilon_{ijk} \tau_{jk} - I \frac{\partial^2 \omega_i^c}{\partial t^2} = 0. \qquad (77)$$

In order to obtain an appropriate linear system allowing a bifurcation analysis, we insert the incremental constitutive relationships (66) and the definitions of the strain and curvature rates (1) and (2) into the perturbations of (76) and (77) to obtain

$$E_{jklm}^{ep}(\tilde{u}_{l,mk} + \epsilon_{lmo}\tilde{\omega}_{o,k}^{c}) + F_{jklm}^{ep}\tilde{\omega}_{l,mk}^{c} - \rho \frac{\partial^2 \tilde{u}_j}{\partial t^2} = 0, \qquad (78)$$

$$K_{jklm}^{ep}(\tilde{u}_{l,mk} + \epsilon_{lmo}\tilde{\omega}_{o,k}^{c}) + M_{jklm}^{ep}\tilde{\omega}_{l,mk}^{c} - \epsilon_{jkl}(E_{klmo}^{ep}(\tilde{u}_{m,o} + \epsilon_{mop}\tilde{\omega}_{p}^{c}) + F_{klmo}^{ep}\tilde{\omega}_{m,o}^{c}) - I \frac{\partial^2 \tilde{\omega}_j^c}{\partial t^2} = 0. \qquad (79)$$

We now have a set of 6 equations with 6 unknowns, $\tilde{u}_j$ and $\tilde{\omega}_j^c$. The above linear system has solutions of the following form:

$$U_j^*(x_k, t) = \bar{U}_j e^{st + \frac{2\pi i}{\Lambda} x_k n_k}, \qquad (80)$$

$$\Omega_j^*(x_k, t) = \bar{\Omega}_j e^{st + \frac{2\pi i}{\Lambda} x_k n_k}, \qquad (81)$$

where $\bar{U}_j$ and $\bar{\Omega}_j$ are vectors of coefficients, $s$ is the growth coefficient in time of the perturbation (also known as the Lyapunov exponent), $i = \sqrt{-1}$ is the imaginary unit, $\Lambda$ is the perturbation wavelength, $n_k$ the polarisation direction and $U_j^*(x_k, t)$ and $\Omega_j^*(x_k, t)$ are vectors that contain the unknown perturbation fields. We can insert (80) and (81)
into (78) and (79) to obtain

$$ \left[ E_{jklm}^{ep} \left( \frac{4\pi^2 i^2}{\Lambda^2} n_m n_k \bar{U}_l + \epsilon_{lmo} \frac{2\pi i}{\Lambda} n_k \bar{\Omega}_o \right) + F_{jklm}^{ep} \frac{4\pi^2 i^2}{\Lambda^2} n_m n_k \bar{\Omega}_l - \rho s^2 \bar{U}_j \right] e^{st + \frac{2\pi i}{\Lambda} x_p n_p} = 0, \quad (82) $$

$$ \left[ K_{jklm}^{ep} \left( \frac{4\pi^2 i^2}{\Lambda^2} n_m n_k \bar{U}_l + \epsilon_{lmo} \frac{2\pi i}{\Lambda} n_k \bar{\Omega}_o \right) + M_{jklm}^{ep} \frac{4\pi^2 i^2}{\Lambda^2} n_m n_k \bar{\Omega}_l - \epsilon_{jkl} \left( E_{klmo}^{ep} \left( \frac{2\pi i}{\Lambda} n_o \bar{U}_l + \epsilon_{mop} \bar{\Omega}_p \right) + F_{klmo}^{ep} \frac{2\pi i}{\Lambda} n_o \bar{\Omega}_m \right) - Is^2 \bar{\Omega}_j \right] e^{st + \frac{2\pi i}{\Lambda} x_q n_q} = 0. \quad (83) $$

We may re-write (82) and (83) as a matrix acting on the coefficients $\bar{U}_j$ and $\bar{\Omega}_j$. We now find when this matrix becomes singular, which results in a characteristic polynomial in $s$.

We explore how the roots of the polynomial change as a function of the changing parameters and state variables of the system, as well as the wavelength of the perturbation $\Lambda$. We may characterise the roots as being one of four types, depending on whether $s$:

1. has negative real part. This indicates that the perturbation decreases in time and the homogeneous system is stable.

2. is complex with a positive real part. This indicates a flutter instability where the system will cycle between loading and unloading. However, if the constitutive model consists of an elastic and an elasto-plastic region, the unloading will result in a return to the elastic region. It is hence unclear if we may conclude anything from a result of this type.

3. is real with positive real part. This indicates that the perturbation grows in time and the homogeneous system is unstable. Hence, we have a bifurcation.

4. has real part equal to zero. We cannot conclude anything in this case.

For a given set of state variables, we may search over a range of wavelengths and orientations to find the largest real positive $s$, which will be the fastest growing perturbation. The wavelength $\Lambda_{max}$ associated with this perturbation is interpreted as the width of the localisation band. $\Lambda_{min}$ is the wavelength where the real part of $s = 0$, and any perturbation with a smaller wavelength than this cannot propagate. These points are illustrated in Figure 3.

<figure><img src="image_3.png" /><figcaption>Figure 3: A typical example of the behaviour of the Lyapunov exponent $s$ versus the perturbation wavelength $\Lambda$: (a) homogeneous deformations are stable, and (b) homogeneous deformations are unstable, i.e. strain localisation occurs. The points defining $\Lambda_{min}$, the minimum wavelength that can propagate, and $\Lambda_{max}$, the fastest growing wavelength, are highlighted.</figcaption></figure>

We consider a reference state corresponding to homogeneous deformation, so that all Cosserat terms are zero. As a consequence, the bifurcation matrix is significantly simplified, as setting $\kappa_{ij}^e = 0$ results in $F_{ijkl}^{ep} = K_{ijkl}^{ep} = 0$, and only the elastic component of $M_{ijkl}^{ep}$ remains. This reduces the equation we seek to solve to the determinant:

$$ \left| E_{jklm}^{ep} \frac{4\pi^2 i^2}{\Lambda^2} n_m n_k - \rho s^2 \delta_{jl} - \epsilon_{jkl} E_{klmo}^{ep} \frac{2\pi i}{\Lambda} n_o \right| M_{jklm}^{ep} \frac{4\pi^2 i^2}{\Lambda^2} n_m n_k \delta_{lp} - \epsilon_{jkl} E_{klmo}^{ep} \epsilon_{mop} - Is^2 \delta_{jp} = 0. \quad (84) $$

It should also be noted that beyond the point of initial bifurcation, the homogeneous solution we assume for the linear stability analysis is no longer valid. However, by continuing to integrate the model (assuming homogeneity) beyond the bifurcation point, we can explore the tendencies of the system close to the homogeneous state.
We then simulate an individual load path using a material model implemented in $h^2$plasticity, and pass the set of state variables at each step in time to the analysis code, which returns whether there is a localisation, and if so the corresponding minimum and maximum wavelengths.

## 6.2 System analysis

In this paper we focus on shearing under constant volume, which corresponds to fast undrained shear that can be seen in seismogenic faults. While our model does not feature any hydraulic coupling, it is still important to understand how the material model behaves under these conditions.

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

Figure 4: A schematic of the considered system. We apply a total displacement of $\Delta u_1$ in the $x_1$ direction that acts as a pre-confinement, marked in blue. We then start our test at $t = 0$ and apply a displacement $\Delta u_2$ in the $x_2$ direction to induce shearing, marked in red. During shearing we maintain constant volume (no further $\Delta u_1$ is applied). The system is invariant in the $x_3$ direction.

The system that we consider is illustrated in Figure 4. The first step is using the LSA to understand the behaviour of the system in a general sense. In order to do this, we examine where on the yield surface can support bifurcation. This is a qualitative description that gives an indication of the behaviour. If we limit our analysis to consider only the pairs of elastic strains $\gamma_{11}^e$ and $\gamma_{21}^e$ that place the system described above on the yield surface, there appears to be no bifurcation supported by the system. This behaviour can ultimately be traced to the fact that $\gamma_{22/33}^e = 0$ during elastic loading.

However, if we were to imagine following an individual load path, $\gamma_{22/33}^e$ will start to accumulate as we enter the plastic region of the load path, which ultimately allows the system to bifurcate towards a state of non-homogeneous deformation. In the example this corresponds to the nucleation of a shear band. To represent this behaviour, when analysing the localisation regions of the yield surface we consider systems with isotropic values of the elastic strains. This allows the determination of which parts of the yield surface can in principle support bifurcation and a shear band parallel to the layer. Thus, the calculation steps are:

1. Fix a value of $B$ (thus determining the yield surface),

2. Set an isotropic value of $\gamma_{11/22/33}^e$ less than or equal to $\gamma_{p_{crit}}^e$, the isotropic elastic strain that corresponds to the critical pressure for that value of $B$,

3. Calculate the shearing strain $\gamma_{21}^e$ that places the system on the yield surface,

4. Check whether for a very large wavelength ($\Lambda \gg 20d_{Max}$) there is a value of $s$ with positive real part.

If the condition in step 4 is fulfilled, there is a bifurcation and the system will form a shear band. We check only a large wavelength as this is computationally more efficient than determining the value of $\Lambda_{max}$ for every stress state. We check only the horizontal orientation $n_i = 1, 0, 0$ in step 4, as for these load conditions this is the fastest growing orientation.
<figure><img src="image_5.png" /><figcaption>Figure 5: The regions of the yield surface that can support localisation parallel to the shear direction when sheared under constant volume conditions, using the calibrated values of the material parameters. The analysis is performed for $B = 0, 0.1, ..., 0.9$, with the corresponding yield surfaces above. We place a marker to highlight the maximum confining stress capable of supporting localisation for each value of $B$.</figcaption></figure>

In Figure 5 we can observe that the region that supports shear band formation for the calibrated values consistently occurs on the low confining pressure side of the yield surfaces, where plastic shear straining is the dominant form of dissipation.

We can now follow an individual load path to observe the favoured localisation width and the predicted post-localisation behaviour. For this analysis we apply the confinement as described in Figure 4, rather than the isotropic confinement we have used to study the localisation regions of the yield surface. Thus, we impose $\gamma_{11}$ while setting $\gamma_{22} = \gamma_{33} = 0$. The reference state of the material is assumed homogeneous and to act as if in the classical continuum. Thus, the applied shear motion in the $u_2$ direction will result in symmetric strains $\gamma_{12}$ and $\gamma_{21}$.

<figure><img src="image_6.png" /><figcaption>Figure 6: The system behaviour at constant volume with increasing shear strain. (a) The shear stress $\tau_{21}$, and (b) the minimum wavelength for perturbation growth $\Lambda_{min}$ and the fastest growing wavelength $\Lambda_{max}$. Both wavelengths have been normalised by the largest grain size $d_{Max}$. We place a marker on the minimum value of $\Lambda_{max}$ for clarity.</figcaption></figure>

In Figure 6 we observe that a shear band forms during shear at constant volume for the calibrated values we have chosen. We interpret the minimum value attained by the $\Lambda_{max}$ curve to be the predicted thickness of the shear band in these conditions. It should be noted that the shear bands form at $\gamma_{21} \approx 0.02$, which is more than twice the strain associated with the peak of the $\tau_{21} - \gamma_{21}$ curve. We can define a ratio $R$, which is the ratio of the confining stress when the system first enters the bifurcation region $p_{yield}$ to the maximum confining stress that can support bifurcation $p_{max}$:

$$R = \frac{p_{yield}}{p_{max}} \qquad (85)$$

In this case $R = 0.509$. While the stress path enters the bifurcation region illustrated in Figure 5 as soon as plasticity begins, a certain amount of elastic straining in the 22 and 33 entries must be accumulated before the system localises into a shear band. These elastic strains are equal and opposite to the accumulating plastic strains in those directions, in order to maintain the constant volume condition.
In the $\tau_{21} - \gamma_{21}$ curve we observe some apparent softening of the system. We emphasise that this softening is only apparent, and that the system is actually hardening in the classical sense (the yield surface is expanding).

<figure><img src="image_7.png" /><figcaption>Figure 7: A comparison of the $p-q$ points that lie on the initial yield surface for $B = 0$, and for $B = 0.217$ (corresponding to $\gamma_{21} = 0.2$), with the load path of the curve in Figure 6. We place a marker to highlight the maximum confining stress capable of supporting localisation for each value of $B$.</figcaption></figure>

In Figure 7 we demonstrate the contrast between the apparent softening and the true hardening. For the sake of clarity we have continued the load path integration to $\gamma_{21} = 0.2$ to show the growth of the yield surface. This behaviour arises as the growth of the breakage variable $B$ contributes to an expanded yield surface (that is to say, the system experiences true hardening), while also reducing the effective stiffness of the material, leading to the apparent softening of the system. We also note that stress paths in $p-q$ space such as shown in Figure 7 are not entirely realistic for typical strong material, and are a legacy of our simplifying choice not to include porosity in our analysis. Within our model, there is a competition between the effect of the rate of breakage and the effect of the rate of elastic strain, in particular in the 11 direction. If we increase the initial value of $B$ the apparent softening decreases as the rate of breakage growth is slower (illustrated in Figure 9). However, as we increase the value of $\omega$, despite the decreasing rate of breakage growth, the system softens faster (illustrated in Figure 11) as the elastic strain in the 11 direction decreases in order to accommodate the increasing plastic strain.

In a model with a more physically accurate porosity-breakage coupling, we would expect a less dramatic apparent softening as the porosity saturates to a minimum value and the system enters a stable critical state. Despite the hardening that can arise from dilative behaviour, this would not necessarily prevent a model with porosity from demonstrating localisation behaviour over as wide a range of pressures as our current model, given Figure 7 shows that our model supports localisation during hardening. This is particularly the case if the model with porosity were to be sheared under constant volume, which would suppress the hardening due to dilation.

As well as the calibrated values of the system, we can explore how the system changes when we change some of the parameters related to Breakage Mechanics. In particular, we are interested in how the behaviour changes as we change the maximum grain size, the initial value of the breakage, and the tendency of the system to dissipate via either grain crushing or plastic volumetric straining. We focus on these three parameters because they reveal important information about the underlying physics of the process. $\omega$ characterises whether grain breakage or pore collapse is the favoured form of energy dissipation in the system. $d_{Max}$ and $B$ tell us about the grain size distribution and how much energy is available to be dissipated in grain crushing. Both of these aspects of the system are particularly salient to our focus on the formation of cataclastic shear bands, making an analysis of the system sensitivity to changes in these parameters potentially illuminating for both our understanding of the physics, and the strengths and limitations of our model.

We first consider changes in the maximum grain size of our model. While these induce changes in the value of $\theta_{\gamma}$, thus changing the yield surfaces slightly, Figure 5 remains representative of which stress states support localisation.
<figure><img src="image_8.png" /><figcaption>Figure 8: The system behaviour at constant volume with increasing shear strain, for six different values of maximum grain size $d_{Max}$. (a) The shear stress $\tau_{21}$, and (b) the fastest growing wavelength $\Lambda_{max}$. The wavelengths have been normalised by the internal length scale of the respective material. Each of these curves were confined to the same initial $\gamma_{11}^e$ before shearing commenced, resulting in values of $R \approx 0.51$. A marker is placed on the minimum values of $\Lambda_{max}$ for clarity.</figcaption></figure>

In Figure 8 we observe that when we normalise the curves by their respective internal lengths $\ell^e$, we observe that the system effectively collapses onto one curve, despite orders of magnitude difference in the absolute size of the predicted shear band. It is important to note here that the normalising length is not constant here, and evolves due to the increase in B. We can see that in this instance the internal length sufficiently captures all of the relevant physics (both the absolute grain size and the relative polydispersity) to serve as an ideal normalisation length. In this case, the polydispersity is characterised purely by the ratio of the smallest to the largest particle size. However, we can also investigate a more complex measure of polydispersity by changing the relative weight of intermediate particle sizes, via changing the initial values of B. In this case, our plot in Figure 5 remains valid to characterise which regions of the yield surface can support shear band formation.
<figure><img src="image_9.png" /><figcaption>Figure 9: The system behaviour at constant volume with increasing shear strain, for four different initial values of B. (a) The shear stress $\tau_{21}$, (b) the fastest growing wavelength $\Lambda_{max}$ normalised by $d_{Max}$, (c) $\Lambda_{max}$ normalised by the respective initial values of $d_{50}$, (d) $\Lambda_{max}$ normalised by the respective values of $d_{50}$ as they evolve, (e) $\Lambda_{max}$ normalised by the respective initial values of $\ell^e$, and (f) $\Lambda_{max}$ normalised by the respective values of $\ell^e$ as they evolve. Each of these curves were confined to the same initial $\gamma_{11}$ before shearing commenced, resulting in values of R of 0.509, 0.362, 0.212 and 0.163 for B = 0, 0.3, 0.6 and 0.7 respectively. A marker is placed on the minimum values of $\Lambda_{max}$ for clarity.</figcaption></figure>

We plot in Figure 9 a number of different means of normalising the width of the predicted shear band. We also continue to shear up to $\gamma_{21} = 0.2$ to observe the longer-term tendencies of the behaviour. Firstly, we observe that the same qualitative behaviour is observed for all of the different initial values of B. However, as the initial particle breakage is increased, the formation of a shear band requires increasingly more shear strain, and the initial band width is wider, but also grows at a slower rate. This last effect is sufficiently pronounced that after $\gamma_{21} \approx 0.1$ the shear band associated with the initially unbroken material is predicted to be thicker relative to its grains than the bands with already broken grains.

In particular, when we study different normalisation strategies, we observe that normalising by the largest grain
size $d_{Max}$ gives the closest grouping of initial widths, and thus the most accurate indication of potential shear band widths, particularly at lower values of $B$, and thus polydispersity. We also note that when we normalise by the evolving internal length scale $\ell^e$, the curves converge faster than if we normalise solely by the respective initial values $\ell_0^e$. At larger values of shear, the normalisation with evolving length scale causes the three curves for $B \ge 0.3$ to converge. Most interestingly from these normalisation comparisons, we observe that normalising either by the initial or evolving mean grain size $d_{50}$, we obtain the largest spread between the initial localisation widths, and no or very slow convergence.

However, we note that the results suggesting that the shear band thickens with increasing $B$ in Figure 9(b) should be taken with caution. Due to the absence of porosity in the model, we do not capture the coupled effects of grain breakage leading to pore collapse, which may well lead to thinner bands with increasing $B$, as might be expected given the overall finer grain size distribution.

Now, by changing the material parameter $\omega$ we may instead examine the tendency of the system to favour dissipation via plastic volumetric straining over dissipation via grain crushing. In this instance the regions of the yield surface which can support shear band formation must be recalculated.

<figure><img src="image_10.png" /><figcaption>Figure 10: The regions of the yield surface that can in principle support localisation parallel to the shear direction when sheared under constant volume conditions, for (a) $\omega = 0°$, (b) $\omega = 45°$, and (c) $\omega = 90°$. The analysis is once again performed for $B = 0, 0.1, ..., 0.9$, with the correspondingly growing yield surfaces shown above. We place a marker to highlight the maximum confining stress capable of supporting localisation for each value of $B$.</figcaption></figure>

In Figure 10 we can see that changing the value of the parameter $\omega$ has a significant effect on the ability of the system to support shear band formation. In the system where $\omega = 0°$ (i.e. where grain crushing is most strongly favoured) only low confining pressures allow the formation of shear bands. The fraction of the yield surface that can support band formation also dramatically reduces as $B$ increases. By contrast, the system where $\omega = 90°$ (i.e. where pore collapse is most strongly favoured) supports shear band formation for any region of the yield surface to the left of the maximum value of $q$. This behaviour does not change with increasing values of $B$. The system where $\omega = 45°$ (i.e. where the two dissipative mechanisms are balanced) is intermediate between the two end members in both the extent of the bifurcation region at $B = 0$, and the behaviour of the bifurcation region as $B$ grows.
<figure><img src="image_11.png" /><figcaption>Figure 11: The system behaviour at constant volume with increasing shear strain, for four different values of $\omega$. (a) The shear stress $\tau_{21}$, and (b) the fastest growing wavelength $\Lambda_{max}$. Each of these curves were confined to the same initial $\gamma_{11}^e$ before shearing commenced, resulting in values of R of 0.617, 0.604, 0.557 and 0.491 for $\omega = 0, 20, 45$ and $90^\circ$ respectively. The wavelengths have been normalised by the largest grain size $d_{Max}$.</figcaption></figure>

In Figure 11 we can see the different behaviour of the system following a given load path, for various values of $\omega$. When considering the $\tau_{21} - \gamma_{21}$ plot we can see that for smaller $\omega$ the apparent softening is smaller. When we consider the localisation behaviour, we observe that increasing the value of $\omega$ results in localisation at a smaller value of shear strain, and the localisation is thinner relative to the constituent grains. The thickness of the band grows faster as $\omega$ increases.

# 7 Finite element method

Next, we use the finite element method (FEM) to more rigorously explore our system post-localisation. We consider the static momentum balance equations, and apply Dirichlet boundary conditions for both the displacements and rotations on the boundaries, corresponding to constant volume shearing.

The system is integrated in the Numerical Geolab framework (Stefanou, 2018), which provides a generalised continuum layer on top of the finite element solver FEniCS (Alnæs et al., 2015; Logg, Mardal and Wells, 2012). We use an implicit backwards Euler method, and select quadratic test functions for the displacements and linear test functions for the rotations (see Godio et al. (2015) for a more thorough discussion of finite element implementation in the Cosserat continuum).

When we vary the maximum grain size $d_{Max}$, or the initial breakage $B$, we observe that the FEM returns similar results to the LSA. However, when we vary the initial value of $\omega$, we observe some distinct behaviours not predicted by the LSA.
<figure><img src="image_12.png" /></figure>

Figure 12: The evolution of *B* for: (a) $\omega = 0^{\circ}$, with increments of homogeneous strain between the blue lines being $\Delta\gamma_{21}^h = 0.0076$, between the black lines $\Delta\gamma_{21}^h = 0.0036$, and between the red lines $\Delta\gamma_{21}^h = 0.1888$, (b) $\omega = 20^{\circ}$, with increments of homogeneous strain between the blue lines being $\Delta\gamma_{21}^h = 0.006$, between the black lines $\Delta\gamma_{21}^h = 0.0005$, and between the red lines $\Delta\gamma_{21}^h = 0.1935$, (c) $\omega = 45^{\circ}$, with increments of homogeneous strain between the blue lines being $\Delta\gamma_{21}^h = 0.0077$, between the black lines $\Delta\gamma_{21}^h = 0.0022$, and between the red lines $\Delta\gamma_{21}^h = 0.1901$. (d) The evolution of $\gamma_{21}^p$ for $\omega = 90^{\circ}$, with increments of homogeneous strain between the blue lines being $\Delta\gamma_{21}^h = 0.0329$.

In Figure 12, we observe that there is a distinct change in behaviour of the system for values of $\omega \le 45^{\circ}$. Rather than one stable shear band that slowly grows as predicted by the LSA (and as happens for our calibrated value $\omega = 70^{\circ}$ as well as $\omega = 90^{\circ}$), we observe instead three phases of behaviour. The first of these is the initial homogeneous deformation that transitions into a single shear band, as predicted. This is subsequently followed by a second phase where two distinct shear bands are growing at the edges of the initial shear band. Finally, one of these shear bands becomes dominant and continues to grow while shifting position in space. The dominant shear band ultimately saturates and reaches a state where grains can essentially no longer be broken.

# 8 Discussion

We observe certain trends across both our analyses using the calibrated values and the sensitivity analyses we have also conducted. Firstly, we observe that starting the test with different initial values of *B* does not have a very strong qualitative effect on the system behaviour, and that increasing initial values of *B* leads to increasing band thickness. However, we also observe that systems with less initial breakage are predicted to feature faster growing shear bands.

We also showed in Figure 9 that $d_{50}$ was not predictive of shear band width. This is a striking result, which can help to explain why researchers have thus far been unable to determine a constant scaling between $d_{50}$ and shear band width, with laboratory observations ranging from 1-2 $d_{50}$ (Bésuelle, Desrues and Raynaud, 2000; Ord, Vardoulakis and Kajewski, 1991) to 17 $d_{50}$ (Hall et al., 2010; Mühlhaus and Vardoulakis, 1987), as well as intermediate values (Andò et al., 2012; Roscoe, 1970). In the field, ultra-cataclastic zones have been observed that are thousands of grains thick (Chester and Chester, 1998). While it is important to note that while these results are obtained under different loading conditions or with additional physical processes involved, rendering a direct comparison with our predicted band widths moot, our results can help to explain the absence of a constant scaling with $d_{50}$. In essence, relying solely on $d_{50}$ fails to capture the complex role that polydispersity plays in shear band formation. When we change the polydispersity solely by varying the ratio between the smallest and largest grains, Figure 8 demonstrates that our
evolving internal length scale $l^e$ is sufficient to fully characterise the shear band width. However, when we introduce a more complex measure of polydispersity by changing the initial value of $B$, no single length that we are able to define is able to fully capture these effects. However, the largest grain size gives the best estimate for the initial band thickness in this case, particularly for less polydisperse systems.

We also see that the value of $\omega$ has a very strong effect on the overall behaviour. Firstly, when applying the LSA to the system as a whole to determine which stress conditions support localisation, we observe that systems which favour grain breakage over pore collapse only localise in regions of low confining stress. However, as we increase the tendency towards pore collapse by increasing $\omega$, we observe that the regions of the yield surface which can in principle support localisation also increase substantially, with the entire left-hand (i.e. shear dominated) side of the yield surface supporting localisation when $\omega = 90^\circ$, shown in Figure 10. This behaviour suggests that pore collapse is the driving physical mechanism behind localisation formation. When we analyse individual load paths, we observe that the analysis predicts diminishing initial thickness with increasing $\omega$. However, this is also accompanied by a more rapid increase in thickness post-localisation. The systems with higher values of $\omega$ also localise at smaller values of shear strain, as they demonstrate stronger apparent softening behaviour. Once again, these behaviours indicate that pore collapse is the key mechanism in localisation. Under the specified load conditions, the system is maximising the dissipation, first by the formation of a narrow band, which then grows as dissipation is maximised by further regions outside the band collapsing into it. This behaviour can be conceptualised as a localisation-delocalisation behaviour, with the further crushing of grains and collapse of pores becoming increasingly difficult in the original band, causing the zone to expand outwards. We observe qualitatively similar behaviour in faults, with the thickness of the fault being observed to increase as it is subjected to increasing shear (Shipton et al., 2006). This tendency has also been observed in experiments on dry sand subjected to arbitrarily large shear (Chambon, Schmittbuhl and Corfdir, 2006), so we can be confident that the observed behaviour is due to changes in the underlying granular skeleton, which our model accurately captures.

In addition, when we simulated our system using FEM, we observed that changing the value of $\omega$ introduced an unexpected behaviour, namely the formation of a double shear band structure.

<figure><img src="image_13.png" /><figcaption>Figure 13: An image of two ultra-cataclastic zones sandwiching a cataclastic zone, from Nicchio et al. (2018). The zone boundaries are highlighted in yellow.</figcaption></figure>

As Figure 13 shows, these double band structures occur in the field, in this case in feldspar-rich conglomerates. Similar observations can be found in Aydin and Johnson (1978), Bürgi, Parriaux and Franciosi (2001) and Wibberley and Shimamoto (2002) amongst others, in a variety of different rock types and locales. Indeed, in examining closely the figure we can see that not only do we observe a double cataclastic band either side of the central zone, but the right hand band appears to have experienced more extensive grain crushing than the left, suggesting that it may have become the dominant band and proceeded towards saturation, as our results suggest. The observation of these structures in the field can allow us to calibrate our models (specifically the value of $\omega$) and infer that at the time that the bands formed, grain crushing was the dominant mode of dissipation.

However, there are certainly improvements that can still be made to our model. The most important is to incorporate a true critical state, i.e. include porosity $\phi$ as an internal variable (Tengattini, Das and Einav, 2016). Introducing porosity to the model would increase the fidelity of the model towards the underlying physics of the system. In reality, the grain breakage and pore collapse processes are strongly coupled together, but in our model they are essentially competing processes, with the relative strength of each determined by $\omega$. Implementing an accurate pore-scale model is also essential in order to be able to expand the model to hydro-mechanical processes, where the collapse of pore space can lead to dramatic increases in pore pressure that in turn cause dramatic decreases in the effective stress
and overall sliding resistance of the system (Rattez, Stefanou and Sulem, 2018; Rattez et al., 2018a,b). We can also include aspects from models that account for cement bridges present in granular rocks (Das et al., 2014; Tengattini et al., 2014), or non-linear pressure-dependant elasticity (Nguyen and Einav, 2009; Tengattini, Das and Einav, 2016).

We also note that the Cosserat continuum cannot regularise problems featuring pure compaction or dilation bands. This could be addressed by moving to a full first-order micromorphic continuum (Stefanou, Sulem and Vardoulakis, 2010) i.e. one that includes the second gradient as well, or by using a viscous regularisation. Recent work has related rate-dependence at the continuum scale to crack growth within the individual grains (Marinelli and Buscarnera, 2019; Zhang and Buscarnera, 2017). A more refined constitutive model including one or both of these regularisation methods, combined with the additional physics discussed above, would seem to be the natural next step in model development. A model of this type would provide all the necessary ingredients to couple Cosserat (or first-order-micromorphic) Breakage Mechanics models to thermo-hydro-chemo-mechanical models so that the benefits of a full description of the GSD and its evolution can be brought to bear on fault modelling.

## 9 Conclusion

In this paper, we first expanded the upscaling procedures of Breakage Mechanics consistent with Einav (2007a) to the Cosserat continuum. This allows us to include the contributions of the elastic curvatures to the total stored energy, by making use of a new material constant $\theta_{\kappa}$ that can be measured without making any constitutive assumptions about the material. This material constant and the upscaling procedure allow us to define a new internal length for the system, $l^e$, that takes into account the entire grain size distribution, and evolves as that distribution evolves. The upscaling techniques were further extended to account for the moment of rotational inertia, allowing the definition of another new granulometry constant $\theta_I$ that allows the entire grain size distribution and its evolution to be taken into account when analysing the momentum balance equations of the Cosserat continuum. An analysis of the limit values of the new granulometry parameters was undertaken, which confirmed the necessity of including all of the parameters in the model in order to retain sufficient information about the grain size distribution.

We then outlined a constitutive model, expanding a previously developed model into the Cosserat continuum using our newly developed upscaled relationships. In doing so, we derived certain relationships between the second deviatoric stress and strain invariants that are slightly more general than in previous works, and consistent with the definition of the deviator stress used in triaxial tests. We also obtained certain relationships between the bulk and shear moduli of the classical continuum and three additional Cosserat stiffness moduli. We then presented methods to calibrate the values of all parameters. We emphasise here that the change to the Cosserat continuum requires no additional calibration over the same model expressed in the classical continuum, as all new material parameters can be calculated with the same information necessary to calculate their classical equivalents, based on typical geotechnical tests.

By applying linear stability analysis to the model under conditions corresponding to fast undrained shear in faults, we demonstrated that the polydispersity of the initial grain size distribution plays an important role in determining the initial thickness of the localisation with respect to its constituent grains. We also showed that the pore collapse process governs the range of stress states which can support localisation into a shear band, as well as the width of the initial localisation and the post-localisation tendency. Using the finite element method, we were able to show that systems in which grain breakage is dominant are able to form double-band structures.

This work has combined two successful families of models, adding additional predictive power to Breakage Mechanics and additional physical fidelity to Cosserat models. The development of this model has clarified the role of polydispersity in shear band formation, and helped to clarify the interaction between the micro-structure and the shear band width in crushable granular media.

## Acknowledgements

This work was supported by Australian Research Council Discovery Projects DP160104310 and DP190103487. The numerical developments (Numerical Geolab) were supported by the European Research Council (ERC) under the European Union Horizon 2020 research and innovation program (Grant agreement 757848 CoQuake). The authors acknowledge the Sydney Informatics Hub and the University of Sydney's high performance computing cluster Artemis which have contributed to the results reported in this paper.

The authors also thank Panos Papanastasiou, Leong Hien Poh and Giang Nguyen for their valuable comments that have helped to strengthen this work.
References

Alnæs, Martin S et al. (2015). "The FEniCS Project Version 1.5". In: *Archive of Numerical Software* 3.100, pp. 9-23.
DOI: https://doi.org/10.11588/ans.2015.100.20553.

Andò, Edward et al. (2012). "Grain-scale experimental investigation of localised deformation in sand: A discrete particle tracking approach". In: *Acta Geotechnica* 7.1, pp. 1-13. ISSN: 18611125. DOI: 10.1007/s11440-011-0151-6.

Aydin, Atilla and Arvid M. Johnson (1978). "Development of Faults as Zones of Deformation Bands and as Slip Surfaces in Sandstone". In: *Pure and Applied Geophysics* 116, pp. 931-942.

Baud, Patrick, Emmanuel Klein and Teng-fong Wong (2004). "Compaction localization in porous sandstones: spatial evolution of damage and acoustic emission activity". In: *Journal of Structural Geology* 26, pp. 603-624. DOI: 10.1016/j.jsg.2003.09.002.

Baud, Patrick, Veronika Vajdova and Teng-fong Wong (2006). "Shear-enhanced compaction and strain localization: Inelastic deformation and constitutive modeling of four porous sandstones". In: *Journal of Geophysical Research: Solid Earth* 111.12, pp. 1-17. ISSN: 21699356. DOI: 10.1029/2005JB004101.

Bauer, Erich (2016). "Simulation of the Influence of Grain Damage on the Evolution of Shear Strain Localization". In: *Continuous Media with Microstructure* 2. Ed. by B. Albers and M. Kuczma. Springer International Publishing, pp. 231-244. ISBN: 9783319282411. DOI: 10.1007/978-3-319-28241-1.

Bauer, Erich, Linke Li and Wenxiong Huang (2015). "Hypoplastic Constitutive Modelling of Grain Damage Under Plane Shearing". In: *Bifurcation and Degradation of Geomaterials in the New Millenium: Proceedings of the 10th International Workshop on Bifurcation and Degradation in Geomaterials*. Ed. by Kam-Tim Chau and Jidong Zhao, pp. 181-187. ISBN: 9783319135052. DOI: 10.1007/978-3-319-13506-9.

Bauer, Erich, Linke Li and Mohammedkeya Khosravi (2017). "Modelling Grain Damage Under Plane Strain Compression Using a Micro-polar Continuum". In: *Bifurcation and Degradation of Geomaterials with Engineering Applications: Proceedings of the 11th International Workshop on Bifurcation and Degradation in Geomaterials dedicated to Hans Muhlhaus, Limassol, Cyprus, 21-25 May 2017*. Ed. by Euripides Papamichos et al. Springer International Publishing, pp. 349-359. ISBN: 9783319563961. DOI: 10.1007/978-3-319-56397-8. URL: http://dx.doi.org/10.1007/978-3-319-56397-8.

Bésuelle, Pierre, Jacques Desrues and Suzanne Raynaud (2000). "Experimental characterisation of the localisation phenomenon inside a Vosges sandstone in a triaxial cell". In: *International Journal of Rock Mechanics and Mining Sciences* 37, pp. 1223-1237. ISSN: 1365-1609. DOI: 10.1016/S1365-1609(00)00057-5.

Brantut, Nicolas and Jean Sulem (2012). "Strain Localization and Slip Instability in a Strain-Rate Hardening, Chemically Weakening Material". In: *Journal of Applied Mechanics* 79.3, p. 031004. ISSN: 00218936. DOI: 10.1115/1.4005880.

Brantut, Nicolas, Jean Sulem and A. Schubnel (2011). "Effect of dehydration reactions on earthquake nucleation: Stable sliding, slow transients, and unstable slip". In: *Journal of Geophysical Research: Solid Earth* 116.5, pp. 1-16. ISSN: 21699356. DOI: 10.1029/2010JB007876.

Bürgi, C., A. Parriaux and G. Franciosi (2001). "Geological characterization of weak cataclastic fault rocks with regards to the assessment of their geomechanical properties". In: *Quarterly Journal of Engineering Geology and Hydrogeology* 34.2, pp. 225-232. ISSN: 14709236. DOI: 10.1144/qjegh.34.2.225.

Buscarnera, Giuseppe and Arghya Das (2016). "Chemo-mechanics of cemented granular solids subjected to precipitation and dissolution of mineral species". In: *International Journal for Numerical and Analytical Methods in Geomechanics* 40, pp. 1295-1320. ISSN: 03639061. DOI: 10.1002/nag.2486.

Buscarnera, Giuseppe and Itai Einav (2012). "The yielding of brittle unsaturated granular soils". In: *Géotechnique* 62.2, pp. 147-160. ISSN: 0016-8505. DOI: 10.1680/geot.10.P.118.

Chambon, Guillaume, Jean Schmittbuhl and Alain Corfdir (2006). "Frictional response of a thick gouge sample: 1. mechanical measurements and microstructures". In: *Journal of Geophysical Research: Solid Earth* 111.9, pp. 1-17. ISSN: 21699356. DOI: 10.1029/2003JB002731.

Chester, Frederick M. and Judith S. Chester (1998). "Ultracataclasite structure and friction processes of the Punchbowl fault, San Andreas system, California". In: *Tectonophysics* 295.1-2, pp. 199-221. ISSN: 00401951. DOI: 10.1016/S0040-1951(98)00121-8.

Cosserat, François and Eugène Cosserat (1909). *Théorie des corps déformables*. Paris: Librairie Scientifique A. Hermann et Fils.

Das, Arghya (2013). "A theoretical study of grain crushing induced compaction localization in porous sandstones". PhD thesis. The University of Sydney. ISBN: 9783642328145.

Das, Arghya, Giang D. Nguyen and Itai Einav (2011). "Compaction bands due to grain crushing in porous rocks: A theoretical approach based on breakage mechanics". In: *Journal of Geophysical Research: Solid Earth* 116.8, pp. 1-14. ISSN: 21699356. DOI: 10.1029/2011JB008265.

— (2013). "The propagation of compaction bands in porous rocks based on breakage mechanics". In: *Journal of Geophysical Research: Solid Earth* 118.5, pp. 2049-2066. ISSN: 21699356. DOI: 10.1002/jgrb.50193.

Das, Arghya et al. (2014). "A thermomechanical constitutive model for cemented granular materials with quantifiable internal variables. Part II - Validation and localization analysis". In: *Journal of the Mechanics and Physics of*
Solids 70.1, pp. 382–405. ISSN: 00225096. DOI: 10.1016/j.jmps.2014.05.022. URL: http://dx.doi.org/10.1016/j.jmps.2014.05.022.

Einav, Itai (2007a). "Breakage mechanics-Part I: Theory". In: *Journal of the Mechanics and Physics of Solids* 55.6, pp. 1274–1297. ISSN: 00225096. DOI: 10.1016/j.jmps.2006.11.003.

— (2007b). "Breakage mechanics-Part II: Modelling granular materials". In: *Journal of the Mechanics and Physics of Solids* 55.6, pp. 1298–1320. ISSN: 00225096. DOI: 10.1016/j.jmps.2006.11.004.

— (2007c). "Fracture propagation in brittle granular matter". In: *Proceedings of the Royal Society A* 463.August, pp. 3021–3035. ISSN: 1364-5021. DOI: 10.1098/rspa.2007.1898.

— (2007d). "Soil mechanics : breaking ground". In: *Philosophical Transactions of The Royal Society A* 365.1861, pp. 2985–3002. ISSN: 1364-503X. DOI: 10.1098/rsta.2007.0009.

— (2012). "The unification of hypo-plastic and elasto-plastic theories". In: *International Journal of Solids and Structures* 49.11-12, pp. 1305–1315. ISSN: 00207683. DOI: 10.1016/j.ijsolstr.2012.02.003. URL: http://dx.doi.org/10.1016/j.ijsolstr.2012.02.003.

Einav, Itai and Giang D. Nguyen (2008). "Hardening plasticity for geomaterials : micromechanical roots and alternative formulations". In: *5th. European Congress on Computational Methods in Applied Sciences and Engineering*.

Einav, Itai and Julio R. Valdes (2008). "On comminution and yield in brittle granular mixtures". In: *Journal of the Mechanics and Physics of Solids* 56.6, pp. 2136–2148. ISSN: 00225096. DOI: 10.1016/j.jmps.2008.02.002.

Eringen, A. Cemal (1967). "Linear theory of micropolar viscoelasticity". In: *International Journal of Engineering Science* 5.2, pp. 191–204. ISSN: 00207225. DOI: 10.1016/0020-7225(67)90004-3.

Germain, P (1973). "The Method of Virtual Power in Continuum Mechanics - Part II: Microstructure". In: *SIAM Journal on Applied Mathematics* 25.3, pp. 556–575.

Godio, Michele et al. (2015). "Dynamic finite element formulation for Cosserat elastic plates". In: *International Journal for Numerical Methods in Engineering* 101, pp. 992–1018. ISSN: 17359260. DOI: 10.1002/nme.

Hall, Stephen A. et al. (2010). "Discrete and continuum analysis of localised deformation in sand using X-ray µCT and volumetric digital image correlation". In: *Géotechnique* 60.5, pp. 315–322. ISSN: 0016-8505. DOI: 10.1680/geot.2010.60.5.315. URL: http://www.icevirtuallibrary.com/doi/10.1680/geot.2010.60.5.315.

Herbold, E. B., M. A. Homel and Miles B. Rubin (2020). "A thermomechanical breakage model for shock-loaded granular media". In: *Journal of the Mechanics and Physics of Solids* 137, pp. 1–21. ISSN: 00225096. DOI: 10.1016/j.jmps.2019.103813.

Houlsby, Guy T. and Alexander M. Puzrin (2007). *Principles of Hyperplasticity*. London: Springer-Verlag London, p. 351. ISBN: 978-1-84628-240-9. DOI: 10.1007/978-1-84628-240-9. arXiv: arXiv:1011.1669v3. URL: http://link.springer.com/10.1007/978-1-84628-240-9.

Kendall, K. (1978). "The impossibility of comminuting small particles by compression". In: *Nature* 272, pp. 710–711.

Logg, Anders, Kent-Andre Mardal and Garth Wells, eds. (2012). *Automated Solution of Differential Equations by the Finite Element Method*. Vol. 84. Lecture Notes in Computational Science and Engineering. Berlin, Heidelberg: Springer Berlin Heidelberg. ISBN: 978-3-642-23099-8. DOI: 10.1007/978-3-642-23099-8. URL: http://link.springer.com/10.1007/978-3-642-23099-8.

Marinelli, Ferdinando and Giuseppe Buscarnera (2019). "A Generalized Backward Euler algorithm for the numerical integration of a viscous breakage model". In: *International Journal for Numerical and Analytical Methods in Geomechanics* 43.1, pp. 3–29. ISSN: 10969853. DOI: 10.1002/nag.2841.

Mühlhaus, Hans-Bernd and Ioannis Vardoulakis (1987). "The thickness of shear bands in granular materials". In: *Géotechnique* 37.3, pp. 271–283.

Nguyen, Giang D. (2011). "A damage model with evolving nonlocal interactions". In: *International Journal of Solids and Structures* 48.10, pp. 1544–1559. ISSN: 00207683. DOI: 10.1016/j.ijsolstr.2011.02.002. arXiv: j.ijsolstr.2011.02.002 [10.1016/j.ijsolstr.2011.02.002].

Nguyen, Giang D. and Itai Einav (2009). "The Energetics of Cataclasis Based on Breakage Mechanics". In: *Pure and Applied Geophysics* 166, pp. 1693–1724. DOI: 10.1007/978-3-0346-0138-2. URL: http://link.springer.com/10.1007/978-3-0346-0138-2.

— (2010). "Nonlocal regularisation of a model based on breakage mechanics for granular materials". In: *International Journal of Solids and Structures* 47.10, pp. 1350–1360. ISSN: 0020-7683. DOI: 10.1016/j.ijsolstr.2010.01.020. URL: http://dx.doi.org/10.1016/j.ijsolstr.2010.01.020.

Nicchio, Matheus A. et al. (2018). "Development of cataclastic foliation in deformation bands in feldspar-rich conglomerates of the Rio do Peixe Basin, NE Brazil". In: *Journal of Structural Geology* 107.July 2017, pp. 132–141. ISSN: 01918141. DOI: 10.1016/j.jsg.2017.12.013. URL: https://doi.org/10.1016/j.jsg.2017.12.013.

Ord, Alison, Ioannis Vardoulakis and R. Kajewski (1991). "Shear band formation in Gosford Sandstone". In: *International Journal of Rock Mechanics and Mining Sciences* 28.5, pp. 397–409. ISSN: 01489062. DOI: 10.1016/0148-9062(91)90078-Z.

Papanastasiou, Panos and Ioannis Vardoulakis (1989). "Bifurcation analysis of deep boreholes: II. Scale effect". In: *International Journal for Numerical and Analytical Methods in Geomechanics* 13.2, pp. 183–198. ISSN: 01489062. DOI: 10.1016/0148-9062(89)90982-0.
Papanastasiou, Panos and Ioannis Vardoulakis (1992). "Numerical treatment of progressive localization in relation to borehole stability". In: *International Journal for Numerical and Analytical Methods in Geomechanics* 16.6, pp. 389-424. ISSN: 10969853. DOI: 10.1002/nag.1610160602.

Pijaudier-Cabot, Gilles, Khalil Haidar and Jean François Dubé (2004). "Non-local damage model with evolving internal length". In: *International Journal for Numerical and Analytical Methods in Geomechanics* 28.7-8, pp. 633-652. ISSN: 03639061. DOI: 10.1002/nag.367.

Platt, John D., John W. Rudnicki and James R. Rice (2014). "Stability and localization of rapid shear in fluid-saturated fault gouge: 2. Localized zone width and strength evolution". In: *Journal of Geophysical Research: Solid Earth* 119.5, pp. 4334-4359. ISSN: 21699356. DOI: 10.1002/2013JB010711.

Poh, Leong Hien and Gang Sun (2017). "Localizing gradient damage model with decreasing interactions". In: *International Journal for Numerical Methods in Engineering* 110.6, pp. 503-522. ISSN: 10970207. DOI: 10.1002/nme.5364.

Rastiello, Giuseppe et al. (2018). "From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length". In: *Computer methods in applied mechanics and engineering* 331, pp. 650-674. DOI: 10.1016/j.cma.2017.12.006.

Rattez, Hadrien (2018). "Couplages Thermo-Hydro-Mécaniques et localisation dans les milieux de Cosserat : Application à l'analyse de stabilité du cisaillement rapide des failles". PhD thesis. Université Paris-Est.

Rattez, Hadrien, Ioannis Stefanou and Jean Sulem (2018). "The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part I: Theory and linear stability analysis". In: *Journal of the Mechanics and Physics of Solids* 115, pp. 54-76. ISSN: 00225096. DOI: 10.1016/j.jmps.2018.03.004. URL: https://doi.org/10.1016/j.jmps.2018.03.004.

Rattez, Hadrien et al. (2018a). "Numerical Analysis of Strain Localization in Rocks with Thermo-hydro-mechanical Couplings Using Cosserat Continuum". In: *Rock Mechanics and Rock Engineering* 51.10, pp. 3295-3311. ISSN: 07232632. DOI: 10.1007/s00603-018-1529-7. URL: http://dx.doi.org/10.1007/s00603-018-1529-7.

— (2018b). "The importance of Thermo-Hydro-Mechanical couplings and microstructure to strain localization in 3D continua with application to seismic faults. Part II: Numerical implementation and post-bifurcation analysis". In: *Journal of the Mechanics and Physics of Solids* 115, pp. 54-76. ISSN: 00225096. DOI: 10.1016/j.jmps.2018.03.004. URL: https://doi.org/10.1016/j.jmps.2018.03.003.

Rice, James R. (2006). "Heating and weakening of faults during earthquake slip". In: *Journal of Geophysical Research: Solid Earth* 111.5, pp. 1-29. ISSN: 21699356. DOI: 10.1029/2005JB004006.

Rice, James R., John W. Rudnicki and John D. Platt (2014). "Stability and localization of rapid shear in fluid-saturated fault gouge: 1. Linearized stability analysis". In: *Journal of Geophysical Research: Solid Earth* 119.5, pp. 4311-4333. ISSN: 21699356. DOI: 10.1002/2013JB010710.

Roscoe, K. H. (1970). "The Influence of Strains in Soil Mechanics". In: *Géotechnique* 20.2, pp. 129-170. ISSN: 0016-8505. DOI: 10.1680/geot.1970.20.2.129. URL: http://www.icevirtuallibrary.com/doi/10.1680/geot.1970.20.2.129.

Rubin, Miles B. (2001). "Physical reasons for abandoning plastic deformation measures in plasticity and viscoplasticity theory". In: *Archives of Mechanics* 53.4-5, pp. 519-539.

Rudnicki, John W. and James R. Rice (1975). "Conditions for the localization of deformation in pressure-sensitive dilatant materials". In: *Journal of the Mechanics and Physics of Solids* 23.6, pp. 371-394. ISSN: 00225096. DOI: 10.1016/0022-5096(75)90001-0.

Sammis, Charles, Geoffrey King and Ronald Biegel (1987). "The kinematics of gouge deformation". In: *Pure and Applied Geophysics* PAGEOPH 125.5, pp. 777-812. ISSN: 00334553. DOI: 10.1007/BF00878033.

Scherer, Jean-Michel et al. (2019). "Strain gradient crystal plasticity with evolving length scale: Application to voided irradiated materials". In: *European Journal of Mechanics - A/Solids*. ISSN: 09977538. DOI: 10.1016/j.euromechsol.2019.04.003. URL: https://linkinghub.elsevier.com/retrieve/pii/S0997753818309458.

Shipton, Zoe K. et al. (2006). "How Thick is a Fault? Fault Displacement-Thickness Scaling Revisited". In: *Earthquakes: Radiated Energy and the Physics of Faulting*. American Geophysical Union, pp. 193-198. ISBN: 978-0-87590-435-1. DOI: 10.1029/170GM19.

Sibson, Richard H. (2003). "Thickness of the seismic slip zone". In: *Bulletin of the Seismological Society of America* 93.3, pp. 1169-1178. ISSN: 00371106. DOI: 10.1785/0120020061.

Stefanou, Ioannis (2018). "Numerical Geolab, Controlling earthQuakes project. Nantes. URL: www.coquake.eu."

— (2019). "Controlling Anthropogenic and Natural Seismicity: Insights From Active Stabilization of the Spring-Slider Model". In: *Journal of Geophysical Research: Solid Earth* 124.8, pp. 8786-8802. ISSN: 21699356. DOI: 10.1029/2019JB017847. URL: https://doi.org/10.1029/2019JB017847.

Stefanou, Ioannis and Jean Sulem (2013). "Formation de bandes de compaction dans les géomatériaux causée par une dégradation chimique". In: *21ème Congrès Français de Mécanique*.

— (2014). "Chemically induced compaction bands: Triggering conditions and band thickness". In: *Journal of Geophysical Research: Solid Earth* 119.100, pp. 8964-8975. ISSN: 21699313. DOI: 10.1002/2014JB011567. Received.
Stefanou, Ioannis, Jean Sulem and Ioannis Vardoulakis (2008). "Three-dimensional Cosserat homogenization of masonry structures: Elasticity". In: *Acta Geotechnica* 3.1, pp. 71–83. ISSN: 18611125. DOI: 10.1007/s11440-007-0051-y.

— (2010). "Homogenization of interlocking masonry structures using a generalized differential expansion technique". In: *International Journal of Solids and Structures* 47.11-12, pp. 1522–1536. ISSN: 00207683. DOI: 10.1016/j.ijsolstr.2010.02.011.

Sulem, Jean (2007). "Stress orientation evaluated from strain localisation analysis in Aigion Fault". In: *Tectonophysics* 442.1-4, pp. 3–13. ISSN: 00401951. DOI: 10.1016/j.tecto.2007.03.005. arXiv: 0808.2012.

— (2010). "Bifurcation theory and localization phenomena". In: *European Journal of Environmental and Civil Engineering* 14.8-9, pp. 989–1009. ISSN: 1964-8189. DOI: 10.1080/19648189.2010.9693276.

Sulem, Jean and V Famin (2009). "Thermal decomposition of carbonates in fault zones: Slip-weakening and temperature-limiting effects". In: *Journal of Geophysical Research* 114.B6, pp. 1–14. ISSN: 0148-0227. DOI: 10.1029/2009JB006576.

Sulem, Jean and Ioannis Stefanou (2016). "Thermal and chemical effects in shear and compaction bands". In: *Geo-mechanics for Energy and the Environment* 6, pp. 4–21. ISSN: 2352-3808. DOI: 10.1016/j.gete.2015.12.004. URL: http://dx.doi.org/10.1016/j.gete.2015.12.004.

Sulem, Jean, Ioannis Stefanou and Emmanuil Veveakis (2011). "Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure". In: *Granular Matter* 13.3, pp. 261–268. ISSN: 14345021. DOI: 10.1007/s10035-010-0244-1.

Sulem, Jean and Ioannis Vardoulakis (1990). "Bifurcation analysis of the triaxial test on rock specimens. A theoretical model for shape and size effect". In: *Acta Mechanica* 83.3-4, pp. 195–212. ISSN: 00015970. DOI: 10.1007/BF01172981.

Tembe, Sheryl, Patrick Baud and Teng-fong Wong (2008). "Stress conditions for the propagation of discrete compaction bands in porous sandstone". In: *Journal of Geophysical Research: Solid Earth* 113.9, pp. 1–16. ISSN: 21699356. DOI: 10.1029/2007JB005439.

Tengattini, Alessandro (2015). "A micro-mechanical study of cemented granular materials". PhD thesis. The University of Sydney.

Tengattini, Alessandro, Arghya Das and Itai Einav (2014). "A theory predicting breakage dependence of critical state in sand". In: *Geomechanics from Micro to Macro*. Vol. 1. 1. London: Taylor & Francis Group, pp. 695–698. ISBN: 9781138027077. DOI: 10.1201/b17395-124. URL: http://www.crcnetbase.com/doi/abs/10.1201/b17395-124.

— (2016). "A constitutive modelling framework predicting critical state in sand undergoing crushing and dilation". In: *Géotechnique* 66.9, pp. 695–710. ISSN: 0016-8505. DOI: 10.1680/jgeot.14.P.164. URL: http://www.icevirtuallibrary.com/doi/10.1680/jgeot.14.P.164.

Tengattini, Alessandro et al. (2014). "A thermomechanical constitutive model for cemented granular materials with quantifiable internal variables. Part I—Theory". In: *Journal of the Mechanics and Physics of Solids* 70.1, pp. 281–296. ISSN: 00225096. DOI: 10.1016/j.jmps.2014.05.021. URL: http://linkinghub.elsevier.com/retrieve/pii/S0022509614001185.

Turcotte, D. L. (1986). "Fractals and Fragmentation". In: *Journal of Geophysical Research* 91.5, pp. 1921–1926.

Vardoulakis, Ioannis (2019). *Cosserat Continuum Mechanics*. Ed. by Eleni Gerolymatou et al. Cham: Springer International Publishing. ISBN: 978-3-319-95155-3. DOI: 10.1007/978-3-319-95156-0. URL: http://link.springer.com/10.1007/978-3-319-95156-0.

Vardoulakis, Ioannis and Panos Papanastasiou (1988). "Bifurcation analysis of deep boreholes: I. Surface instabilities". In: *International Journal for Numerical and Analytical Methods in Geomechanics* 12.4, pp. 379–399. ISSN: 10969853. DOI: 10.1002/nag.1610120404.

Vardoulakis, Ioannis and Jean Sulem (1995). *Bifurcation analysis in geomechanics*. Glasgow: Blackie Academic & Professional, p. 466. ISBN: 0203697030. DOI: 10.1061/(ASCE)0733-9399(2006)132:9(952).

Veveakis, Emmanuil, Sotiris Alevizos and Ioannis Vardoulakis (2010). "Chemical reaction capping of thermal instabilities during shear of frictional faults". In: *Journal of the Mechanics and Physics of Solids* 58.9, pp. 1175–1194. ISSN: 00225096. DOI: 10.1016/j.jmps.2010.06.010.

Veveakis, Emmanuil, Ioannis Stefanou and Jean Sulem (2013). "Failure in shear bands for granular materials: thermo-hydro-chemomechanical effects". In: *Géotechnique Letters* 3.2, pp. 31–26.

Veveakis, Emmanuil, Ioannis Vardoulakis and Giulio Di Toro (2007). "Thermoporomechanics of creeping landslides: The 1963 Vaiont slide, northern Italy". In: *Journal of Geophysical Research* 112.April, pp. 1–21. DOI: 10.1029/2006JF000702.

Viswanath, Parol and Arghya Das (2019). "Modelling the response of chemically degraded carbonate sands". In: *Geomechanics and Geoengineering* 6025. ISSN: 17486033. DOI: 10.1080/17486025.2019.1587178.

Wang, Zhao and Leong Hien Poh (2018). "A homogenized localizing gradient damage model with micro inertia effect". In: *Journal of the Mechanics and Physics of Solids* 116, pp. 370–390. ISSN: 00225096. DOI: 10.1016/j.jmps.2018.04.007. URL: https://doi.org/10.1016/j.jmps.2018.04.007.
Wibberley, Christopher A J and Toshihiko Shimamoto (2002). "Internal structure and permeability of major strike-slip fault zones: The Median Tectonic Line in Mie Prefecture, Southwest Japan". In: *Journal of Structural Geology* 25.1, pp. 59-78. ISSN: 01918141. DOI: 10.1016/S0191-8141(02)00014-7.

Wong, Teng-fong, Patrick Baud and Emmanuelle Klein (2001). "Localized failure modes in a compactant porous rock". In: *Geophysical Research Letters* 28.13, pp. 2521-2524.

Zervos, Antonios, Panos Papanastasiou and Ioannis Vardoulakis (2001). "Modelling of localisation and scale effect in thick-walled cylinders with gradient elastoplasticity". In: *International Journal of Solids and Structures* 38.30-31, pp. 5081-5095. ISSN: 00207683. DOI: 10.1016/S0020-7683(00)00337-1.

Zhang, Chunshun (2012). "Breakage mechanics in large granular flow problems". PhD thesis. The University of Sydney, pp. 1-158.

Zhang, Chunshun, Itai Einav and Giang D. Nguyen (2012). "The end-bearing capacity of piles penetrating into crushable soils". In: *Géotechnique* 63.5, pp. 1-14. ISSN: 0016-8505. DOI: 10.1680/geot.11.P.117.

Zhang, Yida and Giuseppe Buscarnera (2014). "Grainsize dependence of clastic yielding in unsaturated granular soils". In: *Granular Matter* 16, pp. 469-483. DOI: 10.1007/s10035-014-0491-7.

— (2015). "Prediction of breakage-induced couplings in unsaturated granular soils". In: *Géotechnique* 65.2, pp. 135-140.

— (2017). "A rate-dependent breakage model based on the kinetics of crack growth at the grain scale". In: *Géotechnique* 67.11, pp. 953-967. ISSN: 0016-8505. DOI: 10.1680/jgeot.16.p.181.

— (2018). "Breakage mechanics for granular materials in surface-reactive environments". In: *Journal of the Mechanics and Physics of Solids* 112, pp. 89-108. ISSN: 00225096. DOI: 10.1016/j.jmps.2017.11.008. URL: https://doi.org/10.1016/j.jmps.2017.11.008.

Zhang, Yida, Giuseppe Buscarnera and Itai Einav (2016). "Grain size dependence of yielding in granular soils interpreted using fracture mechanics, breakage mechanics and Weibull statistics". In: *Géotechnique* 66.2, pp. 149-160. ISSN: 0016-8505. DOI: 10.1680/jgeot.15.P.119. URL: http://www.icevirtuallibrary.com/doi/10.1680/jgeot.15.P.119.