Granular material

Chain of transmission of stress forces in a granular medium

Coulomb regarded internal forces between granular particles as a friction process, and proposed the friction law, that the force of friction of solid particles is proportional to the normal pressure between them and the static friction coefficient is greater than the kinetic friction coefficient. He studied the collapse of piles of sand and found empirically two critical angles: the maximal stable angle ${\displaystyle \theta _{m}}$ and the minimum angle of repose ${\displaystyle \theta _{r}}$. When the sandpile slope reaches the maximum stable angle, the sand particles on the surface of the pile begin to fall. The process stops when the surface inclination angle is equal to the angle of repose. The difference between these two angles, ${\displaystyle \Delta \theta =\theta _{m}-\theta _{r}}$, is the Bagnold angle, which is a measure of the hysteresis of granular materials. This phenomenon is due to the force chains: stress in a granular solid is not distributed uniformly but is conducted away along so-called force chains which are networks of grains resting on one another. Between these chains are regions of low stress whose grains are shielded for the effects of the grains above by vaulting and arching. When the shear stress reaches a certain value, the force chains can break and the particles at the end of the chains on the surface begin to slide. Then, new force chains form until the shear stress is less than the critical value, and so the sandpile maintains a constant angle of repose.[7]

In 1895, H. A. Janssen discovered that in a vertical cylinder filled with particles, the pressure measured at the base of the cylinder does not depend on the height of the filling, unlike Newtonian fluids at rest which follow Stevin's law. Janssen suggested a simplified model with the following assumptions:

1) The vertical pressure, ${\displaystyle \sigma _{zz}}$, is constant in the horizontal plane;

2) The horizontal pressure, ${\displaystyle \sigma _{rr}}$, is proportional to the vertical pressure ${\displaystyle \sigma _{zz}}$, where ${\displaystyle K={\frac {\sigma _{rr}}{\sigma _{zz}}}}$ is constant in space;

3) The wall friction static coefficient ${\displaystyle \mu ={\frac {\sigma _{rz}}{\sigma _{rr}}}}$ sustains the vertical load at the contact with the wall;

4) The density of the material is constant over all depths.

The pressure in the granular material is then described in a different law, which accounts for saturation:

$${\displaystyle p(z)=p_{\infty }[1-\exp(-z/\lambda )]}$$

where ${\displaystyle \lambda ={\frac {R}{2\mu K}}}$ and ${\displaystyle R}$ is the radius of the cylinder, and at the top of the silo ${\displaystyle z=0}$.

The given pressure equation does not account for boundary conditions, such as the ratio between the particle size to the radius of the silo. Since the internal stress of the material cannot be measured, Janssen's speculations have not been verified by any direct experiment.

In the early 1960s, Rowe studied dilatancy effect on shear strength in shear tests and proposed a relation between them.

The mechanical properties of assembly of mono-dispersed particles in 2D can be analyzed based on the representative elementary volume, with typical lengths, ${\displaystyle \ell _{1},\ell _{2}}$, in vertical and horizontal directions respectively. The geometric characteristics of the system is described by ${\displaystyle \alpha =\arctan({\frac {\ell _{1}}{\ell _{2}}})}$ and the variable ${\displaystyle \beta }$, which describes the angle when the contact points begin the process of sliding. Denote by ${\displaystyle \sigma _{11}}$ the vertical direction, which is the direction of the major principal stress, and by ${\displaystyle \sigma _{22}}$ the horizontal direction, which is the direction of the minor principal stress.

Then stress on the boundary can be expressed as the concentrated force borne by individual particles. Under biaxial loading with uniform stress ${\displaystyle \sigma _{12}=\sigma _{21}=0}$ and therefore ${\displaystyle F_{12}=F_{21}=0}$.

At equilibrium state:

$${\displaystyle {\frac {F_{11}}{F_{22}}}={\frac {\sigma _{11}\ell _{2}}{\sigma _{22}\ell _{1}}}=\tan(\theta +\beta )}$$

where ${\displaystyle \theta }$ , the friction angle, is the angle between the contact force and the contact normal direction.

${\displaystyle \theta _{\mu }}$, which describes the angle that if the tangential force falls within the friction cone the particles would still remain steady. It is determined by the coefficient of friction ${\displaystyle \mu =tg\phi _{u}}$, so ${\displaystyle \theta \leq \theta _{\mu }}$. Once stress is applied to the system then ${\displaystyle \theta }$ gradually increases while ${\displaystyle \alpha ,\beta }$ remains unchanged. When ${\displaystyle \theta \geq \theta _{\mu }}$ then the particles will begin sliding, resulting in changing the structure of the system and creating new force chains. ${\displaystyle \Delta _{1},\Delta _{2}}$,the horizontal and vertical displacements respectively satisfies:

$${\displaystyle {\frac {\dot {\Delta _{2}}}{\dot {\Delta _{1}}}}={\frac {{\dot {\varepsilon _{22}}}\ell _{2}}{{\dot {\varepsilon _{11}}}\ell _{1}}}=-\tan \beta }$$

Consider ${\displaystyle N}$ particles, particle ${\displaystyle i}$ having energy ${\displaystyle \varepsilon _{i}}$. At some constant rate per unit time, randomly choose two particles ${\displaystyle i,j}$ with energies ${\displaystyle \varepsilon _{i},\varepsilon _{j}}$ and compute the sum ${\displaystyle \varepsilon _{i}+\varepsilon _{j}}$. Now, randomly distribute the total energy between the two particles: choose randomly ${\displaystyle z\in \left[0,1\right]}$ so that the first particle, after the collision, has energy ${\displaystyle z\left(\varepsilon _{i}+\varepsilon _{j}\right)}$, and the second ${\displaystyle \left(1-z\right)\left(\varepsilon _{i}+\varepsilon _{j}\right)}$.

The stochastic evolution equation:

$${\displaystyle \varepsilon _{i}(t+dt)={\begin{cases}\varepsilon _{i}(t)&probability:\,1-\Gamma dt\\z\left(\varepsilon _{i}(t)+\varepsilon _{j}(t)\right)&probability:\,\Gamma dt\end{cases}}}$$

where ${\displaystyle \Gamma }$ is the collision rate, ${\displaystyle z}$ is randomly picked from ${\displaystyle \left[0,1\right]}$ (uniform distribution) and j is an index also randomly chosen from a uniform distribution. The average energy per particle:

$${\displaystyle {\begin{aligned}\left\langle \varepsilon (t+dt)\right\rangle &=\left(1-\Gamma dt\right)\left\langle \varepsilon (t)\right\rangle +\Gamma dt\cdot \left\langle z\right\rangle \left(\left\langle \varepsilon _{i}\right\rangle +\left\langle \varepsilon _{j}\right\rangle \right)\\&=\left(1-\Gamma dt\right)\left\langle \varepsilon (t)\right\rangle +\Gamma dt\cdot {\dfrac {1}{2}}\left(\left\langle \varepsilon (t)\right\rangle +\left\langle \varepsilon (t)\right\rangle \right)\\&=\left\langle \varepsilon (t)\right\rangle \end{aligned}}}$$

The second moment:

${\displaystyle {\begin{aligned}\left\langle \varepsilon ^{2}(t+dt)\right\rangle &=\left(1-\Gamma dt\right)\left\langle \varepsilon ^{2}(t)\right\rangle +\Gamma dt\cdot \left\langle z^{2}\right\rangle \left\langle \varepsilon _{i}^{2}+2\varepsilon _{i}\varepsilon _{j}+\varepsilon _{j}^{2}\right\rangle \\&=\left(1-\Gamma dt\right)\left\langle \varepsilon ^{2}(t)\right\rangle +\Gamma dt\cdot {\dfrac {1}{3}}\left(2\left\langle \varepsilon ^{2}(t)\right\rangle +2\left\langle \varepsilon (t)\right\rangle ^{2}\right)\end{aligned}}}$

Now the time derivative of the second moment:

${\displaystyle {\dfrac {d\left\langle \varepsilon ^{2}\right\rangle }{dt}}=lim_{dt\rightarrow 0}{\dfrac {\left\langle \varepsilon ^{2}(t+dt)\right\rangle -\left\langle \varepsilon ^{2}(t)\right\rangle }{dt}}=-{\dfrac {\Gamma }{3}}\left\langle \varepsilon ^{2}\right\rangle +{\dfrac {2\Gamma }{3}}\left\langle \varepsilon \right\rangle ^{2}}$

In steady state:

${\displaystyle {\dfrac {d\left\langle \varepsilon ^{2}\right\rangle }{dt}}=0\Rightarrow \left\langle \varepsilon ^{2}\right\rangle =2\left\langle \varepsilon \right\rangle ^{2}}$

Solving the differential equation for the second moment:

${\displaystyle \left\langle \varepsilon ^{2}\right\rangle -2\left\langle \varepsilon \right\rangle ^{2}=\left(\left\langle \varepsilon ^{2}(0)\right\rangle -2\left\langle \varepsilon (0)\right\rangle ^{2}\right)e^{-{\frac {\Gamma }{3}}t}}$

However, instead of characterizing the moments, we can analytically solve the energy distribution, from the moment generating function. Consider the Laplace transform: ${\displaystyle g(\lambda )=\left\langle e^{-\lambda \varepsilon }\right\rangle =\int _{0}^{\infty }e^{-\lambda \varepsilon }\rho (\varepsilon )d\varepsilon }$.

Where ${\displaystyle g(0)=1}$, and ${\displaystyle {\dfrac {dg}{d\lambda }}=-\int _{0}^{\infty }\varepsilon e^{-\lambda \varepsilon }\rho (\varepsilon )d\varepsilon =-\left\langle \varepsilon \right\rangle }$

the n derivative:

${\displaystyle {\dfrac {d^{n}g}{d\lambda ^{n}}}=\left(-1\right)^{n}\int _{0}^{\infty }\varepsilon ^{n}e^{-\lambda \varepsilon }\rho (\varepsilon )d\varepsilon =\left\langle \varepsilon ^{n}\right\rangle }$

now:

${\displaystyle e^{-\lambda \varepsilon _{i}(t+dt)}={\begin{cases}e^{-\lambda \varepsilon _{i}(t)}&1-\Gamma t\\e^{-\lambda z\left(\varepsilon _{i}(t)+\varepsilon _{j}(t)\right)}&\Gamma t\end{cases}}}$

${\displaystyle \left\langle e^{-\lambda \varepsilon \left(t+dt\right)}\right\rangle =\left(1-\Gamma dt\right)\left\langle e^{-\lambda \varepsilon _{i}(t)}\right\rangle +\Gamma dt\left\langle e^{-\lambda z\left(\varepsilon _{i}(t)+\varepsilon _{j}(t)\right)}\right\rangle }$

${\displaystyle g\left(\lambda ,t+dt\right)=\left(1-\Gamma dt\right)g\left(\lambda ,t\right)+\Gamma dt\int _{0}^{1}{\underset {=g^{2}(\lambda z,t)}{\underbrace {\left\langle e^{-\lambda z\varepsilon _{i}(t)}\right\rangle \left\langle e^{-\lambda z\varepsilon _{j}(t)}\right\rangle } }}dz}$

Solving for ${\displaystyle g(\lambda )}$ with change of variables ${\displaystyle \delta =\lambda z}$:

${\displaystyle \lambda g(\lambda )=\int _{0}^{\lambda }g^{2}(\delta )d\delta \Rightarrow \lambda g'(\lambda )+g(\lambda )=g^{2}(\lambda )\Rightarrow g(\lambda )={\dfrac {1}{\lambda T+1}}}$

We will show that ${\displaystyle \rho (\varepsilon )={\dfrac {1}{T}}e^{-{\frac {\varepsilon }{T}}}}$ (Boltzmann Distribution) by taking its Laplace transform and calculate the generating function:

${\displaystyle \int _{0}^{\infty }{\dfrac {1}{T}}e^{-{\frac {\varepsilon }{T}}}\cdot e^{-\lambda \varepsilon }d\varepsilon ={\dfrac {1}{T}}\int _{0}^{\infty }e^{-\left(\lambda +{\frac {1}{T}}\right)\varepsilon }d\varepsilon =-{\dfrac {1}{T\left(\lambda +{\frac {1}{T}}\right)}}e^{-\left(\lambda +{\frac {1}{T}}\right)\varepsilon }|_{0}^{\infty }={\dfrac {1}{\lambda T+1}}=g(\lambda )}$

Sand dunes

Some beach sands, such as those of the aptly named Squeaky Beach, exhibit squeaking when walked upon. Some desert dunes are known to exhibit booming during avalanching or when their surface is otherwise disturbed. Granular materials discharged from silos produce loud acoustic emissions in a process known as silo honking.