Effect of particle friction and polydispersity on the macroscopic stress–strain relations of granular materials

Abstract

The macroscopic mechanical behavior of granular materials inherently depends on the properties of particles that compose them. Using the discrete element method, the effect of particle contact friction and polydispersity on the macroscopic stress response of 3D sphere packings is studied. The analytical expressions for the pressure, coordination number and fraction of rattlers proposed for isotropically deformed frictionless systems also hold when the interparticle coefficient of friction is finite; however, the numerical values of the parameters such as the jamming volume fraction change with varying microscopic contact and particle properties. The macroscopic response under deviatoric loading is studied with triaxial test simulations. Concerning the shear strength, our results agree with previous studies showing that the deviatoric stress ratio increases with particle coefficient of friction μ starting from a nonzero value for μ = 0 and saturating for large μ. On the other hand, the volumetric strain does not have a monotonic dependence on the particle contact friction. Most notably, maximum compaction is reached at an intermediate value of the coefficient of friction μ ≈ 0.3. The effect of polydispersity on the macroscopic stress–strain relationship cannot be studied independent of initial packing conditions. The shear strength increases with polydispersity when the initial volume fraction is fixed, but the effect of polydispersity is much less pronounced when the initial pressure of the packings is fixed. Finally, a simple hypoplastic constitutive model is calibrated with numerical test results following an established procedure to ascertain the relation between particle properties and material coefficients of the macroscopic model. The calibrated model is in good qualitative agreement with simulation results.

Appendix: Calibration of the hypoplastic constitutive model

The hypoplastic constitutive model given in Eq. (4) can be calibrated for a specific material with the results of a triaxial test [15].

Due to the simple geometry of the test setup, the stress and strain rate tensors are characterized by their principal components:

$$ -{\mathbf{T}}=\left(\begin{array}{ccc}\sigma_1 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 &\sigma_3 \end{array} \right) , \quad {\mathbf{D}}=\left(\begin{array}{ccc} \dot{\varepsilon_1} & 0 & 0 \\ 0 & \dot{\varepsilon_2} & 0\\ 0 & 0 &\dot{\varepsilon_3}\end{array}\right), $$

(5)

where compressive stresses are positive.

As illustrated in Fig. 14, the values of (σ ₁ − σ ₂)_max, the slope ³ E and the angles β _A and β _B at points A and B can be computed from the test results and are related to $\mathbf{T}$ and $\mathbf{D}$:

$$ \beta_{A/B}=\arctan\left(\frac{\dot{\varepsilon}_\mathrm{v}} {\dot{\varepsilon}_1}\right)_{\text{A/B}} =\arctan\left(\frac{\dot{\varepsilon}_1+2\dot{\varepsilon}_2}{ \dot{\varepsilon}_1}\right)_{\text{A/B}} =\arctan\left(1+2\frac{\dot{\varepsilon}_2}{ \dot{\varepsilon}_1}\right)_{\text{A/B}} $$

(6)

Fig. 14

Schematic representation of a triaxial test result for the calibration of the hypoplastic constitutive model

Since the hypoplastic constitutive model is rate independent, the magnitude of the strain rate $|\dot{\varepsilon}_1|$ can be arbitrary. However, the sign of $\dot{\varepsilon}_1$ must be negative due to compression during a conventional triaxial test. Therefore, for simplicity $\dot{\varepsilon}_1=-1$ is chosen, so that the strain rate tensor $\mathbf{D}$ at points A and B is:

$$ {\mathbf{D}}_{A/B}=\left(\begin{array}{ccc} -1 & 0 & 0 \\ 0 & \frac{1}{2}(1-\tan\beta_{A/B}) & 0 \\ 0 & 0 & \frac{1}{2}(1-\tan\beta_{A/B}) \end{array}\right). $$

(7)

The stress tensor $\mathbf{T}$ at points A and B is known:

$$ - {\mathbf{T}}_A=\left(\begin{array}{ccc} \sigma_2 & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 &\sigma_2 \end{array} \right) \quad\hbox{ and }\, - {\mathbf{T}}_B=\left(\begin{array}{ccc} \sigma_2-(\sigma_1-\sigma_2)_{\rm max} & 0 & 0 \\ 0 & \sigma_2 & 0 \\ 0 & 0 &\sigma_2 \end{array}\right), $$

(8)

and the stress rates are given by:

$$ \dot{{\mathbf{T}}}_A=\left(\begin{array}{ccc}- E & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right) , \quad \dot{{\mathbf{T}}}_B=\left(\begin{array}{ccc}0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}\right), $$

(9)

where at the point $A, \dot{\sigma}_1=E\dot{\varepsilon}_1$ since $\dot{\sigma_2}=0$ and $\dot{\varepsilon}_1=-1$.

Substituting $\mathbf{D},\mathbf{T}\hbox{ and }\dot{\mathbf{T}}$ computed at points A and B into Eq. (4), the following system of equations is obtained with the unknowns C ₁, C ₂, C ₃ and C ₄:

$$ \left( \begin{array}{cccc} \hbox{tr}{\mathbf{T}}_AD^A_{1,1} & T^A_{1,1}\frac{\hbox{tr}({\mathbf{T}}{\mathbf{D}})_A}{\hbox{tr}{\mathbf{T}}_A} &||{\mathbf{D}}_A||\frac{(T_{1,1})_A^2}{\hbox{tr}{\mathbf{T}}_A} &||{\mathbf{D}}_A||\frac{(T^{*}_{1,1})_A^2}{\hbox{tr}{\mathbf{T}}_A} \\ \hbox{tr}{\mathbf{T}}_AD^A_{2,2} & T^A_{2,2}\frac{\hbox{tr}({\mathbf{T}}{\mathbf{D}})_A}{\hbox{tr}{\mathbf{T}}_A} &||{\mathbf{D}}_A||\frac{(T_{2,2})_A^2}{\hbox{tr}{\mathbf{T}}_A} &||{\mathbf{D}}_A||\frac{(T^{*}_{2,2})_A^2}{\hbox{tr}{\mathbf{T}}_A}\\ \hbox{tr}{\mathbf{T}}_BD^B_{1,1} &T^B_{1,1}\frac{\hbox{tr}({\mathbf{T}} {\mathbf{D}})_B}{\hbox{tr}{\mathbf{T}}_B}& ||{\mathbf{D}}_B||\frac{(T_{1,1})_B^2}{\hbox{tr}{\mathbf{T}}_B} & ||{\mathbf{D}}_B||\frac{(T^{*}_{1,1})_B^2}{\hbox{tr}{\mathbf{T}}_B} \\ \hbox{tr}{\mathbf{T}}_BD^B_{2,2} &T^B_{2,2}\frac{\hbox{tr} ({\mathbf{T}}{\mathbf{D}})_B}{\hbox{tr}{\mathbf{T}}_B}& ||{\mathbf{D}}_B||\frac{(T_{2,2})_B^2}{\hbox{tr}{\mathbf{T}}_B} & ||{\mathbf{D}}_B||\frac{(T^{*}_{2,2})_B^2}{\hbox{tr}{\mathbf{T}}_B} \end{array} \right) \left( \begin{array}{c} C_1 \\ C_2 \\ C_3 \\ C_4 \\ \end{array} \right) = \left( \begin{array}{c} - E \\ 0 \\ 0 \\ 0 \\ \end{array} \right) $$

(10)

where for clarity, the letters A and B have been switched to superscripts when the indicial notation of the tensors is used. The solution to (10) can be obtained by simple matrix inversion using linear algebra or well-known numerical methods such as Gauss–Seidel.