Discrete Element Method Investigation of Binary Granular Flows with Different Particle Shapes

Abstract

The effects of particle shape differences on binary mixture shear flows are investigated using the Discrete Element Method (DEM). The binary mixtures consist of frictionless rods and disks, which have the same volume but significantly different shapes. In the shear flows, stacking structures of rods and disks are formed. The effects of the composition of the mixture on the stacking are examined. It is found that the number fraction of stacking particles is smaller for the mixtures than for the monodisperse rods and disks. For binary mixtures with small particle shape differences, the mixture stresses are bounded by the stresses of the two corresponding monodisperse systems. However, for binary mixtures with large particle shape differences, the stresses of the mixtures can be larger than the stresses of the monodisperse systems at large solid volume fractions because larger differences in particle shape cause geometrical interference in packing, leading to stronger particle–particle interactions in the flow. The stresses in dense binary mixtures are found to be exponential functions of the order parameter, which is a measure of particle alignment. Based on the simulation results, an empirical expression for the bulk friction coefficient (ratio of the shear stress to normal stress) for dense binary flows is proposed by accounting for the effects of particle alignment and solid volume fraction.

1. Introduction

Granular shear flows are prevalent in natural and industrial processes, such as snow avalanches, landslides, sand storms, mineral transport, chemical and pharmaceutical processing, coal combustion, and food production. Studies on the mechanics and physics of granular flows have been carried out because a good understanding of the granular flow properties is useful to control damage caused by natural disasters and to optimize industrial processes. Theories of granular flow rheology have been developed. By treating particles as gas molecules, the kinetic theory of gases was introduced to statistically describe mechanical and thermal properties of dilute, dry granular flows [1], and the theory was thereafter extended to the dense flow regime [2,3]. Jop et al. [4] found that the ratio of shear to normal stress, $\mu$, was a function of the inertial number, $I=\gamma d_{\mathrm{p}}/\sqrt{p/{\rho }_{\mathrm{p}}}$, in which $\gamma$ is local shear rate, $p$ is local pressure, and $d_{\mathrm{p}}$ and ${\rho }_{\mathrm{p}}$ are particle diameter and density, respectively. Thus, the constitutive relation for dense flows is described by the correlation $\mu =\mu \left(I\right)$, which is well known as the $\mu \left(I\right)$-theory. Non-local effects have been considered in subsequent work [5,6]. Recently, Pähtz et al. [7] observed that the ratio of shear to normal stress, $\mu$, scaled with the square root of the local Péclet number, $\sqrt{P_{\mathrm{e}}}$, for dry, wet, dilute, and dense granular flows, provided that the particle diameter exceeds the particle mean free path.

Most granular materials are mixtures containing particles of different sizes, densities, shapes, and elastic moduli. The flow behavior of mixtures can be distinct from that of monodisperse systems. Shen [8] derived a set of constitutive equations for the stresses in a rapid flow of spherical particles with two sizes based on the assumptions: (i) binary collisions are the major stress-generating mechanism; (ii) stresses are caused by momentum exchange between surface particles and exterior particles; (iii) the work done by shear stresses on the boundaries is balanced by the energy dissipation through particle–particle collisions; and (iv) a larger particle and a smaller particle have the same random fluctuation energy, i.e., equipartition of energy. The predictions by Shen’s theory was good for a diameter ratio close to unity and deteriorated for large particle size differences, as pointed out in [9]. Jenkins and Mancini [9] and Farrell et al. [10] extended the kinetic theory, which was developed for monodisperse systems [1], to describe granular flows of binary mixtures. In the work of Farrell et al. [10], equipartition of energy between different particle species was also assumed and a radial distribution function at contact was defined as a function of the sizes and solid concentrations of the two particle species. Compared with experimental results, the prediction of Farrell et al. [10] was in good agreement at small solid volume fractions, but significant discrepancies at large solid volume fractions were observed due to the dominance of enduring contacts, which violated the binary collision assumption. In the kinetic theory of Jenkins and Mancini [9], it was assumed that the pair distribution function for two 2D disks is the product of Maxwellian velocity distributions and a factor that accounts for the effects of excluded area and collisional shielding, and non-equipartition of energy is considered. Compared to Shen’s work [8], the theory of Jenkins and Mancini can correctly predict that the stresses decrease with an increase in small particle concentration over a much wider range of particle diameter ratios. By incorporating the Revised Enskog Theory [11] and using non-Maxwellian velocity distributions, Willits and Arnarson [12] developed another version of kinetic theory, which showed better agreement with the discrete particle simulation results than the Jenkins-Mancini theory [9]. Shear flow simulations of binary mixtures [13] showed that equipartition of energy between large and small particles rapidly deteriorated as the coefficient of restitution decreased and the particle size difference increased. Galvin et al. [14] evaluated the effect of non-equipartition of energy on rapid flows of granular mixtures. It was found that for simple shear flows without segregation, the prediction from kinetic theory was insensitive to the equipartition versus non-equipartition treatment, while for segregating flows the presence of non-equipartition gave rise to additional components of the driving forces for the size segregation, resulting in better predictions. Thus, non-equipartition of energy should be considered in kinetic theory for describing a wide range of polydisperse granular flows. Iddir and Arastoopour [15] considered non-Maxwellian velocity distributions and non-equipartition in their kinetic theory. By comparing four different kinetic theories of granular mixtures, Benyahia [16] demonstrated that the Iddir–Arastoopour theory gave the best prediction of the discrete particle simulations.

Dense flows of granular mixtures were computationally investigated using the Discrete Element Method (DEM) by Tripathi and Khakhar [17] and Gu et al. [18]. It was found that like monodisperse dense flows, polydisperse dense flows also followed a $\mu \left(I\right)$-like rheological correlation, in which the shear to normal stress ratio $\mu$ is a function of the inertial number $I$. However, the yielding stress ratio in the $\mu \left(I\right)$ model depended on the polydispersity and skewness of the particle size distribution [18]. Dahl et al. [19,20] explored the effects of continuous size distributions on rapid granular flows using numerical simulations of simple shear flows. It was observed that the stresses with Gaussian and lognormal size distributions were consistent with the stresses of monodisperse systems when the root-mean-square (2D) or root-mean-cubed (3D) particle diameter of the mixture was equal to the diameter of the monodisperse system. A similar conclusion was also obtained from the kinetic theory of Farrell et al. [10].

Previous work of granular mixtures mainly focused on differences in particle size and density, and the effects of particle shape differences in a mixture were much less studied. Yang et al. [21] performed DEM simulations of shear flows of binary mixtures with large oval and rod-like particles and small spheres. They found that the particle projected area in the plane perpendicular to the flow direction played an important role. The stress–solid volume fraction curves of binary mixtures with different volume ratios of small-to-large particle species collapsed on a single master curve by normalizing the stresses using the effective particle projected area and root-mean-cubed diameter. Shear flows of mixtures of cylindrical particles were studied using DEM by Hao et al. [22]. In these mixtures, the cylindrical particles had the same diameter but different lengths. It was observed that the interaction between different particle species increased the alignment of shorter particles and reduced the alignment of longer particles. The total stress tensor of a mixture could be expressed as a sum of the stress tensors of the monodispersed particle species, which were linearly scaled by the concentrations of the corresponding species.

For an improved understanding of the effects of particle shape differences, shear flows of binary mixtures of elongated rods and flat disks are simulated using DEM in the present work. In these binary mixtures, a rod has the same volume as a disk, and thus all the particles have the same equivalent volume diameter $d_{\mathrm{v}}$. Therefore, the effects of particle size differences are eliminated and only particle shape differences exist in the mixtures. In the shear flows, the rods and disks exhibit significant stacking and alignment. The effects of the composition in a mixture on the stacking and alignment behaviors and the particle-phase stresses are then explored. Lastly, a correlation between the microstructure and the stresses is investigated.

2. Methodology

2.1. Cylindrical Particle Discrete Element Method

The Discrete Element Method (DEM), which was first proposed by Cundall and Strack [23], is a popular numerical tool for studying particulate systems. In DEM, the motion of each particle is incrementally solved by time integration, and the contact forces between the particles need to be determined in each time step. Based on the algorithms of contact detection for cylindrical objects presented by Kodam et al. [24] and Guo et al. [25], an in-house DEM code is developed and utilized in this work for the modeling of cylindrical particles. In this method, the translational and rotational motion of an individual particle i is described based on Newton’s law of motion,

$$m_{\mathrm{i}}\frac{d{\mathit{v}}_{\mathrm{i}}}{dt}={\mathit{F}}_{\mathrm{ci}}$$

(1)

$${\mathit{I}}_{\mathrm{i}}\frac{d{\mathit{\omega }}_{\mathrm{i}}}{dt}-\left({\mathit{I}}_{\mathrm{i}}·{\mathit{\omega }}_{\mathrm{i}}\right)\times {\mathit{\omega }}_{\mathrm{i}}={\mathit{T}}_{\mathrm{i}}$$

(2)

where v_i and ω_i are the translational and angular velocities of the particle i, respectively. The translational movement of the particle of the mass m_i is driven by the contact force ${\mathit{F}}_{\mathrm{ci}}$, and the rotation is induced by the torque T_i generated from the contact forces (including both normal and tangential forces for a cylindrical particle). For 3D non-spherical particles, the parameter I_i is the moment of inertia tensor. The velocities and positions of the particles can be updated by the time integration of Equations (1) and (2) with a fixed time step. The orientation of a cylindrical particle is described by a rotation matrix and updated by the current angular velocity at each time step [26].

The geometries of cylindrical particles in the DEM simulations are illustrated in Figure 1. The aspect ratio (AR) of a cylinder is defined as the ratio of particle length to diameter, i.e., $AR=l/d$. Thus, a rod and a disk correspond to AR > 1 and AR < 1, respectively. Six typical contact types exist for the contact between two cylindrical particles [24]: face–face, face–band, face–edge, band–band (parallel and skewed), band–edge, and edge–edge. In addition, several special contact types were recognized [25]. The algorithms to determine overlap size, contact normal direction, and contact point position for each contact type, which are provided in [24,25], are summarized for convenience in Appendix A. The overlap size ${\delta }_{\mathrm{n}}$ is used to calculate the normal contact force $F_{\mathrm{n}}$ based on the Hertz model [27],

$$F_{\mathrm{n}}=\frac{4}{3}{E}^{*}{r}^{*\frac{1}{2}}{\delta }_{\mathrm{n}}^{\frac{3}{2}}$$

(3)

where

$$\frac{1}{E^{*}}=\frac{1-{\zeta }_1^2}{E_1}+\frac{1-{\zeta }_2^2}{E_2}$$

(4)

$$\frac{1}{r^{*}}=\frac{1}{r_1}+\frac{1}{r_2}$$

(5)

In Equation (4), E₁, E₂ and ζ₁, ζ₂ are the Young’s moduli and Poisson’s ratios of the two particles in contact, respectively. In Equation (5), $r_1$ and $r_2$ represent the radii of the two cylindrical particles in contact. To consider the dissipation of kinetic energy in the particle–particle contact process, a damping force $F_{\mathrm{nd}}$ is added to the normal contact force, which has the expression,

$$F_{\mathrm{nd}}=2\sqrt{\frac{5}{6}}\beta \sqrt{K_{\mathrm{n}}{m}^{*}}{v}_{\mathrm{n}}^{\mathrm{r}}$$

(6)

where

$$K_{\mathrm{n}}=\frac{d{F}_{\mathrm{n}}}{d{\delta }_{\mathrm{n}}}=2E^{*}{R}^{*\frac{1}{2}}{\delta }_{\mathrm{n}}^{\frac{1}{2}}$$

(7)

$$m^{*}=\frac{m_1{m}_2}{m_1+m_2}$$

(8)

and $v_{\mathrm{n}}^{\mathrm{r}}$ is the normal component of the relative velocity at the contact point between two particles. In Equation (8), $m_1$ and $m_2$ are the masses of two contacting particles. The damping coefficient β in Equation (6) is a function of coefficient of normal restitution $e$ (i.e., the ratio of post-collisional to pre-collisional the relative normal velocity at the contact point).

$$\beta =-\frac{\mathit{ln}e}{\sqrt{{\mathsf{\pi }}^2+{\mathit{ln}}^{2}e}}$$

(9)

In this study, the cylindrical particles are frictionless without tangential contact forces. It should be noted that the same contact force models of Equations (3) and (6) are used for each contact type without considering the difference in contact geometry. This simplification works well for a number of packing and flow processes as discussed in the following section.

2.2. DEM Code Validation

To validate the present DEM code, particle packing, hopper discharge, and granular shear flow with the cylindrical particles were simulated in our previous studies. In the particle packing [28], cylindrical particles were dropped into a cylindrical container one-by-one, and the packing density was evaluated after all the particles settled with nearly zero velocities. The packing densities of AR = 6 cylinders obtained from the DEM simulations were consistent with experimental results for two different dropping heights. In the hopper discharge of cylindrical particles with different particle aspect ratios [29], similar particle discharge rates were obtained for the DEM simulations and experiments. Simple shear flows of monodispersed cylinders were also simulated using the present code [30]. It was observed that the stresses of the cylindrical particles with AR = 1 were close to the prediction of kinetic theory [1], which usually predicts the stresses of spheres. This consistency is reasonable as the difference in sphericity is small between the spheres and the cylinders with AR = 1. It was also observed that the stress discrepancy between the spheres and the cylinders increased as the sphericity of the cylinders decreased (i.e., as AR increased). This dependence of stresses on particle sphericity was qualitatively consistent with previous findings [31]. As a result, the algorithms used in the present cylindrical particle DEM code are verified as the agreement is achieved between the simulation and experimental results.

2.3. Computational Set-up of Shear Flow

Shear flows of binary mixtures of rods and disks in the absence of an interstitial fluid are simulated using the DEM code. The gravitational force is switched off in order to achieve a uniform field of stress and particle distribution in the shear domain. As shown in Figure 2, a well-mixed binary mixture of rods and disks is uniformly generated in a rectangular domain of dimensions L_x = 0.0072 m, L_y = 0.0072 m, and L_z = 0.0036 m (i.e., L_x/d_v = 30.25, L_y/d_v = 30.25, and L_z/d_v = 15.13). The properties of the rods (AR = 6) and disks (AR = 0.1, 0.15, 0.2, and 0.3) are presented in

Table 1. Particle properties and simulation parameters.

AR	dv (m)	ρ (kg/m3)	E (Pa)	ζ	$\mathit{e}$	γ (s−1)
0.1, 0.15, 0.2, 0.3, and 6	2.38 × 10−4	2500	8.7 × 109	0.3	0.9	100

. All of the particles have the same volume and thus the same equivalent volume diameter d_v. An approach to generate a uniform, densely packed bed of cylindrical particles in a specified domain was proposed in [25].

In the simulations, shear flow is initiated by applying a linear x-velocity profile of particles in the y-direction,

$$V_{\mathrm{ix}}=\gamma ·\left(y_{\mathrm{i}}-Y_{\mathrm{c}}\right)$$

(10)

where $V_{\mathrm{ix}}$ is the x-component of the velocity of particle i, $\gamma$ is the assigned shear rate, $y_{\mathrm{i}}$ and $Y_{\mathrm{c}}$ are the y-coordinate of the mass center of particle i and the y-coordinate of the center of the shear domain, respectively. Periodic boundary conditions are used in the x and z directions, and a Lees–Edwards boundary condition [32] is specified in the y-direction. According to the Lees–Edwards boundary condition, when the mass center of a particle moves out of the upper (or lower) boundary, it is reintroduced into the shear domain from the lower (or upper) boundary with the particle x-velocity minus (or plus) $\gamma L_{\mathrm{y}}$, by which a constant average shear rate is maintained in the domain. Meanwhile, the x-coordinate of the reintroduced particle is modified to account for the cumulative shear deformation [32], the y-coordinate is minus (or plus) $L_{\mathrm{y}}$, and the z-coordinate remains unchanged. Thus, a granular shear flow occurs in the x-y plane. It is noted that during the shear flow process, the particle x-velocity profile always maintains approximately linear as described by Equation (10).

In the shear flow simulations, the average stress tensor of the particle phase within the shear domain can be measured. The total stress tensor ${\mathit{\sigma }}^{\mathrm{tot}}$ consists of a kinetic component ${\mathit{\sigma }}^{\mathrm{kin}}$ and a collisional component ${\mathit{\sigma }}^{\mathrm{col}}$, i.e., ${\mathit{\sigma }}^{\mathrm{tot}}={\mathit{\sigma }}^{\mathrm{kin}}+{\mathit{\sigma }}^{\mathrm{col}}$. The kinetic component ${\mathit{\sigma }}_{\mathrm{kin}}$ arises from the particle velocity fluctuation and is written as,

$${\mathit{\sigma }}^{\mathrm{kin}}={\displaystyle \sum }_{i=1}^N\left(m_{\mathrm{i}}{\mathit{V}}_{\mathrm{i}}\otimes {\mathit{V}}_{\mathrm{i}}\right)/V_{\mathrm{domain}}$$

(11)

in which ${\mathit{V}}_{\mathrm{i}}={\mathit{v}}_{\mathrm{i}}-\langle \mathit{v}\rangle$ is the fluctuating velocity of particle i, and $\langle \mathit{v}\rangle$ represents average particle velocity in a local steady state, which can be determined using Equation (10) because the average velocity profile follows such a profile at the steady state. The parameter N is the total number of particles in the domain of total volume $V_{\mathrm{domain}}$ and ⊗ represents the dyadic product operator.

The collisional stress tensor, which arises from the particle–particle contact forces, is given by

$${\mathit{\sigma }}^{\mathrm{col}}={\displaystyle \sum }_{\mathrm{j}=1}^{N_{\mathrm{c}}}\left({\mathit{F}}_{\mathrm{cj}}\otimes {\mathit{l}}_{\mathrm{cj}}\right)/V_{\mathrm{domain}}$$

(12)

in which ${\mathit{F}}_{\mathrm{cj}}$ is the contact force, ${\mathit{l}}_{\mathrm{cj}}$ is a vector connecting the centers of mass of two contacting particles, and $N_{\mathrm{c}}$ is the number of contacts. A steady state of shear flow occurs when a balance is reached between shear production (from the upper and lower boundaries) and energy dissipation (via particle–particle collisions). At the steady state, the physical quantities, such as stresses and number of contacts, fluctuate slightly around constant values.

2.4. Sensitivity to Time Step and Domain Size

In previous DEM simulations with spherical particles, a critical time step was determined based on the Rayleigh wave propagation on a particle [33]. In this work, a critical time step for DEM modeling of cylindrical particles is proposed by a simple modification to the expression of critical time step for spherical particles,

$$\Delta t_{\mathrm{c}}=\frac{\pi l_{\mathrm{c}}}{0.3262+1.7532}\sqrt{\frac{\rho }{G}}$$

(13)

where the characteristic particle size $l_{\mathrm{c}}$ is the minimum value between the cylindrical particle diameter $d$ and length $l$, and $G$ is the shear modulus ($G=0.5\cdot E/\left(1+\zeta \right)$). The real time step $\Delta t$ used in the simulations is usually a small fraction of the critical time step, i.e., $\Delta t=frac·\Delta t_{\mathrm{c}}$, in which $frac$ is between 0 and 1. Different time steps of $frac$ = 0.025, 0.05, 0.075, 0.1, 0.3, and 0.6 are used to check the sensitivity in the shear flow simulations of the binary mixture of AR = 0.3, AR = 6 and C = 0.5 at the solid volume fraction $\nu$ = 0.4. The rods fraction C is defined as the number fraction of rods, i.e., $C=N_{\mathrm{rod}}/\left(N_{\mathrm{rod}}+N_{\mathrm{disk}}\right)$, in which $N_{\mathrm{rod}}$ and $N_{\mathrm{disk}}$ are the numbers of rods and disks, respectively, in a binary mixture. As shown in Figure 3a, the total shear stresses ${\sigma }_{\mathrm{xy}}^{\mathrm{tot}}$, normalized by $\rho d_{\mathrm{v}}^2{\gamma }^2$, are nearly unaffected by the four different time steps, indicating that numerical convergence has been achieved. Thus, a relatively larger time step with $frac$ = 0.1 is chosen for the rest of simulations in order to reduce the computational cost.

The base shear flow domain has the dimensions [L_x = 30.25 d_v, L_y = 30.25 d_v, L_z = 15.13 d_v]. To examine the effect of domain size, two larger domains are considered: [1.5L_x, 1.5L_y, 1.5L_z] and [2L_x, 2L_y, 2L_z]. A comparison of simulation results with the three different domains is made in Figure 3b. It can be seen that a similar evolution of normalized shear stress is obtained and the average stresses at steady state are nearly the same with the three different domains. Therefore, the base flow domain is sufficiently large to produce domain size-independent results. The stress results of a dilute flow at $\nu$ = 0.05 are shown in Figure 3b. Nevertheless, the domain of [L_x = 30.25 d_v, L_y = 30.25 d_v, L_z = 15.13 d_v] can also produce the domain size-independent results for the dense flows, due to the fact that a larger number of particles (in a dense flow than in a dilute flow) leads to a better statistical evaluation. As a result, the base domain of dimensions [L_x = 30.25 d_v, L_y = 30.25 d_v, L_z = 15.13 d_v] has been used for all the shear flow simulations reported in the following sections. In the shear flow process, the volume of the shear domain is kept constant.

3. Particle Stacking and Alignment

In dense shear flows, rods and disks exhibit significant stacking and ordering behavior. By introducing parameters to quantify the stacking and ordering, the effect of microstructure on bulk stresses is explored in this section.

Stacking occurs to the particles of the same shape in binary flows. A stacking structure can be defined as an assembly of rods or disks in which two neighboring particles are close to each other and have similar orientation. As shown in Figure 4a, a mathematical criterion is proposed to determine whether two neighboring rods or disks are within a stacking structure,

$${\left(x_{\mathrm{i}}-x_{\mathrm{j}}\right)}^2+{\left(y_{\mathrm{i}}-y_{\mathrm{j}}\right)}^2+{\left(z_{\mathrm{i}}-z_{\mathrm{j}}\right)}^2\le D^2$$

(14)

$$\arccos \left|m_1\cdot n_1+m_2\cdot n_2+m_3\cdot n_3\right|\le \varphi$$

(15)

where (x_i, y_i, z_i) and (x_j, y_j, z_j) are coordinates of the mass centers of two neighboring particles, and (m₁, m₂, m₃) and (n₁, n₂, n₃) are the unit vectors of the major axes of the two particles. According to Equations (14) and (15), two rods or disks are within a stacking when the distance between their centers is smaller than or equal to D and the angle between their major axes is smaller than or equal to $\varphi$. In the present study, the threshold values are chosen as D = 1.25d for rods or 1.5l for disks and $\varphi$ = 10° for both rods and disks. Based on the above criterion, stacking structures are determined and shown in Figure 4b for a shear flow of a binary mixture of rods (AR = 6) and disks (AR = 0.3) at the solid volume fraction of $\nu$ = 0.4. The particles that are not in the stackings are omitted in the image of Figure 4b.

To quantify and compare the degree of particle stacking, average stack size R_ave and fraction of stacking particles P are introduced. The average stack size R_ave is written as,

$$R_{\mathrm{ave}}=\frac{{{\displaystyle \sum }}_{\mathrm{i}=1}^{N_{\mathrm{st}}^{\mathrm{tot}}}{R}_{\mathrm{i}}}{N_{\mathrm{st}}^{\mathrm{tot}}}$$

(16)

in which R_i is the number of particles in the stack i and $N_{\mathrm{st}}^{\mathrm{tot}}$ is the total number of stacks. The fraction of stacking particles P is defined as,

$$P=\frac{N_p^{st}}{N_p^{tot}}$$

(17)

where $N_p^{\mathrm{st}}$ is the number of particles that are included in the stacks and $N_p^{\mathrm{tot}}$ is the total number of particles in the flow domain.

In dense flows of rods or disks, particles exhibit strong orientational preference with all the particle major axes aligned in a similar direction. An order parameter, S, was used by Börzsönyi [34] and Wegner [35] to quantify the degree of this orientational preference or particle alignment. The components of a second-order particle orientation tensor are written as,

$$Q_{ij}=\frac{3}{2N}{\displaystyle \sum }_{k=1}^N\left({\lambda }_i^k{\lambda }_j^k-\frac{1}{3}{\delta }_{ij}\right)$$

(18)

in which ${\lambda }_i^k$ and ${\lambda }_j^k$ are the i and j components of the unit vector of the major axis of particle $k$, ${\delta }_{ij}$ is the Kronecker delta, and $N$ is the total number of particles of interest. The order parameter S is obtained as the largest eigenvalue of the 3 × 3 matrix represented by the tensor $Q_{ij}$. The rods or disks in the domain are randomly oriented with the order parameter S = 0 and completely aligned in the same direction with S = 1. The parameter S usually runs between 0 and 1, and a larger value of S indicates a larger degree of particle alignment.

3.1. Monodisperse flows

Particle stacking and ordering are first discussed for monodisperse flows. Shear flows of monodisperse disks and rods at the solid volume fraction of ν = 0.5 are simulated. At steady state, stacking structures of disks and rods are generated as shown in Figure 5. The disks tend to stack face-to-face in the direction perpendicular to the flow direction, and the rods tend to stack band-by-band into bundles oriented in the flow direction.

The evolution of normalized shear stress and average stack size $R_{\mathrm{ave}}$ with cumulative shear strain for monodisperse shear flows is shown in Figure 6a,b. For disks with AR = 0.3, both the stress and average stack size fluctuate slightly around constant values in the flow process. For rods with AR = 6, the stress initially decreases and then remains nearly unchanged after a cumulative shear strain of 20, and the average stack size $R_{\mathrm{ave}}$ initially increases and then fluctuates slightly around a constant value after a strain of 20. Thus, it shows that the stress decreases as the average stack size increases. Similar trends are observed for the fraction of stacking particles P, as shown in Figure 6c,d. The stress and P remain constant for AR = 0.3 disks, and the stress decreases as P increases before they reach the steady-state values for AR = 6 rods. Consistent with the variation of shear stress with the stack size and fraction of stacking particles, the shear stress decreases as the order parameter S increases, as shown in Figure 6e,f.

It is noted that the normal stress shows a similar dependence on the stack size, fraction of stacking particles, and order parameter as the shear stress. Thus, the formation of stacking structures and increased particle alignment tends to reduce the stresses.

3.2. Binary Mixture Flows

Shear flows of binary mixtures of rods and disks, which have the same particle volume, are simulated to examine the effect of particle shape difference on the stacking and ordering behavior. Figure 7 shows the evolution of normalized shear stress and average stack size $R_{\mathrm{ave}}$ with cumulative shear strain for binary mixture flows with various rods fractions $C$. In the binary mixtures, rods and disks have aspect ratios of AR = 6 and AR = 0.3, respectively, and the solid volume fraction is set to ν = 0.5. A stack is formed by particles of the same shape, as shown in Figure 4b. Thus, the stack sizes are measured for the rods and disks separately. As shown in Figure 7a, for a binary mixture with C = 0.25, the average shear stress initially decreases and then approaches a steady-state value. The average stack sizes, which are the same for the rods and disks, remain constant in the flow process. As the rods fraction C increases (Figure 7b,c), an initial decrease in the shear stress is also obtained; however, a larger magnitude of the stress decrease is observed. The average stack size of the disks with AR = 0.3 still remains constant, while the average stack size of the rods with AR = 6 shows a slight increase with the cumulative shear strain.

The variation of normalized shear stress and fraction of stacking particles P with cumulative shear strain is plotted in Figure 8 for the binary mixture flows with various rods fractions C. The fraction of stacking particles is calculated for each particle species

$$P_{\mathrm{i}}=\frac{N_{\mathrm{pi}}^{\mathrm{st}}}{N_{\mathrm{pi}}^{\mathrm{tot}}}$$

(19)

where $N_{pi}^{\mathrm{st}}$ is the number of particles of species i that are included in the stacks and $N_{pi}^{\mathrm{tot}}$ is the total number of particles of species i in the mixture. It can be seen from Figure 8 that as the rods fraction $C$ increases, the fraction of stacking rods with AR = 6 increases and the fraction of stacking disks with AR = 0.3 decreases. Thus, the concentration of a particle species can promote stacking of that species.

For a binary mixture, Equation (17) can be used to calculate the overall fraction of stacking particles $P_{\mathrm{tot}}=N_{\mathrm{p}}^{\mathrm{st}}/N_{\mathrm{p}}^{\mathrm{tot}}$, in which $N_{\mathrm{p}}^{\mathrm{st}}$ is the total number of particles (including rods and disks) in the stacks and $N_{\mathrm{p}}^{\mathrm{tot}}$ is the total number of all the particles in the mixture. The average normalized shear stress and the average value of $P_{\mathrm{tot}}$ at the flow’s steady state are plotted as a function of rods fraction C in Figure 9. The rods fraction C = 0 corresponds to monodisperse disk flow and C = 1 corresponds to monodisperse rod flow. Figure 9 shows that the shear stress decreases monotonically with an increase of C, while the overall fraction of stacking particles $P_{\mathrm{tot}}$ shows a U-shape dependence on C. A stack of particles is a structure with geometric tessellation and mechanical stability. The face-on-face stacking of disks and band-on-band stacking of rods can achieve such compacted tessellation with good mechanical stability. In binary mixtures, due to the interaction between the particles of different shapes (i.e., disks and rods), stacking is interrupted. As a result, the overall fractions of stacking particles $P_{\mathrm{tot}}$ for binary mixtures are smaller than those of monodisperse rods and disks. The worst stacking occurs in the well-mixed binary mixture with 50% rods and 50% disks (i.e., C = 0.5).

Evolution of the normalized shear stress and order parameter S with cumulative shear strain is plotted in Figure 10 for binary mixture flows with various rods fractions. In Figure 10, the order parameter S is calculated for each particle species in the mixture. At the early stage of the flow process, as the stress declines, the order parameter S for the rods increases. The stress and order parameter reach steady state at the same cumulative shear strain. This observation is consistent with that of the monodisperse rod flow (see Figure 6f). The order parameter S for the disks fluctuates between 0.5 and 0.6 in the flow process, which is similar to the observation of monodisperse disk flow (see Figure 6e).

In a binary mixture, given the order parameter of rods $S_{\mathrm{rod}}$ and order parameter of disks $S_{\mathrm{disk}}$, a volume-averaged order parameter S_ave can be defined to quantify the degree of particle alignment for the whole mixture

$$S_{\mathrm{ave}}=C\cdot S_{\mathrm{rod}}+\left(1-C\right)\cdot S_{\mathrm{disk}}$$

(20)

The average normalized shear stress and the average order parameter $S_{\mathrm{ave}}$ are plotted in Figure 11 as a function of rods fraction C for the binary mixtures. The stress declines monotonically with increasing C and the order parameter $S_{\mathrm{ave}}$ increases monotonically with increasing C. This indicates that the stress decreases as the order parameter $S_{\mathrm{ave}}$ increases. The correlation between the system stresses and order parameter is further discussed in the following sections.

4. Effect of Particle Shape Polydispersity on Stresses

Different binary mixtures with the same rods (AR = 6) but different disks have been modeled. Disks of AR = 0.3, 0.2, 0.15 and 0.1 (the smaller AR corresponds to the flatter disks) have been considered in this work. Average normalized shear stress at steady state is plotted as a function of solid volume fraction $\nu$ for various binary mixtures in Figure 12. In general, the stress curves of binary mixtures exhibit a W-shape, similar to those of monodisperse cylindrical particle flows [30]. For binary mixtures of AR = 0.3 and AR = 6 (Figure 12a), the stresses of the mixtures (with $0<C<1$) are bounded by those of the two monodisperse flows (with $C=0$ and $C=1$). When the disks become flatter by reducing the AR, as shown in Figure 12b–d, the stresses of the mixtures are sandwiched between those of the two monodisperse systems at low solid volume fractions ($\nu$ < 0.4). Nevertheless, it is interesting to observe that the stresses of the mixtures can be larger than the stresses of both corresponding monodisperse systems at large solid volume fractions ($\nu$ > 0.4).

For a clearer comparison between the stresses of binary and monodisperse systems, the average normalized shear stresses are plotted as a function of rods fraction C in Figure 13. The stresses show a monotonic change with rods fraction at each solid volume fraction considered for the mixtures of AR = 0.3 and AR = 6 (Figure 13a). For the mixtures of AR = 0.2 and AR = 6 (Figure 13b), the stress decreases monotonically with increasing C when $\nu$ ≤ 0.4 and the stress curve exhibits a hump shape when $\nu$ = 0.5, demonstrating that the stresses of the mixtures can be greater than those of the monodisperse systems. For the mixtures of AR = 0.15 and AR = 6 (Figure 13c), the stress is independent of rods fraction C at solid volume fractions of $\nu$ ≤ 0.4, due to the fact that two monodisperse flows (C = 0 and 1) have the same stress. The humped stress curves are observed when the solid volume fraction is $\nu$ > 0.4, and the stress difference between the mixtures and monodisperse systems increases as $\nu$ increases. For the mixtures of AR = 0.1 and AR = 6 (Figure 13d), the stress increases monotonically with increasing C when $\nu$ ≤ 0.4, because the stress of monodisperse rods (C = 1) is larger than that of the monodisperse disks (C = 0). The humped stress curve occurs when $\nu$ = 0.5.

As discussed previously, the results of Figure 13 show that increased particle shape difference and increased solid volume fraction lead to larger stresses for the binary mixtures when compared to the monodisperse systems. It is presumed that the large particle shape difference makes the particles more difficultly tessellate in a confined space for the dense flows, thus promoting particle–particle interaction. As a result, larger stresses are obtained for the mixtures than for the monodisperse systems. Evidence for the enhanced particle–particle interaction is provided by considering the particle rotational kinetic energies. As shown in Figure 14a, the average rotational kinetic energy of particles generally decreases with increasing rods fraction $C$ for the binary mixtures of AR = 0.3 and 6 at the solid volume fraction ν = 0.5. This trend is consistent with the variation of the stress with $C$, as shown in Figure 13a. For the binary mixtures of AR = 0.15 and 6 at ν = 0.55, larger average rotational kinetic energies of particles are obtained for the mixtures (Figure 14b), and consistently larger stresses are obtained for these mixtures (Figure 13c). Thus, larger stresses of the mixtures are caused by the enhanced particle–particle interaction, resulting in larger particle rotational kinetic energies. The present simulation results also show that the average translational particle kinetic energy is independent of the rods fraction $C$ (not shown here) and therefore, the translational motion of particles is mainly determined by the assigned mean flow field. In addition, for the binary mixtures of AR = 0.15 and 6 at ν = 0.55, smaller average order parameters $S_{\mathrm{ave}}$ are observed compared to those of the monodisperse flows, as shown in Figure 15. The reduction in the order parameters $S_{\mathrm{ave}}$, i.e., the degree of particle alignment, is consistent with the enhanced particle rotation (Figure 14b).

5. Correlation between Stresses and Order Parameter

As discussed in the previous sections, the particle-phase stresses exhibit a strong dependence on the order parameter $S$ in dense granular flows ($\nu$ ≥ 0.4). In the shear flow simulations of monodispersed cylindrical particle flows by Berzi et al. [36], they found that the viscosity, i.e., the ratio of shear stress to shear rate, was strongly dependent on the order parameter at large solid volume fractions, for which particle alignment was significant. In this section, the correlation between the particle-phase stresses and the order parameter for binary mixtures is discussed.

Different binary mixtures of different combinations of particle aspect ratios (AR = 0.1 and 6, AR = 0.15 and 6, AR = 0.2 and 6, and AR = 0.3 and 6) with different rods fractions (C = 0.25, 0.5, and 0.75) are modeled in the present simulations. Variations of the normalized shear stress, normalized pressure, and ratio of shear stress to pressure with average order parameter $S_{\mathrm{ave}}$ is plotted in Figure 16 for various binary mixtures at a solid volume fraction of $\nu$ = 0.5. The pressure $P_{\mathrm{tot}}$ is the average value of three normal stress components, i.e., $P_{\mathrm{tot}}=\left({\sigma }_{\mathrm{xx}}^{\mathrm{tot}}+{\sigma }_{\mathrm{yy}}^{\mathrm{tot}}+{\sigma }_{\mathrm{zz}}^{\mathrm{tot}}\right)/3$. As shown in Figure 16a,b, the data for shear stresses and pressures collapse on master curves, especially for $S_{\mathrm{ave}}$ > 0.5. Some scattering data are observed at smaller order parameters. These results indicate that the particle-phase stresses have a strong dependence on the order parameter when the particle alignment becomes significant with $S_{\mathrm{ave}}$ > 0.5, regardless of particle aspect ratios and compositions of the binary mixtures; While the particle shape and the composition of a mixture also play a role in determining the stresses of a randomly oriented particle bed with a small order parameter. In dense flows, the collisional stress component ${\mathit{\sigma }}^{\mathrm{col}}$ determines the total stress [30]. The collisional stress ${\mathit{\sigma }}^{\mathrm{col}}$ depends on the particle–particle contact forces ${\mathit{F}}_{\mathrm{cj}}$ (see Equation (11)). The alignment of particles can reduce the projected area of a particle in the plane perpendicular to the flow, reducing the interference of neighboring particles and thus decreasing the particle–particle contact forces. As a result, the particle alignment-dependence of the particle–particle contact forces ${\mathit{F}}_{\mathrm{cj}}$ may lead to the order parameter-dependence of the particle phase stresses. Since the shear stresses and pressures are functions of the order parameter $S_{\mathrm{ave}}$, their ratio, which is known as the bulk friction coefficient, also depends on $S_{\mathrm{ave}}$, as shown in Figure 16c.

The normalized shear stress and pressure are assumed to follow the exponential correlations,

$$\frac{{\sigma }_{\mathrm{xy}}^{\mathrm{tot}}}{\rho d_{\mathrm{v}}^2{\gamma }^2}=\exp \left(a_1\cdot S_{\mathrm{ave}}+b_1\right)$$

(21)

$$\frac{P_{\mathrm{tot}}}{\rho d_{\mathrm{v}}^2{\gamma }^2}=\exp \left(a_2\cdot S_{\mathrm{ave}}+b_2\right)$$

(22)

The coefficients $a_1$ and $b_1$ can be determined by fitting to the data in Figure 16a, and $a_2$ and $b_2$ can be determined by fitting to the data in Figure 16b. The values of $a_1$, $b_1$, $a_2$, and $b_2$ are provided in

Table 2. Coefficients in Equations (21) and (22) for binary flows at the solid volume fraction of ν = 0.5.

Coefficients	${\mathit{a}}_1$	${\mathit{b}}_1$	${\mathit{a}}_2$	${\mathit{b}}_2$
Values	−5.0174	2.2389	−3.3785	2.8122

. Dividing Equation (21) by Equation (22) gives,

$$\frac{{\sigma }_{\mathrm{xy}}^{\mathrm{tot}}}{P_{\mathrm{tot}}}=\frac{\exp \left(a_1\cdot S_{\mathrm{ave}}+b_1\right)}{\exp \left(a_2\cdot S_{\mathrm{ave}}+b_2\right)}=\exp \left[\left(a_1-a_2\right)S_{\mathrm{ave}}+b_1-b_2\right]$$

(23)

The curve on Figure 16c represents Equation (23).

The solid volume fraction also has a significant impact on the stresses. Figure 17 shows the variation of normalized shear stress, normalized pressure, and ratio of shear stress to pressure with average order parameter at various solid volume fractions. The binary mixtures consist of the disks with AR = 0.15 and the rods with AR = 6. It can be seen that different exponential expressions with different sets of coefficients ($a_1$, $b_1$, $a_2$, and $b_2$) can be determined for the different solid volume fractions. Thus, the four coefficients, $a_1$, $b_1$, $a_2$, and $b_2$, can be expressed as functions of solid volume fraction $\nu$, and polynomial functions are assumed as follows,

$$a_1, b_1, a_2, or b_2=A_1{\nu }^3+A_2{\nu }^2+A_3\nu +A_4$$

(24)

By fitting to the obtained data, the coefficients in Equation (24) can be determined for $a_1$, $b_1$, $a_2$, and $b_2$, as presented in

Table 3. Coefficients in Equation (24).

	${\mathit{A}}_1$	${\mathit{A}}_2$	${\mathit{A}}_3$	${\mathit{A}}_4$
$a_1$	3235.0	−4730.0	2258.6	−356.60
$b_1$	2243.7	−3343.9	1617.2	−257.18
$a_2$	−2439.4	3610.5	−1731.9	270.80
$b_2$	−1643.9	2506.1	−1222.6	193.67

.

Equation (23) can be rewritten as,

$${\left[\frac{{\sigma }_{\mathrm{xy}}^{\mathrm{tot}}}{P_{\mathrm{tot}}}\exp \left(b_2-b_1\right)\right]}^{\frac{1}{a_1-a_2}}=\exp \left(S_{\mathrm{ave}}\right)$$

(25)

The term on the left-hand side of the above equation is the bulk friction coefficient modified by the parameters which can be expressed as functions of the solid volume fraction ν (Equation (24)). The modified bulk friction coefficient is plotted against the order parameter $S_{\mathrm{ave}}$ in Figure 18. The data for the various binary mixtures, which are presented in Figure 17c, tend to collapse on an exponential curve in Figure 18. The effects of particle alignment, compositions of the mixtures, and solid volume fractions $\nu$ on the bulk friction coefficient have been considered in the empirical correlation shown in Figure 18. It should be noted that the results presented in Figure 16, Figure 17 and Figure 18 include not only the steady-state data, but also the data before the steady flow is achieved, illustrating the strong order parameter-dependence of the stresses.

6. Conclusion

Granular shear flows of binary mixtures of frictionless rods and disks are simulated using the Discrete Element Method (DEM). In the binary mixtures, rods and disks have the same volume but significantly different shapes. Thus, the effects of particle shape differences on the binary mixture flows are investigated in this work. In shear flows, the rods and disks show stacking and ordering behavior, which is enhanced as the flow develops before reaching a steady state. By changing the rods fraction in the mixture, the average stack size changes slightly. The number fraction of the particles that are included in the stacks is smaller for the binary mixtures than for the monodisperse rods or disks, because the probability to form a stacking structure is reduced by the contacts between the particles of different shapes.

A summary of the main results and conclusions on the stresses of the binary mixtures of rods and disks are presented as follows:

♦

For the binary mixtures of the rods with AR = 6 and the disks with AR = 0.3, the stresses of the mixtures are bounded between the stresses of the monodisperse rods (AR = 6) and the monodisperse disks (AR = 0.3). At a given solid volume fraction, as the rods fraction increases from 0 to 1, the stress of the mixture changes monotonically from the stress of the monodisperse disks (AR = 0.3) to the stress of the monodisperse rods (AR = 6).

♦

For the binary mixtures of the rods with AR = 6 and the disks with AR = 0.2, 0.15 and 0.1, the stresses of the mixtures are still bounded by the stresses of corresponding monodisperse systems at the low solid volume fractions ($\nu$ < 0.4), while larger stresses are obtained for the mixtures than for the monodisperse systems at the high solid volume fractions ($\nu$ ≥ 0.4). It is found that the average particle rotational kinetic energies of these mixtures are greater than those of the monodisperse rods and disks in the dense flows, because larger difference in particle shape causes difficulties for the dense packing, leading to stronger particle–particle interaction in the flow. Due to the enhanced particle rotation, smaller order parameters, which is a measure of the degree of particle alignment, are obtained for the mixtures than the monodisperse systems.

♦

The data of the stresses versus the order parameter for dense binary flows with different particle shapes and different compositions can collapse on exponential function curves, indicating the strong dependence of the stresses on the particle alignment. Based on the simulation results, an empirical expression of bulk friction coefficient (ratio of the shear stress to pressure) for the dense binary flows is proposed to account for the effects of the particle alignment and solid volume fraction.

♦

The previous dense granular flow models, such as the $\mu \left(I\right)$-theory [4], have been successfully applied to the spherical particle flows. However, the effect of particle shape is difficult to be taken into account in these models. The present results (Equations (21)–(25)) may give a hint that the order parameter S should be included in the rheological laws for the elongated and flat particle flows.