Micro-mechanics of fabric and failure in granular materials

TL;DR

This paper develops a micro-mechanical theory for fabric evolution and failure in granular materials, validated by DEM simulations, to predict how anisotropy affects soil strength at large strains.

Contribution

It introduces a mathematical framework linking fabric anisotropy to micro-mechanics, verified through DEM simulations, for predicting granular material behavior under loading.

Findings

Fabric anisotropy evolves via convection, contact generation, and diffusion.

DEM simulations confirm the proposed fabric evolution model.

The theory predicts the influence of intermediate principal stress on strength.

Abstract

The paper addresses the underlying source of two forms of induced anisotropy in granular materials: contact orientation anisotropy and contact force anisotropy. A rational, mathematical structure is reviewed for the manner in which fabric anisotropy emerges and evolves during loading. Fabric is expressed as an orientation density, and transport phenomena such as convection, contact generation, and diffusion control the rate of fabric evolution during loading. The paper proposes specific measurable forms for all terms, based upon the micro-mechanics of particle interactions. Discrete element (DEM) simulations are used to verify and quantify these terms, so that the theory can be applied to general loading conditions. The DEM simulations are of densely packed durable spheres, and the emphasis is on soil behavior at large strains, specifically on fabric and strength at the critical state. Once the theory has been developed and quantified, it is applied to predict the effect of the intermediate principal stress on strength.

Tables

Table 1. Table 1: Characteristics of DEM assemblies

Characteristic	Value
Assemblies	20
Assembly shape	cube
Assembly particles	4096
Assembly dimension	13.4$D_{50}$
Assembly boundaries	periodic
Particle shape	spherical
Particle size range	0.4$D_{50}$ – 1.2$D_{50}$
Particle shear modulus, $G$	29 GPa
Particle Poisson ratio, $\nu$	0.15
Inter-particle friction ratio	0.50
Initial particle arrangement	dense, isotropic
Initial void ratio, solids fraction	0.510, 0.662
Initial avg. pressure, $p_{\text{o}}=-{\sigma }_{ii}/3$	320 kPa
Initial avg. coord. no., $2M/N$	5.48

Equations

$$M/N={\int }_{\Omega }\hat{g}(\mathbf{n})d\Omega =4\pi \widehat{g_{\text{avg}}}$$

$$\overline{\mathbf{f}}(\mathbf{n})=\hat{\mathbf{f}}(\mathbf{n})/\hat{g}(\mathbf{n})$$

$$\mathbf{\sigma }=\frac{1}{V}\sum _{m=1}^M{\mathbf{l}}^m\otimes {\mathbf{f}}^m$$

$$\mathbf{\sigma }=\frac{1}{V}\sum _{m=1}^M{\ell }^m{\mathbf{n}}^m\otimes {\mathbf{f}}^m\approx \frac{\overline{\ell }}{V}\sum _{m=1}^M{\mathbf{n}}^m\otimes {\mathbf{f}}^m$$

$$\mathbf{\sigma }\approx \frac{\overline{\ell }N}{V}{\int }_{\Omega }\mathbf{n}\otimes \hat{\mathbf{f}}(\mathbf{n})d\Omega =\frac{\overline{\ell }}{\overline{v}}{\int }_{\Omega }\mathbf{n}\otimes \hat{\mathbf{f}}(\mathbf{n})d\Omega$$

$$\dot{\mathbf{\sigma }}=-\frac{\dot{v}}{\overline{v}}\mathbf{\sigma }+\frac{\overline{\ell }}{\overline{v}}{\int }_{\Omega }\mathbf{n}\otimes {\frac{\partial }{\partial t}\hat{\mathbf{f}}(\mathbf{n})|}_{\mathbf{n}}d\Omega$$

$${\frac{\partial \hat{\mathbf{f}}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}=\frac{d\hat{\mathbf{f}}(\mathbf{n})}{dt}-{\nabla }_{{}_{\Omega }}\cdot \left(\dot{\mathbf{n}}\otimes \hat{\mathbf{f}}(\mathbf{n})\right)+{\left(\frac{\partial \hat{\mathbf{f}}(\mathbf{n})}{\partial t}\right)}_{\mathrm{diff}}$$

$${\frac{\partial \hat{g}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}=\frac{d\hat{g}(\mathbf{n})}{dt}-{\nabla }_{{}_{\Omega }}\cdot \left(\dot{\mathbf{n}}\hat{g}(\mathbf{n})\right)+{\left(\frac{\partial \hat{g}(\mathbf{n})}{\partial t}\right)}_{\mathrm{diff}}$$

$$\begin{array}{rl}d{\mathbf{f}}^m=& {\left(d{\mathbf{f}}^m\right)}_{\text{matl}}+{\mathbf{f}}^m\times (d{\mathbf{n}}^m\times {\mathbf{n}}^m)\\ & +\frac{1}{2}\left[\left(d{\mathbf{\theta }}^i+d{\mathbf{\theta }}^j\right)\cdot {\mathbf{n}}^m\right]\left({\mathbf{n}}^m\times {\mathbf{f}}^m\right)\end{array}$$

$$\hat{\mathbf{f}}(\mathbf{n})=-\widehat{f^{\text{n}}}(\mathbf{n})\mathbf{n}+\widehat{f^{\text{t}}}(\mathbf{n})\mathbf{t}(\mathbf{n})$$

$$\frac{d\hat{\mathbf{f}}(\mathbf{n})}{dt}=\frac{d\widehat{{\mathbf{f}}^{\text{n}}}(\mathbf{n})}{dt}+\frac{d\widehat{{\mathbf{f}}^{\text{t}}}(\mathbf{n})}{dt}$$

$$\begin{array}{rl}\frac{d\widehat{{\mathbf{f}}^{\text{n}}}(\mathbf{n})}{dt}=& -{\left(\frac{\partial \widehat{f^{\text{n}}}(\mathbf{n})}{\partial t}\right)}_{\text{matl}}\mathbf{n}\\ & -\widehat{f^{\text{n}}}(\mathbf{n})\mathbf{n}\times \left(\dot{\mathbf{n}}(\mathbf{n})\times \mathbf{n}\right)\end{array}$$

$$\begin{gathered} \frac{d\widehat{{\mathbf{f}}^{\text{t}}}(\mathbf{n})}{dt}={\left(\frac{\partial \widehat{f^{\text{t}}}(\mathbf{n})}{\partial t}\right)}_{\text{matl}}\mathbf{s}(\mathbf{n}) \\ +\widehat{f^{\text{t}}}(\mathbf{n})\left[\mathbf{t}(\mathbf{n})\times \left(\dot{\mathbf{n}}(\mathbf{n})\times \mathbf{n}\right)\right]+{\left(\frac{d\widehat{{\mathbf{f}}^{\text{t}}}\left(\mathbf{n}\right)}{dt}\right)}_{\text{twirl}} \end{gathered}$$

$$\begin{array}{rl}{\nabla }_{{}_{\Omega }}\cdot \left(\dot{\mathbf{n}}\otimes \hat{\mathbf{f}}(\mathbf{n})\right)=& -{\dot{n}}_{j,j}\widehat{f^{\text{n}}}{n}_i-{\dot{n}}_j{\widehat{f^{\text{n}}}}_{,j}{n}_i-{\dot{n}}_j\widehat{f^{\text{n}}}{n}_{i,j}\\ & +{\dot{n}}_{j,j}\widehat{f^{\text{t}}}{t}_i+{\dot{n}}_j{\widehat{f^{\text{t}}}}_{,j}{t}_i+{\dot{n}}_j\widehat{f^{\text{t}}}{t}_{i,j}\end{array}$$

$$-{\dot{n}}_j\widehat{f^{\text{n}}}{n}_{i,j}=$$

$${\dot{n}}_j\widehat{f^{\text{t}}}{t}_{i,j}=$$

$$\widehat{f^{\text{t}}}(\mathbf{n})\dot{\mathbf{t}}(\mathbf{n})=-\widehat{f^{\text{t}}}(\mathbf{n})\left[\mathbf{t}\cdot \dot{\mathbf{n}}(\mathbf{n})\right]\mathbf{n}+{\left(\frac{d\widehat{{\mathbf{f}}^{\text{t}}}\left(\mathbf{n}\right)}{dt}\right)}_{\text{veer}}$$

$${\frac{\partial \hat{\mathbf{f}}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}={\frac{\partial \widehat{{\mathbf{f}}^{\text{n}}}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}+{\frac{\partial \widehat{{\mathbf{f}}^{\text{t}}}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}$$

$$\begin{gathered} {\frac{\partial \widehat{{\mathbf{f}}^{\text{n}}}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}=-[{\left(\frac{\partial \widehat{f^{\text{n}}}(\mathbf{n})}{\partial t}\right)}_{\text{matl}}-\widehat{f^{\text{n}}}(\mathbf{n})\left({\nabla }_{{}_{\Omega }}\cdot \dot{\mathbf{n}}(\mathbf{n})\right) \\ -\dot{\mathbf{n}}(\mathbf{n})\cdot \left({\nabla }_{{}_{\Omega }}\widehat{f^{\text{n}}}(\mathbf{n})\right)+{\left(\frac{\partial \widehat{f^{\text{n}}}(\mathbf{n})}{\partial t}\right)}_{\mathrm{diff}}]\mathbf{n} \end{gathered}$$

$$\begin{array}{rl}{\frac{\partial \widehat{{\mathbf{f}}^{\text{t}}}(\mathbf{n})}{\partial t}|}_{\mathbf{n}}=& {\left(\frac{\partial \widehat{f^{\text{t}}}(\mathbf{n})}{\partial t}\right)}_{\text{matl}}\mathbf{s}(\mathbf{n})\\ & -\widehat{f^{\text{t}}}(\mathbf{n})\left[{\nabla }_{{}_{\Omega }}\cdot \dot{\mathbf{n}}(\mathbf{n})\right]\mathbf{t}(\mathbf{n})\\ & -\dot{\mathbf{n}}(\mathbf{n})\cdot \left({\nabla }_{{}_{\Omega }}\widehat{f^{\text{t}}}(\mathbf{n})\right)\mathbf{t}(\mathbf{n})+{\left(\frac{d\widehat{{\mathbf{f}}^{\text{t}}}\left(\mathbf{n}\right)}{dt}\right)}_{\text{twirl}}\\ & -{\left(\frac{d\widehat{{\mathbf{f}}^{\text{t}}}\left(\mathbf{n}\right)}{dt}\right)}_{\text{veer}}+{\left(\frac{\partial \widehat{f^{\text{t}}}(\mathbf{n})}{\partial t}\right)}_{\mathrm{diff}}\mathbf{t}(\mathbf{n})\end{array}$$

$${\dot{\mathbf{n}}}^m=\dot{\mathbf{n}}({\mathbf{n}}^m)+\delta {\dot{\mathbf{n}}}^m$$

$${\left(\frac{\partial \hat{g}(\mathbf{n},t)}{\partial t}\right)}_{\text{diff}}=D_g{\nabla }^2\hat{g}(\mathbf{n},t){\dot{ϵ}}_{\text{oct}}$$

$${\left(\frac{\partial \widehat{f^{\mathrm{n}}}(\mathbf{n},t)}{\partial t}\right)}_{\text{diff}}=D_g\overline{f^{\text{n}}}\left(\mathbf{n},t\right){\nabla }^2\hat{g}(\mathbf{n},t){\dot{ϵ}}_{\mathrm{oct}}$$

$${\left(\frac{\partial \widehat{f^{\mathrm{t}}}(\mathbf{n},t)}{\partial t}\right)}_{\text{diff}}=D_g\overline{f^{\text{t}}}\left(\mathbf{n},t\right){\nabla }^2\hat{g}(\mathbf{n},t){\dot{ϵ}}_{\mathrm{oct}}$$

$$\begin{array}{rl}\dot{\mathbf{n}}(\mathbf{n})\stackrel{\stackrel{\mathrm{aligned}}{\mathrm{with}}}{⟷}{\dot{\mathbf{n}}}^{\mathrm{proj}}(\mathbf{n})& \equiv {\mathbf{P}}^{\text{n}}(\mathbf{n})\cdot (\mathbf{D}\cdot \mathbf{n})\\ & =(\mathbf{I}-\mathbf{n}\otimes \mathbf{n})\cdot (\mathbf{D}\cdot \mathbf{n})\end{array}$$

$$\dot{\mathbf{n}}\approx \alpha {\dot{\mathbf{n}}}^{\mathrm{proj}}$$

$${\nabla }_{{}_{\Omega }}\cdot \dot{\mathbf{n}}\approx {\nabla }_{{}_{\Omega }}\cdot (\alpha {\dot{\mathbf{n}}}^{\mathrm{proj}})=-3\alpha \mathbf{n}\cdot ({\mathbf{D}}^{\prime }\cdot \mathbf{n})$$

$${\left(\frac{\partial \overline{f^{\text{n}}}(\mathbf{n})}{\partial t}\right)}_{\text{matl}}=-\overline{k^{\text{n}}}(\mathbf{n})\frac{\partial \ell (\mathbf{n})}{\partial t}$$

$$\frac{\partial \ell (\mathbf{n})}{\partial t}\approx {\beta }^{\text{n}}\overline{\ell }\mathbf{n}\cdot (\mathbf{D}\cdot \mathbf{n})$$

$$k^{\text{n},m}={\left[\frac{3G^2{R}^m{f}^{\text{n},m}}{(1-\nu {)}^2}\right]}^{1/3}$$

$$\overline{k^{\text{n}}}(\mathbf{n})\approx {\left[\frac{3G^2\overline{\ell }\overline{f^{\text{n}}}(\mathbf{n})}{2(1-\nu {)}^2}\right]}^{1/3}$$

Full text

††journal: Mechanics of Materials

1 Introduction

A granular material, such as a soil, possesses an inherent fabric that is imprinted during its formation. This initial fabric usually imparts a directional, anisotropic character, so that strength and stiffness depend upon the direction of loading relative to the original deposition (Arthur and Menzies, 1972). Oda (1972b) defined fabric as the spatial arrangement of particles and voids, and he showed that the microscopic fabric can exhibit two forms of anisotropy: a preferred orientation of elongated or flat particles, and a prevalence of inter-particle contacts in preferred directions. In pioneering experiments, he also found that the subsequent loading and deformation of a soil alters the arrangement of its particles and causes contacts to become increasingly aligned in the direction of the major principal stress — a phenomenon now termed stress-induced fabric anisotropy (Oda, 1972a, c). These early studies, which derived from a geometric view of fabric, were later augmented to include a kinetic or statical aspect of anisotropy. Photo-elastic experiments have shown that forces between particles become largest at contacts that are aligned with the loading direction; moreover, Cundall and others have shown that deviatoric stress is largely an expression of such stress-induced force anisotropy (Cundall and Strack, 1983; Thornton and Barnes, 1986).

The current study addresses the induced anisotropies of contact orientation and contact force. The Paper develops a systematic means of tracking and analyzing the evolution of anisotropy. That is, the focus is on the underlying source and rate of induced anisotropy — on the rate of its evolution rather than its state at a particular instant. Section 2 begins with a brief presentation conventional measures of fabric and stress and then presents a rational, mathematical approach to the manner in which fabric anisotropy emerges and evolves during loading. The basis is an idea recently proposed by Ma and Zhang (2006), and the Paper goes beyond this original idea by supplying additional terms and developing expressions for all terms. Although the first part of the paper is primarily analytical, in Section 3 the Paper proposes specific forms for the various mathematical terms and uses discrete element (DEM) simulations to verify and quantify these terms. The DEM simulations are of densely packed durable spheres, and the data analysis concentrates on soil behavior at large, failure strains, specifically on fabric and strength at the critical state. Once the theory has been developed and quantified, Section 4 applies the theory to predict the effect of the intermediate principal stress on soil strength at the critical state. The theory provides an explanation for the shape of the critical state yield surface, and the predictions are compared with published results and with DEM simulations.

2 Fabric rate equations

The Paper pursues two anisotropies, of contact orientation and of contact force, and models their evolution during bulk loading. Fabric anisotropy will be expressed in terms of the contact density $\widehat{g}(\mathbf{n})$, a function that describes contact orientation within an assembly of particles. The individual orientation $\mathbf{n}^{m}$ of a single $m^{\mathrm{th}}$ contact is the unit vector normal to the surfaces of two particles at their contact point (Fig. 1).

The density $\widehat{g}(\mathbf{n})$ within an entire assembly is the average number of contacts having a common orientation $\mathbf{n}$, as expressed per particle and per unit of area on the unit sphere (i.e., per steradian of solid angle). If a granular assembly contains $N$ particles having $M$ contacts, then the bulk average of contacts per particle, $M/N$, is the integral of density $\widehat{g}(\mathbf{n})$ across the unit sphere,

[TABLE]

where $\Omega$ is the unit sphere surface and $\widehat{g_{\text{avg}}}$ is the average density. Note that the more commonly used average coordination number is simply twice $M/N$. An initially isotropic particle arrangement, with no preferred direction of contact orientation, has a uniform density $\widehat{g}(\mathbf{n})=M/N/(4\pi)$ across all of $\Omega$, measured in contacts per particle per steradian of solid angle. Upon loading, the density becomes anisotropic, with larger values in the direction of compressive loading.

In a similar manner, we can consider a “$\>\widehat{\>}\>$” density of contact force, $\widehat{\mathbf{f}}(\mathbf{n})$, also a function of contact orientation $\mathbf{n}$. This vector density has units of force per particle per steradian on the unit sphere. The average force among those contacts having a particular orientation $\mathbf{n}$ will be written as $\overline{\mathbf{f}}(\mathbf{n})$ and is simply the force density divided by the contact density:

[TABLE]

To relate these densities to stress, we begin with the Cauchy formula for the average stress $\boldsymbol{\sigma}$ within a granular region

[TABLE]

a discrete sum of dyadic products $\mathbf{l}^{m}\otimes\mathbf{f}^{m}$ (= $l_{i}^{m}f_{j}^{m}$) for the $M$ contacts within the region’s volume $V$. Vector $\mathbf{f}^{m}$ is the contact force, and $\mathbf{l}^{m}$ is the branch vector between the centers of two contacting particles (Fig. 1). For spherical particles, vector $\mathbf{l}^{m}$ is the product $\ell^{m}\mathbf{n}^{m}$ of the branch length $\ell^{m}$ and the contact’s unit normal vector $\mathbf{n}^{m}$, so that

[TABLE]

where the approximation includes the average branch length $\overline{\ell}$ among all contacts. In making this approximation, we ignore the small correlation between individual lengths $\ell^{m}$ and the corresponding contact directions and forces. Noting that $\widehat{\mathbf{f}}(\mathbf{n})$ is the density of contact force per particle, the average stress within a sphere assembly can be expressed as

[TABLE]

with $\mathbf{n}\otimes\widehat{\mathbf{f}}(\mathbf{n})=n_{i}\widehat{f}_{j}(\mathbf{n})$. In this form, we replace the discrete sum in Eq. (2) with an integral of density on the unit sphere $\Omega$, and we introduce $\overline{v}$, the average volume of a particle and its associated void space ($\overline{v}=V/N$). Because the force density $\widehat{\mathbf{f}}(\mathbf{n})$ is the product $\widehat{g}(\mathbf{n})\overline{\mathbf{f}}(\mathbf{n})$, deviatoric stress is the result of anisotropies in both contact orientation $\widehat{g}(\mathbf{n})$ and average contact force $\overline{\mathbf{f}}(\mathbf{n})$ (see Rothenburg and Bathurst, 1989).

The above principles apply to the state of anisotropy at any instant and are well established in the literature (for reviews, see Oda and Iwashita 1999; Nemat-Nasser 2004). Equation (3) will henceforth be used to investigate the rates of fabric and stress evolution during loading, the primary intent of the Paper. The stress rate will depend, in part, upon the rate at which the contact force density $\widehat{\mathbf{f}}(\mathbf{n})$ evolves, which will be written as $\left.\partial\widehat{\mathbf{f}}(\mathbf{n},t)/\partial t\right|_{\mathbf{n}}$, or simply $\partial\widehat{\mathbf{f}}(\mathbf{n})/\partial t$. The partial derivatives emphasize that the rate of density $\widehat{\mathbf{f}}(\mathbf{n})$ is measured at fixed orientations $\mathbf{n}$, even though the motions of individual particles will cause contacts to pass through any given orientation. Focusing attention on small fixed portions $d\Omega$ of the unit sphere, the stress rate $\dot{\boldsymbol{\sigma}}$ is computed from the rate at which force density $\widehat{\mathbf{f}}(\mathbf{n})$ changes within these regions:

[TABLE]

again assuming spherical particles. This rate expression includes a possible volume change during the loading process ($\dot{v}=d(V/N)/dt=d\overline{v}/dt$) but assumes that the average branch length $\overline{\ell}$ remains constant, an assumption that is appropriate for hard particles having a ratio of stress and shear modulus $\sigma/G$ that is small.

Particles will roll and slide across each other during loading, and the orientation of any particular contact, say $\mathbf{n}^{m}$, can shift to a new orientation $\mathbf{n}^{m}+\dot{\mathbf{n}}^{m}dt$. These contact movements can be highly erratic, but when tracked by the thousands, contacts are observed to migrate from directions of bulk compression toward directions of extension. This average, prevailing migration rate, denoted as $\dot{\mathbf{n}}(\mathbf{n})$, is a vector field tangent to (and on) the unit sphere and has units of radians per time. The migration $\dot{\mathbf{n}}(\mathbf{n})$ is a function of orientation $\mathbf{n}$ and will depend upon the loading process. As an example, Fig. 2 illustrates the average migration vectors among over 800,000 contacts, as measured in DEM simulations of sustained flow during biaxial plane-strain compression at the critical (steady) state (see Sections 3.1–3.2).

As permitted by the symmetry of these loading conditions, the results have been folded into a single octant of the unit sphere. Contacts are seen to migrate (flow) from the compression direction, $\mathbf{n=\mathbf{e}}_{1}$, toward the zero-strain direction $\mathbf{\mathbf{e}}_{2}$ and toward the extension direction $\mathbf{\mathbf{e}}_{3}$. When considered alone, this migration will be seen to have a softening effect: by transporting contacts and their forces from the direction of compression toward the direction of extension, migration usually diminishes both fabric anisotropy and deviatoric stress. A functional form of the migration $\dot{\mathbf{n}}(\mathbf{n})$ is proposed in Section 3.2.

The stress evolution in Eq. (4) depends on the rate of force density $\partial\widehat{\mathbf{f}}(\mathbf{n})/\partial t$, a rate that will depend upon the interactions of the particles and also upon the prevailing contact migration $\dot{\mathbf{n}}(\mathbf{n})$. The density rate can be viewed as a transport problem on the surface of the unit sphere. In this sense, the rate $\partial\widehat{\mathbf{f}}(\mathbf{n},t)/\partial t$ at a given, fixed orientation $\mathbf{n}$ is described by the Fokker-Planck equation with an additional source density:

[TABLE]

where $\boldsymbol{\nabla}_{\!{}_{\Omega}}$ is the gradient $\partial(\bullet)/\partial\mathbf{n}$ on the unit sphere, and $\boldsymbol{\nabla}_{\!{}_{\Omega}}\cdot(\bullet)$ is the corresponding divergence, $\boldsymbol{\nabla}_{\!{}_{\Omega}}\cdot\mathbf{a}=a_{k,k}$. This general form, except for the final (diffusion) term and certain details among the other terms, was recently proposed by Ma and Zhang (2006). We describe the nature of each term in the remainder of this section, noting that all are amenable to rational analysis and to experimental measurement. In Section 3, we propose more specific forms for each term, verify these forms with DEM data, and quantify the relevant material parameters that appear within the specific forms. As a complement to the force rate in Eq. (5), we must also consider the rate of contact density, which is also the result of three terms:

[TABLE]

where the gradient $\boldsymbol{\nabla}_{\!{}_{\Omega}}\cdot(\dot{\mathbf{n}}\,\widehat{g})=(\dot{n}_{i}\widehat{g})_{,i}$. Didwania et al. (2001) proposed an expression similar to Eq. (6) for the evolution contact density, although the mean-field rotation of particle pairs was used in place of the more general migration $\dot{\mathbf{n}}$.

The first terms on the right of Eqs. (5) and (6), both $d\widehat{\mathbf{f}}/dt$ and $d\widehat{g}/dt$, are source densities, representing the rates of force and contact density that are generated by the particle interactions. The other terms account for convection and diffusion processes which are associated with contact migration, as will be discussed later.

Because contact density $\widehat{g}(\mathbf{n})$ is a scalar field, its “$d\,$” source density $d\widehat{g}(\mathbf{n})/dt$ is somewhat easier to describe and measure than is the vector force rate $d\widehat{\mathbf{f}}/dt$. Contacts are continually created and broken during deviatoric loading, producing the net rate $d\widehat{g}(\mathbf{n})/dt$ at any given orientation $\mathbf{n}$. Unlike the force rate, the contact rate is entirely a material rate, $d\widehat{g}(\mathbf{n})/dt=(\partial\widehat{g}(\mathbf{n})/\partial t)_{\text{matl}}$, and will be expressed as such in the remainder of the paper. As an example, in our DEM simulations of dense particle packings, 4096 spherical particles would typically touch at about 8190 contacts while the assembly was flowing at the critical (steady) state in sustained biaxial plane-strain compression (Section 3.1). During such flow, over 6000 contacts were created and about the same number were broken with each 1% of continued deviatoric strain, yet maintaining the nearly constant 8190 contacts at any instant. Contacts were predominantly created in the compressive strain direction — at orientation $\mathbf{n=\mathbf{e}}_{1}$ in the simulations (Fig. 2) — and were predominantly extinguished in the direction of extension, $\mathbf{n=\mathbf{e}}_{3}$. This pattern of contact activity produces the net material rate and is responsible for induced fabric anisotropy. In Section 3.5, we quantify the material rate, which is shown to be a function of orientation $\mathbf{n}$, the loading conditions, and the material characteristics of the particles.

The source density of force, the $d\widehat{\mathbf{f}}/dt$ rate in Eq. (5) is similar to the contact rate $d\widehat{g}/dt$ but also involves rotation of a vector field — a rotation that will always accompany contact migration. For example, a single contact force $\mathbf{f}^{m}$ between two particles will rotate as the particles roll across each other, inducing a certain force increment during the contact rotation $d\mathbf{n}^{m}$ (see Kuhn and Chang, 2006). This induced rotational increment must be added to any material increment in the force that might result from a changing indentation at the contact:

[TABLE]

using the cross product $\mathbf{a}\times\mathbf{b}=e_{ijk}a_{j}b_{k}$ and inner product $\mathbf{a}\cdot\mathbf{b}=a_{i}b_{i}$. The first term is the material change produced by the indentation process; the second term is an induced increment produced by any tilting increment $d\mathbf{n}^{m}$ of the particle pair; and the final term is an induced twirling increment which accompanies any rigid rotation of the two particles, $d\boldsymbol{\theta}^{i}$ and $d\boldsymbol{\theta}^{j}$, about their common normal axis $\mathbf{n}^{m}$. The source density $d\widehat{\mathbf{f}}(\mathbf{n})/dt$ can likewise be viewed as the sum of a material rate and an induced rotation. This analysis is made easier by first treating force density as the sum of two parts: a normal force density and a tangential force density,

[TABLE]

where $\widehat{f^{\textrm{n}}}(\mathbf{n})$ and $\widehat{f^{\textrm{t}}}(\mathbf{n})$ are scalar densities on the unit sphere, and $\mathbf{t}(\mathbf{n})$ is the average unit direction of tangential force at the particular orientation $\mathbf{n}$. A compressive normal force $\widehat{f^{\textrm{n}}}$ is considered positive. With this approach, the source density in Eq. (5) is the sum of two rates,

[TABLE]

and each vector rate is the sum of a material rate and an induced rotation, as in Eq. (7),

[TABLE]

[TABLE]

The vector field $\mathbf{s}(\mathbf{n})$ is the unit direction of the tangential material rate, a direction that might differ from the current direction $\mathbf{t}(\mathbf{n})$ of tangential force $\widehat{\mathbf{f}^{\text{t}}}(\mathbf{n})$. The final term in Eq. (11) is the net effect of the contacts’ twirling upon the force density. Twirling takes place within the tangent plane and does not alter the normal force density, so this effect is absent in Eq. (10).

We now consider the second terms on the right of Eqs. (5) and (6), which include the gradient $\boldsymbol{\nabla}_{\!{}_{\Omega}}$ applied on the curved surface of the unit sphere. These terms account for divergence and convection effects. For example, the gradient of contact density $\widehat{g}(\mathbf{n})$ in Eq. (6), as introduced by Ma and Zhang (2006), can be expanded as two sub-terms: $\boldsymbol{\nabla}_{\!{}_{\Omega}}\cdot(\dot{\mathbf{n}}\,\widehat{g}(\mathbf{n}))=(\boldsymbol{\nabla}_{\!{}_{\Omega}}\cdot\dot{\mathbf{n}})\,\widehat{g}(\mathbf{n})+\dot{\mathbf{n}}\cdot\left(\boldsymbol{\nabla}_{\!{}_{\Omega}}\widehat{g}(\mathbf{n})\right)$. The first of these sub-terms, the divergence rate $(\boldsymbol{\nabla}_{\!{}_{\Omega}}\cdot\dot{\mathbf{n}})\,\widehat{g}(\mathbf{n})=\dot{n}_{i,i}\widehat{g}(\mathbf{n})$, accounts for either a spreading or converging migration that will rarefy or accumulate contacts at particular orientations $\mathbf{n}$. In Fig. 2, for example, migration spreads (diverges) from the compressive $\mathbf{\mathbf{e}}_{1}$ direction and converges into the extensional $\mathbf{\mathbf{e}}_{3}$ direction. The other sub-term, the convective rate $\dot{\mathbf{n}}\cdot\boldsymbol{\nabla}_{\!{}_{\Omega}}\widehat{g}(\mathbf{n})=\dot{n}_{i}\widehat{g}_{,i}$, addresses the drift of contact density: for example, the convection of a higher contact density, moving at rate $\dot{\mathbf{n}}$, toward a lower density, thereby increasing contact density at the latter orientation.

The gradient of force density that appears in Eq. (5) is similar to the gradient of contact density in Eq. (6), except that force density is a vector field, whose gradient can be analyzed by separating force into its normal and tangential parts, as in Eq. (8),

[TABLE]

The third and sixth terms can be written as

[TABLE]

Because $\dot{\mathbf{n}}$ is normal to $\mathbf{n}$, the quantity $\widehat{f^{\textrm{n}}}(\mathbf{n})\,\dot{\mathbf{n}}$ represents a tangential rotation of normal force along a migration path $\dot{\mathbf{n}}$. The vector field $\dot{\mathbf{t}}(\mathbf{n})$ in Eq. (12) is the change in direction of the tangential force density along contact migration paths. This vector rate is orthogonal to $\mathbf{t}(\mathbf{n})$, so that $\dot{\mathbf{t}}(\mathbf{n})$ will have a normal component and, possibly, a tangential component. A tangential component of $\dot{\mathbf{t}}(\mathbf{n})$ occurs at any orientation where $\mathbf{t}$ veers in direction along a path $\dot{\mathbf{n}}$: in Fig. 2, for example, migration arrows are seen to veer toward the east as they approach the $\mathbf{e}_{2}$–$\mathbf{e}_{3}$ equator. The two parts of $\widehat{f^{\textrm{t}}}(\mathbf{n})\dot{\mathbf{t}}(\mathbf{n})$, normal and veering, can be written as

[TABLE]

The final, veering part lies in the tangent plane but is orthogonal to the tangential direction $\mathbf{t}(\mathbf{n})$, as will be demonstrated in an analysis of DEM data in Section 3.4.1.

The full rate of force density that produces the stress rate in Eq. (4) can be separated into normal and tangential parts,

[TABLE]

The two parts are expanded by combining Eqs. (5) and (9)–(15):

[TABLE]

and

[TABLE]

in which we have applied the identity $\mathbf{a}\times(\mathbf{b}\times\mathbf{c})=\mathbf{b}(\mathbf{a}\cdot\mathbf{c})-\mathbf{c}(\mathbf{a}\cdot\mathbf{b})$, noting that $\mathbf{t}\cdot\mathbf{n}=0$, $\dot{\mathbf{n}}\cdot\mathbf{n}=0$, and $\mathbf{n}\cdot\mathbf{n}=1$. The forms of the material force rates $\left(\bullet\right)_{\mathrm{matl}}$ and of the twirling and veering rates are settled by analyzing DEM results in Sections 3.3 and 3.4. Equations (17) and (18) differ from those of Ma and Zhang (2006), with the inclusion of the twirling, veering, and diffusions terms and the mutual cancelling of the tilting terms that arise in Eqs. (9)–(15).

Until now, we have considered the average contact migration $\dot{\mathbf{n}}$ and its effect on fabric and stress rates. We must also consider the random fluctuations among individual contact motions — fluctuations that produce the diffusion terms in Eqs. (6), (17), and (18). The tangential rate $\dot{\mathbf{n}}^{m}$ of an individual contact $m$, oriented in the direction $\mathbf{n}^{m}$, is the sum of the prevailing (mean) migration $\dot{\mathbf{n}}$ and the individual’s fluctuation $\delta\dot{\mathbf{n}}^{m}$ from the mean:

[TABLE]

These fluctuations are quite large and can produce a diffusion of the contact density $\widehat{g}$ during bulk deformation. Contact diffusion at an orientation $\mathbf{n}$ — the final term in Eq. (6) — is driven by the ongoing deformation of the bulk material, causing contacts to diffuse (disperse) from orientations of high contact concentration toward orientations of lower concentration. This phenomenon is distinct from convection and divergence — the “$\boldsymbol{\nabla}_{\!{}_{\Omega}}$” term in Eq. (6) — in which contacts are swept along by the prevailing rate $\dot{\mathbf{n}}$. Contact diffusion can be modelled with the classical diffusion equation, in the form

[TABLE]

where $\nabla^{2}\widehat{g}=\widehat{g}_{,kk}$ is the Laplacian of $\widehat{g}$ on the surface of the unit sphere and $D_{g}$ is the diffusion coefficient which quantifies the process. Because bulk deformation drives the diffusion, we use the instantaneous octahedral strain rate $\dot{\epsilon}_{\text{oct}}$ as a scalar measure of the distortional strain ($\dot{\epsilon}_{\text{oct}}=\sqrt{D_{ij}^{\prime}D_{ij}^{\prime}/3}$; $D_{ij}^{\prime}=D_{ij}-\frac{1}{3}D_{kk}\delta_{ij}$), although other measures could be used as well. In Section 3.6, we provide further reasoning for the form (20) and present the means of extracting the diffusion coefficient $D_{g}$ from DEM simulations.

A diffusion of force density $\widehat{\mathbf{f}}(\mathbf{n})$ will accompany the diffusion of contact density $\widehat{g}(\mathbf{n})$: contact forces are dispersed with their contacts during deformation. The following forms of force diffusion accrue from Eq. (1) and complement the contact diffusion of Eq. (20):

[TABLE]

[TABLE]

which define the scalar diffusion rates that appear in Eqs. (17) and (18).

To summarize this section, the rates of contact and force densities include material rates, combined with divergence, twirling, and diffusion effects (Eqs. 5 and 6). As will be seen, a sort of competition exists between the material rate and the other effects. Anisotropy is usually reduced by divergence, twirling, and diffusion, as these effects tend to diminish $\widehat{g}$ and $\widehat{\mathbf{f}}$ at those orientations $\mathbf{n}$ where the densities are large. At the same orientations, the densities are usually replenished by the material rates, which represent a generation of contacts and of contact force. The material rate will be seen to dominate at the start of loading, inducing fabric anisotropy and deviatoric stress. During failure, the material rate becomes quite small and is counteracted by the other rate effects, eventually leading to a steady, critical state of fabric and stress.

3 Quantifying fabric evolution

In the previous section, we found that the stress rate $\dot{\boldsymbol{\sigma}}$ results from changes in the contact force density $\widehat{\mathbf{f}}(\mathbf{n})$ — changes that can be separated into several rate fields. In this section, specific forms of the various fields are adopted, guided by DEM observations of granular behavior. Taken together, these forms comprise the rudiments of a constitutive model for soils and other granular materials, a model based on micro-mechanics and informed by DEM results. The Paper focuses on behavior at the critical state, when granular materials reach a stationary condition: during sustained steady state flow, the volume, the stress, the fabric, and, most notably, the density functions $\widehat{g}(\mathbf{n})$, $\widehat{f^{\text{n}}}(\mathbf{n})$, and $\widehat{\mathbf{f}^{\text{t}}}(\mathbf{n})$ remain constant. To aid understanding of the critical state behavior, we also consider the rate of $\widehat{f^{\text{n}}}(\mathbf{n})$ at the other extreme of deformation — at the start of loading — which will provide a basis for quantifying the corresponding critical state terms. We begin with a brief description of the simulations.

3.1 DEM simulations

DEM simulations were performed on twenty small assemblies of spherical particles. The simulations permitted the observation of a sufficiently large number of contacts to attain their motions, force rates, and net creation rates across the entire unit sphere of orientations. The assemblies contained the same set of 4096 particles that were randomly packed into cube containers having periodic boundaries on all sides and the dense initial conditions listed in Table 1.

The spheres were polydisperse with diameters ranging from 0.4$D{}_{50}$ to 1.5$D{}_{50}$, where $D{}_{50}$ is the median diameter. The hard non-breaking particles interact at Hertz contacts having a frictional limit $\mu=0.5$ and a modified Mindlin tangential stiffness, as described by Lin and Ng (1997). The cubic assemblies had dimensions of about $13.4\times 13.4\times 13.4$ particle diameters: small enough to prevent shear bands, yet large enough to capture the average, bulk material behavior. Because a small assembly of only 4096 particles will exhibit substantial spikes in stress during deviatoric loading, twenty different assemblies were randomly created and then loaded. Their averaged behavior is reported herein. Although several loading paths are studied, the primary loading was slow biaxial plane-strain compression: the assemblies were compressed by continually reducing their $x_{1}$ dimension at a constant rate ($\dot{\epsilon}_{11}=\textrm{constant}$) while maintaining a constant normal stress in the $x_{3}$ direction and a constant assembly width in the $x_{2}$ direction ($\sigma_{33}$ constant $=p_{\text{o}}=320\>\text{kPa}$, and $\epsilon_{22}=0$; see inset in Fig. 2). Because expansion was permitted against a constant stress $\sigma_{33}$, the simulations can be considered as drained tests in the geotechnical sense. The averaged results are displayed in Fig. 3, which shows the normalized deviator stress $(\sigma_{11}-\sigma_{33})/p_{\text{o}}$ and the fabric anisotropy $F_{11}-F_{33}$ during the course of plane-strain biaxial compression (the fabric measure $F_{ij}$ is defined as $\int_{\Omega}n_{i}n_{j}\widehat{g}(\mathbf{\mathbf{n})\,}d\Omega$, as in Nemat-Nasser 2004).

Because we are interested in the rates of fabric and stress evolution at both the macro (bulk) and micro scales, contact statistics were compiled during loading at the critical state (at $-\epsilon_{11}>0.25)$. “Snapshots” of the contacts and their rates were taken at several such strains, and by doing the same for all twenty assemblies, we were able to analyze over 800,000 contacts, as discussed below.

3.2 Migration rate $\dot{\mathbf{n}}$,

convection, and divergence

During deviatoric loading, particles roll and slide across each other, causing contact orientations to migrate. DEM simulations were used to determine a functional form of the average migration $\dot{\mathbf{n}}(\mathbf{n})$ in relation to the bulk deformation rate. As an example, the migrations of 800,000 contacts were measured at the critical state in simulations of plane-strain biaxial compression, and the average rates are depicted in Fig. 2. Although somewhat obscured, two arrows emanate from each grid point. The heavier arrows are the actual, measured rates. These arrows are closely aligned with lighter arrows that represent a certain projection $\dot{\mathbf{n}}^{\mathrm{proj}}$ of the instantaneous deformation rate $\mathbf{D}$ onto the unit sphere:

[TABLE]

with $\mathbf{P}^{\text{n}}\cdot\mathbf{D}\cdot\mathbf{n}=P_{ij}^{\text{n}}D_{jk}n_{k}$. In this equation, $\mathbf{D}$ is the rate of deformation tensor (instantaneous strain rate), and the matrix $\mathbf{P}^{\text{n}}(\mathbf{n})=\delta_{ij}-n_{i}n_{j}$ projects the rate vector $\mathbf{D}\cdot\mathbf{n}$ ($=D_{jk}n_{k}$) onto the tangent plane. The rate $\dot{\mathbf{n}}^{\mathrm{proj}}$ represents the ideal tangential motion of two contacting spheres whose centers move in perfect accord with the mean, bulk deformation $\mathbf{D}$. Although individual contacts migrate in widely varying directions and rates, the alignment of the average rate $\dot{\mathbf{n}}$ and the ideal rate $\dot{\mathbf{n}}^{\mathrm{proj}}$ is quite close: the two directions differ, on average, by less than $3^{\circ}$.

Although the observed and ideal rates are aligned, they differ in magnitude. Different scales have been applied in Fig. 2 for displaying the measured rates $\dot{\mathbf{n}}$ and the projected rates $\dot{\mathbf{n}}^{\mathrm{proj}}$, and although obscured in the small figure, the lengths of $\dot{\mathbf{n}}$ and $\dot{\mathbf{n}}^{\mathrm{proj}}$ are consistently in about the same proportion: the two vector fields are nearly aligned and proportional. This observation suggests that the migration rate $\dot{\mathbf{n}}$ can be approximated as

[TABLE]

a condition that is closely held throughout the loading process. At the start of loading, the simulations show that the factor $\alpha$ is equal to 1.0, so that the measured and ideal rates are about equal (the assumption of $\dot{\mathbf{n}}=\dot{\mathbf{n}}^{\mathrm{proj}}$ being made by Didwania et al. 2001). As deformation progresses to larger strains, the actual rate exceeds the ideal projected rate, with $\alpha=1.7$ at the critical state.

The approximation of Eq. (24) can be used to estimate the convection and divergence rates in Eqs. (6), (17), and (18). For example, the scalar divergence of $\dot{\mathbf{n}}$ is

[TABLE]

or $-3\alpha D_{ij}^{\prime}n_{i}n_{j}$, where $\mathbf{D}^{\prime}$ ($=D_{ij}-D_{kk}\delta_{ij}/3$) is the deviatoric part of deformation.

3.3 Rate of normal force

Particles press against each other with changing force as a granular material undergoes bulk deformation. The material rate of normal force density, $(\partial\widehat{f^{\textrm{n}}}(\mathbf{n})/\partial t)_{\textrm{matl}}$, is the net effect of these many changes at particular orientations $\mathbf{n}$. To quantify this material rate, several hundreds of thousands of contacts were observed in DEM simulations of biaxial plane-strain loading at two extremes of deformation: at the start of loading and during sustained flow at the critical state. In general, the material rate exhibits the following characteristic: the rate is usually positive (increasingly compressive) at orientations $\mathbf{n}$ in which the bulk strain produces compression; whereas, the material rate is negative (tensile) in directions of extension. This observation, although not surprising, suggests that the average “$\,\overline{\rule{0.0pt}{4.30554pt}\;\;}\,$” compressive rate, $(\partial\overline{f^{\textrm{n}}}(\mathbf{n})/\partial t)_{\textrm{matl}}$, might be approximated as the product of an average normal stiffness $\overline{k^{\textrm{n}}}$ and the average rate of approach between the centers of contacting particles:

[TABLE]

If the particle motions were to conform to the mean, bulk rate of deformation field $\mathbf{D}$, the rate of approach would be $\partial\ell(\mathbf{n})/\partial t=\overline{\ell}(\mathbf{n})\,\mathbf{n}\cdot(\mathbf{D}\cdot\mathbf{n})=\overline{\ell}(\mathbf{n})\,n_{i}D_{ij}n_{j}$, where $\overline{\ell}(\mathbf{n})$ is the average distance between the centers of particles oriented in direction $\mathbf{n}$. Numerous studies have sought expressions for the bulk elastic moduli of granular media by starting with this mean-field assumption. This approach over-estimates the moduli and is usually amended by considering motion fluctuations from the mean (e.g. Jenkins et al. 2005). Successful estimates are only achieved, however, at small strains and while behavior is elastic. Observations have shown that the normal motions between particles, $\partial\ell(\mathbf{n})/\partial t$, are suppressed during deformation — especially at large strain — as particles tend to roll and slide in a manner that minimizes such motion (Kuhn and Bagi, 2004). In this regard, we will introduce a factor $\beta^{\text{n}}$ to reduce the indentation rate between particles. With this change, the scalar rate $\partial\ell(\mathbf{n})/\partial t$ is approximated as

[TABLE]

In this approximation, we also use the average branch length $\overline{\ell}$ in place of the function $\overline{\ell}(\mathbf{n})$, ignoring the small correlation between branch length and orientation.

The average contact stiffness $\overline{k^{\text{n}}}$ in Eq. (26) will depend upon the stiffness characteristics of the particles themselves. For a single $m$th contact between two isotropic elastic spheres of equal size, the Hertz stiffness is

[TABLE]

where $f^{\textrm{n},m}$ is the pair’s current normal (compressive) force, $G$ is the particle shear modulus, $\nu$ is the Poisson’s ratio, and $R^{m}$ is the shared radius. This contact stiffness depends upon the contact force $f^{\textrm{n},m}$. Because the average normal force $\overline{f^{\textrm{n}}}(\mathbf{n})$ within an entire assembly is known to be anisotropic, we would expect the average contact stiffness $\overline{k^{\text{n}}}$ to depend on orientation. The average stiffness among all contacts that share an orientation $\mathbf{n}$ can be approximated as

[TABLE]

using $\overline{\ell}/2$ in place of $R$. By combining Eqs. (26), (27), and (29) and twice applying Eq. (1), we can estimate the material rate of normal force density as

[TABLE]

with the stiffness density

[TABLE]

Coefficient $\beta^{\text{n}}$ was investigated for two extremes of loading: during initial loading and during sustained plastic flow at the critical state.

3.3.1 Rate of normal force at small strains

Upon initial loading of dense DEM sphere packings, the value $\beta^{\text{n}}$ was measured as 0.94: the actual stiffness was slightly smaller than the ideal stiffness that would apply if each particle moved in perfect accord with the mean deformation field. Figure 4 compares Eq. (30) and the DEM data, with both plotted in a dimensionless form.

Four meridians of the unit sphere are shown in the figure, corresponding to the angle $\xi$ in Fig. 2. Equation (30) closely fits the DEM data.

With a $\beta^{\text{n}}=0.94$ at the start of loading, the material rate of normal force density is more than two hundred times larger than the other rates that contribute to the force density: the divergence, convection, and diffusion terms in Eq. (17). The rapid evolutions of fabric and stress at the start of loading are, therefore, dominated by the material rate $(\partial\widehat{f^{\text{n}}}(\mathbf{n})/\partial t)_{\textrm{matl}}$, with the other rates nearly inconsequential. The situation changes, however, upon further loading. Because $\beta^{\text{n}}$ decreases with increasing strain, its hardening effect is progressively diminished, and the influences of divergence, convection, and diffusion become increasingly more significant — nearly dominant — as will be seen in the next paragraphs.

3.3.2 Rate of normal force at the critical state

DEM simulations were also used to measure the material parameter $\beta^{\text{n}}$ during failure at the critical state. Because the total rate of normal force $\left.\partial\widehat{\mathbf{f}^{\text{n}}}/\partial t\right|_{\mathbf{n}}$ is zero at the critical state, the corresponding material rate can be readily computed from the remaining terms in Eq. (17). A $\beta^{\text{n}}=0.0037$ fits the DEM data, although this value must be slightly amended, as described below. Such a small value of $\beta^{\text{n}}$ indicates that the average normal motions between particles, the rate $\partial\ell(\mathbf{n})/\partial t$ in Eq. (27), is much smaller than would be anticipated by assuming that the particle motions conform to a uniform, affine deformation field. On the other hand, we had previously noted that the average tangential motions become somewhat greater than those of uniform deformation, with an $\alpha$ of 1.7 in Eq. (24). These contrasting results are consistent with other evidence that the large-strain motions of particles are dominated by the tangential rolling of particle pairs, but with minimal average normal motions at the contacts (Kuhn and Bagi, 2004).

The DEM simulations also reveal an unexpected aspect of the material rate $(\partial\widehat{f^{\textrm{n}}}(\mathbf{n})/\partial t)_{\textrm{matl}}$: although one might expect a force rate of zero at those neutral orientations $\mathbf{n}$ where particles neither approach nor withdraw (where $\partial\ell(\mathbf{n})/\partial t=0$), we find instead that the material rate is usually slightly negative (depletive or tensile) at these orientations. This anomalous situation is most noticeable at neutral orientations that also have large migration rates $\dot{\mathbf{n}}(\mathbf{n})$. To account for this observation, we apply a small adjustment to the material rate of Eq. (30), replacing that equation as follows:

[TABLE]

in which the new, subtracted $\beta_{2}^{\text{n}}$ term produces a small depletive bias in the material rate. The magnitude $\left|\dot{\mathbf{n}}(\mathbf{n})\right|$ is approximated with Eq. (24) and is normalized by dividing by the average rate $\dot{\epsilon}_{\mathrm{oct}}=\sqrt{D_{ij}^{\prime}D_{ij}^{\prime}/3}$. The factor $\beta_{2}^{\text{n}}$ was measured as 0.65.

We must also consider the role of the mean stress in generating contact force density. For spherical particles, the mean stress depends exclusively on the normal components of the contact forces (Cundall and Strack, 1983), such that the bulk pressure $p$ is proportional to the average normal force density $\widehat{f^{\text{n}}}(\mathbf{n})$. At small strains, the total rate of normal force density, the rate $\partial\widehat{f^{\text{n}}}(\mathbf{n})/\partial t$ on the left of Eq. (17), is dominated by the material rate, and its approximation with Eqs. (31) and (32) would suggest that the bulk stiffness is proportional to $G^{2/3}$ and to $p^{1/3}$. This scaling at small strains is in fair agreement with small-strain vibrational experiments which show that the elastic moduli are proportional to $p^{n}$ with an exponent $n$ between 1/3 and 1/2 (see Goddard 1990 for a review). Granular behavior at large strains scales quite differently. At the critical state, strength is proportional to the confining pressure, $p^{1}$, and is nearly independent of the particle stiffness $G$. If left unmodified, the stiffness density in Eq. (31) would produce a strength proportional to $p^{1/3}$ and to $G^{2/3}$, contrary to the observed behavior of soils and other granular materials. We should expect, therefore, that at large strains, the factor $\beta^{\text{n}}$ will depend upon $p$ and upon the particles’ elastic properties in the following manner:

[TABLE]

where the dimensionless factor $B^{\text{n}}$ effects a proper scaling of granular strength at large strains. We tested this hypothesis by running DEM simulations of biaxial compression that were identical to those previously described, except that the initial confining pressure was increased about six-fold. At the critical state, the mean stress increased from 490 kPa to 3200 MPa and the strength was found to increase by the same factor, but $\beta^{\text{n}}$ had only increased from 0.0037 to 0.012 — not a six-fold increase, but roughly in accord with Eq. (33) and a $B^{\text{n}}=6.3$. The second parameter in Eq. (32), $\beta_{2}^{\text{n}}$, remained about the same for the two confining pressures, as would be expected, since the last term in Eq. (32) is proportional to $\widehat{f^{\text{n}}}$, a form that is consistent with strength being proportional to mean stress.

In Fig. 5, the combination of Eqs. (32) and (33) is compared with data from 800,000 contacts in DEM simulations.

The simulations are of plane-strain biaxial compression at the critical state, and data is presented along four meridians (angle $\xi$, Fig. 2). All results are reported in a dimensionless form: force density is divided by the stress $p$, and its time rate is normalized with respect to the average octahedral rate $\dot{\epsilon}_{\mathrm{oct}}$. Equations (32) and (33) are in close agreement with the DEM data.

3.4 Rates of tangential force

DEM simulations can be used to resolve a realistic form for the material rate of tangential force at the critical state — the rate $(\partial\widehat{f^{\text{t}}}(\mathbf{n})/\partial t)_{\textrm{matl}}$ in Eq. (18). We will quantify this material rate, after first considering the tangential force rates of twirling and veering.

3.4.1 Twirling and veering rates at the critical

state

During sustained flow at the critical state, the DEM simulations show that the tangential force density $\widehat{f^{\text{t}}}(\mathbf{n})\mathbf{t}(\mathbf{n})$ becomes closely aligned with the direction of contact migration $\dot{\mathbf{n}}(\mathbf{n})$, such that the unit direction is approximated as

[TABLE]

where $\dot{\mathbf{n}}(\mathbf{n})$ is the field depicted in Fig. 2 and given by Eqs. (23) and (24). The total rate of tangential force that appears on the left of Eq. (18) results from various effects: a material rate combined with convection, twirling, veering, and diffusion effects. The density rates of twirling and veering involve rotations of tangential force within the tangent plane. We can isolate and directly measure the twirling rate by considering other DEM simulations that eliminate the veering effect: by using triaxial rather than plane-strain loading. With triaxial compression, assemblies are compressed in the $\mathbf{e}_{1}$ direction while constant stress is maintained in the $\mathbf{e}_{2}$ and $\mathbf{e}_{3}$ directions (see inset, Fig. 6). This symmetric loading condition, in which two principal strains are equal, produces contact migrations $\dot{\mathbf{n}}$, as defined in Eq. (23), along meridians (geodesics) that emanate from the $\mathbf{e}_{1}$ pole and approach the $\mathbf{e}_{2}$–$\mathbf{e}_{3}$ equator. As stated above, $\mathbf{t}(\mathbf{n})$ and $\dot{\mathbf{n}}(\mathbf{n})$ are aligned at the critical state, and because neither direction veers under triaxial loading, we can use triaxial DEM data to directly measure any bulk twirling of the tangential forces.

The twirling rates (i.e., from the final term in Eq. 7) of over 800,000 contacts were averaged, and Fig. 6 illustrates these average density rates.

Although somewhat obscured in the small monochrome figure, two vector fields are displayed: thinner arrows correspond to the twirling density $(\partial\widehat{\mathbf{f}^{\text{t}}}(\mathbf{n})/\partial t)_{\textrm{twirl}}$; thicker arrows are the contact migration field $\dot{\mathbf{n}}(\mathbf{n})$. Some scatter is apparent in the data, but the results indicate that the twirling rate vectors are consistently in a direction opposite the migration field. For example, in the northern octant of Fig. 6, tangential forces tend to rotate toward the north, even as their contacts are migrating toward the south. The magnitude of the twirling density was also found to be roughly proportional to the product of the tangential force density $\widehat{f^{\text{t}}}(\mathbf{n})$ and the magnitude of $\dot{\mathbf{n}}(\mathbf{n})$. These observations suggest the following form of the twirling density field:

[TABLE]

The factor $\gamma_{\text{twirl}}$ was about 1.0 in the DEM simulations. Figure 7 compares this equation with the twirling data for triaxial compression at the critical state.

Equation (35) and a $\gamma_{\text{twirl}}=1.0$ fit the data, although the figure also indicates that the twirling rate of force is relatively small when compared with the normal force rates that are shown in Fig. 5.

Having resolved the twirling effect, we can now investigate possible veering of the tangential force direction $\mathbf{t}(\mathbf{n})$. To this end, we return to biaxial plane-strain compression simulations. The veering rate of tangential force is the tangential component of the rate $\widehat{f^{\text{t}}}\dot{\mathbf{t}}$ along migration paths $\dot{\mathbf{n}}$ (from Eqs. 14 and 15):

[TABLE]

where the projection matrix $\mathbf{P}^{\text{n}}(\mathbf{n})$ extracts the tangential component of this rate (see Eq. 23). As an intermediate step in deriving the veering rate, the unit vector $\mathbf{t}(\mathbf{n})$ can be expressed as the product of a matrix $\mathbf{Q}(\mathbf{n})$ and the unit normal $\mathbf{n}$, as $\mathbf{t}(\mathbf{n})=\mathbf{Q}(\mathbf{n})\cdot\mathbf{n}/|\mathbf{Q}(\mathbf{n})\cdot\mathbf{n}|$, in which $\mathbf{Q}\cdot\mathbf{n}=Q_{ij}n_{j}$. After differentiating $\mathbf{t}(\mathbf{n})$ with respect to $\mathbf{n}$,

[TABLE]

where the new operator $\mathbf{P}^{\text{t}}(\mathbf{n})=\mathbf{I}-\mathbf{t}\otimes\mathbf{t}=\delta_{ij}-t_{i}t_{j}$ projects vectors onto a plane that is perpendicular to the tangent direction $\mathbf{t}$. Combined with $\mathbf{P}^{\text{n}}(\mathbf{n})$, the matrix product $\mathbf{P}^{\text{n}}\cdot\mathbf{P}^{\text{t}}=P_{ik}^{\text{n}}P_{kj}^{\text{t}}$ produces a veering rate that is orthogonal to both $\mathbf{n}$ and $\mathbf{t}$.

As expressed in Eq. (34), the simulations show that during sustained flow at the critical state, the tangential force density is closely aligned with the migration direction $\dot{\mathbf{n}}=\mathbf{P}^{\text{n}}\cdot\mathbf{D}\cdot\mathbf{n}$, defined in Eqs. (23) and (24): that is, matrix $\mathbf{Q}(\mathbf{n})$ is equal to the matrix product $\mathbf{P}^{\text{n}}(\mathbf{n})\cdot\mathbf{D}$. The corresponding veering rate in Eq. (36) is, therefore,

[TABLE]

The plane-strain DEM data reveals behavior that is similar to the triaxial simulations that were discussed earlier. The twirling rate $(\partial\widehat{\mathbf{f}^{\text{t}}}(\mathbf{n})/\partial t)_{\textrm{twirl}}$ for plane-strain biaxial compression is also directed roughly opposite $\dot{\mathbf{n}}$, but with one difference. Although the twirling rate and $\dot{\mathbf{n}}(\mathbf{n})$ are roughly opposed in the biaxial simulations, they are not perfectly counter-aligned: the DEM data show that the twirling rate consistently includes a small tangential component that is orthogonal to the $\mathbf{t}(\mathbf{n})$ direction — a small component that is aligned with the veering rate given in Eq. (37). With this observation, we speculate that the original approximation in Eq. (35) can be applied to non-triaxial loading by making the following adjustment:

[TABLE]

Allowing for considerable scatter in the DEM data, this approximation, with $\gamma_{\text{twirl}}=1.0$, gives a reasonable fit to the DEM data of both triaxial and plane-strain simulations.

3.4.2 Rate of tangential force at the critical state

In a previous section, DEM simulations were used to extract an approximation of the material rate of normal force during sustained deformation at the critical state (Section 3.3.2). The process was aided by the use of Eq. (17) and by the presumption of a vanishing rate $\left.\partial\widehat{\mathbf{f}^{\text{n}}}/\partial t\right|_{\text{n}}$ during steady, critical state loading. We now apply a similar approach to tangential force, by using Eq. (18) to find the tangential vector rate $(\partial\widehat{\mathbf{f}^{\text{t}}}(\mathbf{n})/\partial t)_{\text{matl}}$. Because the other terms in Eq. (18) are aligned with the migration field $\dot{\mathbf{n}}(\mathbf{n})$, the unit direction $\mathbf{s}(\mathbf{n})$ of the tangential material rate must also be aligned with $\dot{\mathbf{n}}(\mathbf{n})$. The DEM experiments show that at certain orientations, the directions $\mathbf{s}(\mathbf{n})$ and $\dot{\mathbf{n}}(\mathbf{n})$ do indeed coincide; however, at other orientations, the simulations show that they are collinear but in opposite directions. This paradoxical observation is reconciled by considering another observation: at those orientations $\mathbf{n}$ where the particles approach each other — when the normal rate $(\partial\widehat{f^{\text{n}}}(\mathbf{n})/\partial t)_{\textrm{matl}}$ is compressive — the direction of the tangential rate coincides with that of $\dot{\mathbf{n}}(\mathbf{n})$. Contrarily, for orientations $\mathbf{n}$ at which particles tend to withdraw from each other, the tangential rate and $\dot{\mathbf{n}}(\mathbf{n})$ are counter-aligned. These observations resemble the behavior of a frictional block system that is loaded with normal and tangential forces, such as that described by Baz̆ant and Cedolin (1991), §10.7. The DEM simulations also show that the magnitude of the tangential material force rate correlates with the magnitude of the migration rate, $|\dot{\mathbf{n}}(\mathbf{n})|$: the material rate is larger at orientations where the particles are migrating more rapidly across each other. When considered together, these observations suggest the following form for the tangential material rate at large strains:

[TABLE]

Although Eq. (39) is consistent with a frictional system, the factor $\zeta^{\text{t}}$ was measured as only 0.067 in the simulations, a value much smaller than the 0.50 friction coefficient between particles. A larger factor might have been realized if frictional sliding were to occur at all contacts, but at the critical state, the DEM simulations show that only 18% of contacts are sliding (see also Thornton 2000), as particles tend to roll rather than slide at their contacts (Kuhn and Bagi 2004).

Equation (39) is compared with data from DEM simulations in Fig. 8, which shows results along two meridians of the unit sphere. Equation (39) is in general agreement with the experimental data.

3.5 Rate of contact creation at the

critical state

The material rate of contact density, ${(\partial\widehat{g}(\mathbf{n})/\partial t)}_{\mathrm{matl}}$, is the net rate at which contacts are created or extinguished during deformation, as in Eq. (6). The DEM simulations show that contacts are predominantly created at orientations $\mathbf{n}$ in which the deformation $\mathbf{D}$ produces compression between particle pairs, whereas contacts are predominantly broken at orientations of extension. These trends resemble the situation with normal force density, so we begin by assuming that the net contact creation rate is proportional to the rate at which particles approach (or withdraw from) each other — the rate $\partial\ell(\mathbf{n})/\partial t$ in Eq. (27). In variance with the force rate, however, we now use the deviatoric rate $\mathbf{D}^{\prime}$ ($=D_{ij}-D_{kk}\delta_{ij}/3$) instead of the full rate $\mathbf{D}$, nullifying the influence of bulk volume change on the rate of contact creation. This modification is justified by two observations. The contact network among particles is most effectively rearranged by bulk distortion: although pure isotropic compression will increase the contact forces, it does not greatly alter the number or orientations of the contacts. Furthermore, the dilation that usually accompanies the shearing of dense granular materials does not appreciably disengage contacts: the number of contacts will typically remain nearly constant during post-peak deformation, even as the material is vigorously dilating (Thornton, 2000). We should expect, however, that the contact density rate $(\partial\widehat{g}(\mathbf{n})/\partial t){}_{\mathrm{matl}}$ will depend on the average normal force $\overline{f^{\textrm{n}}}(\mathbf{n})$: larger contact forces imply greater contact indentations and require a greater motion $d\ell$ to disengage the contacts. In this regard, we introduce a reference movement for the Hertz contact between an $m$th pair of elastic spheres (see Eq. 28),

[TABLE]

This reference movement would cause two spheres, initially pressed together with force $f^{\textrm{n},m}$, to withdraw and disengage, if the original stiffness $k^{\textrm{n},m}$ was active throughout the withdrawal process. The material rate of contact density is approximated as the particle withdrawal rate $\partial\ell(\mathbf{n})/\partial t$ divided by the average reference movement $\overline{f^{\textrm{n}}}(\mathbf{n})/\overline{k^{\text{n}}}$ and multiplied by the current density $\widehat{g}(\mathbf{n})$. Applying Eqs. (1), (27), and (31),

[TABLE]

where $\mathbf{D}^{\prime}$ is used in place of $\mathbf{D}$. This material rate is negative at orientations $\mathbf{n}$ where particles withdraw from each other — when $\mathbf{n}\cdot(\mathbf{D}^{\prime}\cdot\mathbf{n})$ is positive (tensile).

Figure 9 compares Eq. (41) with DEM data of biaxial plane-strain compression at the critical state, showing that the equation closely fits the data.

If we compare the material rates of contact force creation and of contact creation (Eqs. 32 and 41), we see that the material rate of force is proportional to $\widehat{g}^{2/3}\widehat{f^{\text{n}}}^{1/3}$; whereas the material rate of contact creation is proportional to $\widehat{g}^{5/3}\widehat{f^{\text{n}}}^{-2/3}$. Despite the different scaling of these two rates, the DEM simulations reveal that the same $\beta^{\text{n}}$ value applies to both rates: the generation of contacts and the generation of contact force apparently share a common origin.

3.6 Diffusion rate at the critical state

Classical diffusion theory explains the diffusion of a molecular species, either within itself (self-diffusion) or through other molecular species (e.g., Jeans, 1962). The process is driven by random fluctuations among the molecules’ velocities, displacing them from their original positions at time $t=0$. With these random motions, individual displacements $\delta\mathbf{r}^{m}(t)$ increase with time, such that the collective mean-square displacement is roughly proportional to time:

[TABLE]

In this form of Einstein-Smoluchowski diffusion, $d$ is the spatial dimension; the $\delta r_{j}$ are the separate $j$-components of the randomly advancing displacements; and the diffusion coefficient $D$ is a measure of the time rate of these growing displacements. In our application, the displacements $\delta\mathbf{r}^{m}$ are not of molecules or particles; rather, they are the tangential, angular displacements of individual contact orientations $\delta\dot{\mathbf{n}}^{m}$ on the unit sphere — the fluctuations in Eq. (19) — that occur as particles slide or roll across each other during bulk deformation.

Equation (42) provides the means for experimentally measuring a diffusion coefficient — in particular, the coefficient $D_{g}$ of contact self-diffusion in Eqs. (20)–(22).

The tangential movements of several contacts are plotted in Fig. 10a over the course of 0.045 strain at the critical state. The figure illustrates the erratic, zigzag nature of contact movement (i.e., a small mean free path). To measure the coefficient of contact diffusion $D_{g}$, we tracked the long-term movements, $\int\delta\dot{\mathbf{n}}^{m}\,dt$, of over 150,000 contacts as the assembly was being deformed at the critical state. The mean-square cumulative contact displacements are plotted in Fig. 10b as a function of the advancing octahedral strain $\epsilon_{\mathrm{oct}}$. The result is a nearly linear relation, consistent with the conceptual Eq. (42), with cumulative strain replacing time. The slope of 0.03 is the contact diffusion coefficient $D_{g}$ that appears in Eqs. (20)–(22), having units of radians2 per unit of strain $\epsilon_{\mathrm{oct}}$. This value of $D_{g}$ is fairly small when compared with the material rate of $\widehat{g}$ shown in Fig. 9. Didwania et al. (2001) also measured a relatively small rate for a kinematic diffusion derived from the relative translational velocities between neighboring (but not necessarily contacting) particle pairs.

4 Effect of the intermediate principal stress

The conventional measure of soil strength is the friction angle $\phi$, based upon the major and minor principal stresses at failure: $\phi=\sin^{-1}(\sigma_{1}-\sigma_{3})/(\sigma_{1}+\sigma_{3})$. Tests using advanced true-triaxial and hollow ring torsion apparatus demonstrate that strength depends on the intermediate principal stress $\sigma_{2}$, whose relative magnitude is usually represented by the $b$-value, defined as $b=(\sigma_{2}-\sigma_{3})/(\sigma_{1}-\sigma_{3})$. Certain phenomena, however, can obscure the influence of the intermediate stress in laboratory tests. Soils exhibit a propensity for inhomogeneous deformation in the form of shear bands. Shear bands usually appear near the peak stress, and their emergence can alter the subsequent stress-strain behavior and the measured strength. The emergence of shear bands can be either suppressed or promoted by the particular specimen dimensions and boundary conditions, so that the measured influence of the intermediate principal stress is subject to the vagaries of the testing equipment (see Lade 2006 for a review).

DEM simulations were conducted with cubical assemblies measuring about 13.4 particle diameters between periodic boundaries (Section 3.1). Although deformation within an assembly is non-uniform, large-scale localization, such as shear banding, is unable to develop within such limited assembly dimensions. The behavior observed in the simulations can, therefore, be considered close to the underlying material behavior that would prevail during homogeneous deformation.

Figure 11 shows results of these DEM simulations, conducted with different intermediate stresses $\sigma_{2}$.

Ten sets of simulations were conducted: triaxial compression ($b=0$), triaxial extension ($b=1$), and eight sets with intermediate conditions. Each set involved loading the same twenty randomly generated assemblies and averaging the results (see Section 3.1). During loading, mixed boundary conditions were applied in the following manner. A constant lateral stress $\sigma_{3}$ was maintained while the specimens were compressed at a constant rate in the $x_{1}$ direction ($\dot{\epsilon}_{11}=\text{constant}$). Strain control was also maintained in the intermediate direction, with $\epsilon_{22}$ advanced in fixed proportions of $\epsilon_{11}$ (i.e., engineering strain rates $\dot{\epsilon}_{22}=a\dot{\epsilon}_{11}$, with the constant $a=-0.5,$ $-$0.375, $-$0.25, $-$0.125, 0, 0.25, 0.5, 0.75, and 1). Tests were stopped well after the critical state was attained and the $b$-value had become stationary, usually at a compressive engineering strain $-\epsilon_{11}$ of $0.30$ to $0.35$.

Figure 11a shows the failure conditions at the peak and critical states. The principal stresses are plotted on the deviatoric pi-plane and are normalized by dividing by the current mean stress $p$. Each smooth, curved envelope is a spline fit of ten data points. Because the initial fabric was isotropic, an entire failure envelope is symmetrically produced from the ten results that are shown as dots in the first sextant. Flattened Tresca hexagons are also shown, since these envelopes would apply if strength were independent of the intermediate principal stress. Neither the peak nor critical state strengths agree with this idealized condition. The critical state results, however, form a slightly less rotund envelope, lying closer to the Tresca condition.

The results are shown in more detail in Fig. 11b, which gives the friction angles $\phi$ for the ten $b$-values. The figure also shows two commonly used methods for fitting a strength envelope to experimental soil data: the Lade and Matsuoka methods (Lade and Duncan, 1975; Matsuoka, 1976), which have both been calibrated to the triaxial compression strengths ($b=0$). The Lade method closely fits the DEM strength data at the peak state, although neither fitting method captures strength at the critical state.

The micro-mechanical model in Sections 2 and 3 can be used for predicting strength at the critical state. The three Eqs. (6), (17), and (18) are general rate expressions for the densities $\widehat{g}(\mathbf{n})$, $\widehat{\mathbf{f}^{\text{n}}}(\mathbf{n})$, and $\widehat{\mathbf{f}^{\text{t}}}(\mathbf{n})$. At the critical state, the three densities are stationary, and the total rates are zero:

[TABLE]

Each of these differential equations can be expanded and expressed with the functional forms that were developed in Section 3 (e.g. $\dot{\mathbf{n}}(\mathbf{n})$, $(\partial\widehat{f^{\text{n}}}(\mathbf{n})/\partial t)_{\text{matl}}$, etc.). In principle, these three equations can be solved for the scalar density functions $\widehat{g}(\mathbf{n})$, $\widehat{f^{\text{n}}}(\mathbf{n})$, and $\widehat{f^{\text{t}}}(\mathbf{n})$, which then can be used in Eqs. (8) and (3) to find the stress tensor at the critical state. Equations (43) are complex: they are coupled non-linear partial differential equations on the surface of the unit sphere. Looking beyond this difficulty, each of the three equations involve the deformation rate $\mathbf{D}$, which can be considered an input parameter. The solution of the three equations and the resulting stress tensor will, therefore, depend upon the direction of deformation.

Equations (43) were solved for several constant-volume deformation rates $\mathbf{D}$, as would apply at the critical state ($D_{kk}=0$). The coupled non-linear equations were solved using the method of weighted residuals, by approximating each density function as a series of even-powered spherical harmonics, for example,

[TABLE]

where the $\widehat{G}_{\cdots}$ are scalar coefficients, and where the direction index $i_{j}\in\{1,2,3\}$. Approximations for $\widehat{g}(\mathbf{n})$, $\widehat{f^{\text{n}}}(\mathbf{n})$, and $\widehat{f^{\text{t}}}(\mathbf{n})$ were substituted into the appropriate forms of all terms on the right sides of expressions (6), (17), and (18), and the harmonic coefficients were sought to minimize these expressions over the unit sphere:

[TABLE]

thus approximating a simultaneous solution of all three Eqs. (43).

The solution results are shown in Fig. 12. The differential equations were quantified with the particle properties $E$ and $\nu$ and with the micro-mechanical transport properties that were measured at the critical state: $\alpha$, $B^{\text{n}}$, $\beta_{2}^{\text{n}}$, $\gamma_{\text{twirl}}$, $\zeta^{\text{t}}$, and $D_{g}$. These values had been extracted from the single set of DEM simulations of plane-strain biaxial loading ($a=0$), and we simply used the same values in solving the differential equations for other deformation directions $\mathbf{D}$. A common mean stress and average coordination number were imposed as auxiliary constraints ($p=490$ kPa, $2M/N=4.0$). The solutions of Eqs. (43) imply stress tensors that would, in theory, produce critical state flow under various deformation directions $\mathbf{D}$. Figure 12 shows that these solutions compare reasonably well with data from the ten sets of true-triaxial simulation experiments.

5 Conclusion

The Paper presents a model for the evolution of fabric and stress in granular materials. A distinction is made between evolution effects produced by the interactions of particles — termed “material” effects — and effects that are due to an en masse shifting of contacts and forces from one orientation to another. The latter include convection, diffusion, and twirling effects. The DEM simulations show that the material rate has a consistent hardening influence during loading: by itself, the material rate increases the deviatoric stress and induces fabric anisotropy. The other rates have a consistent softening influence, reducing the anisotropies of contact force and contact orientation. The model is most useful in assessing failure at large strains, when the two competing changes are of roughly equal magnitude and both must be tracked during the loading process. At the critical state, the effects are in balance, producing stationary stress and fabric. With this observation, the Paper develops one application of the model: a prediction of the effect of the intermediate principal stress on strength at the critical state.

Although the model’s predictions compare favorably with DEM simulations of small assemblies of 4096 particles, such results may seem questionable when applied to the scale of problems that are encountered in industrial and geotechnical situations. In these large-scale problems, deformation and failure are not homogeneous, but are instead concentrated within shear bands and other localization zones. The Paper’s model is based upon an averaging of trends among the many contacts within a small representative volume, and this averaged behavior can be thought to apply to continuum points or to small regions within a shear band, rather than encompassing the entire band thickness. By pursuing a point-based continuum view, the model might be used to develop advanced, comprehensive continuum models that incorporate a length scale and can predict the emergence and evolution of a shear band. Fulfilling this promise requires work beyond that of the Paper. The effect of the intermediate principal stress, suitably predicted by the model, is but one phenomenon that is known to influence the possibility and nature of shear bands. At least three other phenomena are also known to be relevant, and these phenomena are listed as possible future applications of the model, which might culminate in a more comprehensive continuum description.

As shearing progresses within a granular material, shear bands begin to form before the peak stress is attained, but these early bands are usually transitory and will form, grow, and then disappear. Once the peak stress is reached and the soil begins to soften, the numerous bands coalesce into a few persistent shear bands — perhaps a single shear band — and further deformation becomes concentrated and captive within these persistent features. The theoretical framework in the Paper could be used to model the onset and evolution of softening. The divergence, convection, and diffusion processes reduce the pronounced anisotropies of fabric and force that are attained at the peak state, which might account, at least partially, for the strength softening. 2. 2.

Incremental strains within a shear band are often not aligned with the stress increments: granular materials exhibit a non-coaxiality of the principal strain and stress increments. Such non-coaxiality can favor an abrupt change in the direction of deformation, as at the start of shear banding. The non-coaxiality of stress and strain increments is likely due to induced fabric and force anisotropies, which favor stress increments in certain directions, regardless of the direction of deformation. The model is naturally suited to account for such effects of fabric anisotropy. 3. 3.

Shear bands have a characteristic thickness that is related to particle size. A comprehensive and verifiable explanation of this characteristic remains an open problem in granular mechanics. Several continuum theories have been proposed as a rationale for band thickness, but to the author’s knowledge, only one such explanation has been confirmed: a dependence of the stress increment upon the spatial gradients of strain (Kuhn, 2005). The model in the Paper can, perhaps, be extended to understand and explain the effect of strain gradients, since such gradients will alter the contact migration pattern of Eqs. (23) and (24), changing the evolution of fabric and stress.

An application of the model to these three phenomena will likely require an understanding of density rates during unloading, a matter not addressed in the Paper. Although extending the model to include unloading rates presents a further challenge, once achieved, the model might also provide a micro-mechanical alternative to the notion of a yield surface.