Micro-macro correlations and anisotropy in granular assemblies under uniaxial loading and unloading

The influence of contact friction on the behavior of dense, polydisperse granular assemblies under uniaxial (oedometric) loading and unloading deformation is studied using discrete element simulations. Even though the uniaxial deformation protocol is one of the “simplest” element tests possible, the evolution of the structural anisotropy necessitates its careful analysis and understanding, since it is the source of interesting and unexpected observations. On the macroscopic, homogenized, continuum scale, the deviatoric stress ratio and the deviatoric fabric, i.e., the microstructure behave in a different fashion during uniaxial loading and unloading. The maximal stress ratio and strain increase with increasing contact friction. In contrast, the deviatoric fabric reaches its maximum at a unique strain level independent of friction, with the maximal value decreasing with friction. For unloading, both stress and fabric respond to unloading strain with a friction-dependent delay but at different strains. On the micro-level, a friction-dependent non-symmetry of the proportion of weak (strong) and sliding (sticking) contacts with respect to the total contacts during loading and unloading is observed. Coupled to this, from the directional probability distribution, the “memory” and history-dependent behavior of granular systems is confirmed. Surprisingly, while a rank-2 tensor is sufficient to describe the evolution of the normal force directions, a sixth order harmonic approximation is necessary to describe the probability distribution of contacts, tangential force, and mobilized friction. We conclude that the simple uniaxial deformation activates microscopic phenomena not only in the active Cartesian directions, but also at intermediate orientations, with the tilt angle being dependent on friction, so that this microstructural features cause the interesting, nontrivial macroscopic behavior.

Figure 1

Evolution of volume fraction as a function of time. Region A represents the initial isotropic compression below the jamming volume fraction. B represents relaxation of the system to fully dissipate the systems potential and kinetic energy and C represents the subsequent isotropic compression up to ${\nu }_{\max }=0.820$ and then subsequent decompression. Cyan dots represent some of the initial configurations, at different ${\nu }_i$, during the loading cycle; blue stars, for the same ${\nu }_i$, are different configurations, since obtained during the unloading cycle; both can be chosen for further study.

Figure 2

The nondimensional pressure plotted as function of volume fraction under uniaxial deformation for different friction coefficients during (a) loading and (b) unloading and coordination number (excluding rattlers) for the same dataset during (a) loading and (b) unloading. The vertical arrows show increasing (and decreasing) $\mu$ while the tilted arrows show the loading and unloading direction.

Figure 3

The deviatoric stress ratio plotted as a function of deviatoric strain during uniaxial (a) loading and (b) unloading. The corresponding plots of the deviatoric fabric plotted during uniaxial (c) loading and (d) unloading, for different microscopic friction coefficients.

Figure 4

Trend of the peak deviatoric stress and peak deviatoric fabric with increasing microscopic friction coefficient $\mu$ under uniaxial loading, given a maximal strain ${\varepsilon }_{\mathrm{dev}}^{\max }=0.1549.$ $s_{\mathrm{dev}}^{\max }$ values for $\mu >0.1$ are taken at ${\varepsilon }_{\mathrm{dev}}^{\max }$ since no clear maximum is achieved. Dashed line indicates ${\mu }^{\mathrm{macro}}=\mu$.

Figure 5

Eigenvalues of deviatoric stress for $\mu =0.1$ plotted as functions of the deviatoric strain for (a) loading and (b) unloading along with their corresponding orientations with respect to the compressive direction during uniaxial (c) loading and (d) unloading.

Figure 6

Eigenvalues of the deviatoric fabric for $\mu =0.1$ plotted as functions of the deviatoric strain for (a) loading and (b) unloading along with their corresponding orientations with respect to the compressive direction during uniaxial (c) loading and (d) unloading. Dotted lines indicated the slope of the path, identical before and after strain reversal.

Figure 7

Shape factors of (a) stress and (b) fabric as function of the deviatoric strain for some exemplary friction coefficients.

Figure 8

Orientation of the largest positive (a) stress eigenvector for all contacts, (b) fabric eigenvector for all contacts, (c) stress eigenvector for strong contacts, (d) fabric eigenvector for strong contacts, (e) stress eigenvector for weak contacts, (f) fabric eigenvector for weak contacts plotted against dimensionless time for different coefficient of friction.

Figure 9

Strains at which the orientations of the stress eigenvectors cross $\theta =$ ${45}^{\circ }$ and at which the deviatoric stress ratio, deviatoric fabric, and the stress shape factor cross zero for frictions $\mu =0$, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, and 1.0.

Figure 10

Tracking $f_t/\mu f_n$ and $f^{*}=f_n/\langle f_n\rangle$ for two single particle pairs randomly selected from the ensemble during compression and decompression where $\tau$ is the dimensionless time. (a) Particle pair 1. (b) Particle pair 2.

Figure 11

Proportion of (a) weak contacts (b) sliding contacts with respect to the total number of contacts during uniaxial loading and unloading cycle for different friction coefficients.

Figure 12

Fraction of (a) weak sliding contacts ($wsl$) and (b) strong sticking contacts ($sst$) with respect to the total number of weak $(\sum w)$ and strong $(\sum s)$ contacts, respectively, during uniaxial loading and unloading for different friction coefficients.

Figure 13

Fraction of (a) weak sliding contacts ($wsl$) and (b) strong sticking contacts ($sst$) with respect to the total number of sliding $(\sum sl)$ and sticking $(\sum st)$ contacts, respectively, during uniaxial loading and unloading for different friction coefficients.

Figure 14

Normalized probability density of the normal force $P(f/\langle f_{\mathrm{all}}\rangle )$ for the three reference directions and for all contact forces plotted against the normalized force $f_{*}=f/\langle f_{\mathrm{all}}\rangle$ for $\mu =0.1$ and cutoff $\chi =0.8$. Three snapshots are shown at (a) initial state $\tau =0$, (b) maximum, $\tau =0.5$, and (c) final state $\tau =1.0$.

Figure 15

Comparison of the sixth order anisotropy state descriptors (a) ${\epsilon }_2$, (b) ${\epsilon }_4$, (c) ${\epsilon }_6$ as a function of the deviatoric strain for different friction coefficients during loading (left) and unloading (right).

Figure 16

Snapshots of the polar distribution of the (a) contacts $P_6\left(\theta \right)$, (b) normal force $f_n\left(\theta \right)$, (c) tangential force $f_t\left(\theta \right)$, (d) mobilized friction ${\psi }_t\left(\theta \right)$ at different dimensionless time ($\tau$) during uniaxial loading and unloading for friction $\mu =0.1$.

Figure 17

Polar distribution of the (a) contacts $P_6\left(\theta \right)$, (b) normal force $f_n\left(\theta \right)$, (c) tangential force $f_t\left(\theta \right)$, (d) mobilized friction ${\psi }_t\left(\theta \right)$ at dimensionless time ($\tau =0.5$) for friction $\mu =0.1$.

Figure 18

Schematic representation of the angles of the contact unit vector (green arrow) for (a) the constant bin width $\Delta \theta$ method and (b) the constant height method. The angles $\theta$ and $\phi$ are the azimuthal angle and polar angle of the system, respectively.

Figure 19

Sixth order fits of methods 1, 2, and 3 to the $P\left(\theta \right)=({\varphi }^{\theta }/\Delta \theta sin\theta )$ data at dimensionless time $\tau =0.076$. The solid red symbols represent data obtained using the constant height method ($h$), while the triangles are those obtained with the constant bin width method ($b$).

Figure 20

Evolution of the anisotropy state descriptors ${\epsilon }_2$, ${\epsilon }_4$, and ${\epsilon }_6$ of the sixth order expansion Eq. (25) as function of deviatoric strain using the three methods for $\mu =0$.