Discrete Element Method Investigation of Binary Granular Flows with Di ﬀ erent Particle Shapes

: The e ﬀ ects of particle shape di ﬀ erences on binary mixture shear ﬂows are investigated using the Discrete Element Method (DEM). The binary mixtures consist of frictionless rods and disks, which have the same volume but signiﬁcantly di ﬀ erent shapes. In the shear ﬂows, stacking structures of rods and disks are formed. The e ﬀ ects 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 di ﬀ erences, the mixture stresses are bounded by the stresses of the two corresponding monodisperse systems. However, for binary mixtures with large particle shape di ﬀ erences, the stresses of the mixtures can be larger than the stresses of the monodisperse systems at large solid volume fractions because larger di ﬀ erences in particle shape cause geometrical interference in packing, leading to stronger particle–particle interactions in the ﬂow. 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 coe ﬃ cient (ratio of the shear stress to normal stress) for dense binary ﬂows is proposed by accounting for the e ﬀ ects of particle alignment and solid volume fraction.


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, µ, was a function of the inertial number, I = γd p / p/ρ p , in which γ is local shear rate, p is local pressure, and d p and ρ p are particle diameter and density, respectively. Thus, the constitutive relation for dense flows is described by the correlation µ = µ(I), which is well known as the µ(I)-theory. Non-local effects have been considered in subsequent work [5,6]. Recently, Pähtz 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 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.

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, 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 F 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 Energies 2020, 13, 1841 4 of 25 each contact type, which are provided in [24,25], are summarized for convenience in Appendix A. The overlap size δ n is used to calculate the normal contact force F n based on the Hertz model [27], Energies 2020, 13,   In Equation (4), E1, E2 and ζ1, ζ2 are the Young's moduli and Poisson's ratios of the two particles in contact, respectively. In Equation (5), and 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 is added to the normal contact force, which has the expression, where = = 2 * * δ (7) * = + (8) and is the normal component of the relative velocity at the contact point between two particles. In Equation (8), and are the masses of two contacting particles. The damping coefficient β in Equation (6) is a function of coefficient of normal restitution (i.e., the ratio of post-collisional to precollisional the relative normal velocity at the contact point).
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.

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 = In Equation (4), E 1 , E 2 and ζ 1 , ζ 2 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 nd is added to the normal contact force, which has the expression, F nd = 2 5 6 β K n m * v r n (6) where K n = dF n dδ n = 2E * R * 1 2 δ 1 2 n (7) and v r n 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).
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.

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], Energies 2020, 13, 1841 5 of 25 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.

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. 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].
Energies 2020, 13, x FOR PEER REVIEW 5 of 26 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.

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 Lx = 0.0072 m, Ly = 0.0072 m, and Lz = 0.0036 m (i.e., Lx/dv = 30.25, Ly/dv = 30.25, and Lz/dv = 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. All of the particles have the same volume and thus the same equivalent volume diameter dv. 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, where is the x-component of the velocity of particle i, is the assigned shear rate, and 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) , 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  In the simulations, shear flow is initiated by applying a linear x-velocity profile of particles in the y-direction, where V ix is the x-component of the velocity of particle i, γ is the assigned shear rate, y i and Y 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) γL 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 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 σ tot consists of a kinetic component σ kin and a collisional component σ col , i.e., σ tot = σ kin + σ col . The kinetic component σ kin arises from the particle velocity fluctuation and is written as, (11) in which V i = v i − v is the fluctuating velocity of particle i, and v 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 domain and ⊗ represents the dyadic product operator. The collisional stress tensor, which arises from the particle-particle contact forces, is given by (12) in which F cj is the contact force, l cj is a vector connecting the centers of mass of two contacting particles, and N 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.

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, where the characteristic particle size l c is the minimum value between the cylindrical particle diameter d and length l, and G is the shear modulus (G = 0.5 · E/(1 + ζ)). The real time step ∆t used in the simulations is usually a small fraction of the critical time step, i.e., ∆t = f rac·∆t c , in which f rac is between 0 and 1. Different time steps of f rac = 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 ν = 0.4. The rods fraction C is defined as the number fraction of rods, i.e., C = N rod /(N rod + N disk ), in which N rod and N disk are the numbers of rods and disks, respectively, in a binary mixture. As shown in Figure 3a, the total shear stresses σ tot xy , normalized by ρd 2 v γ 2 , are nearly unaffected by the four different time steps, indicating that numerical convergence has been achieved. Thus, a relatively larger time step with f rac = 0.1 is chosen for the rest of simulations in order to reduce the computational cost.

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, arccos | ⋅ + ⋅ + ⋅ | ≤ where (xi, yi, zi) and (xj, yj, zj) are coordinates of the mass centers of two neighboring particles, and (m1, m2, m3) and (n1, n2, n3) 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 . In the present study, the threshold values are chosen as D = 1.25d for rods or 1.5l for disks and = 10° for both rods and disks. Based on the above criterion, stacking structures are determined and shown in Figure 4b

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, arccos 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 1 , m 2 , m 3 ) and (n 1 , n 2 , n 3 ) 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 φ.
In the present study, the threshold values are chosen as D = 1.25d for rods or 1.5l for disks and φ = 10 • for both rods and disks. Based on the above criterion, stacking structures are determined and shown in Figure   To quantify and compare the degree of particle stacking, average stack size Rave and fraction of stacking particles P are introduced. The average stack size Rave is written as, = ∑ (16) in which Ri is the number of particles in the stack i and st tot is the total number of stacks. The fraction of stacking particles P is defined as, = where st is the number of particles that are included in the stacks and 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, in which and are the i and j components of the unit vector of the major axis of particle , is the Kronecker delta, and 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 . 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.

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 faceto-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. 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, (16) in which R i is the number of particles in the stack i and N tot st is the total number of stacks. The fraction of stacking particles P is defined as, where N st p is the number of particles that are included in the stacks and N tot p 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, in which λ k i and λ k j are the i and j components of the unit vector of the major axis of particle k, δ 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.

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. 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. The evolution of normalized shear stress and average stack size R 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 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.

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 with cumulative shear strain for binary mixture flows with various rods fractions . 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 and 7c), 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. 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.

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 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 (19) where N st pi is the number of particles of species i that are included in the stacks and N tot pi 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.
Energies 2020, 13, 1841 11 of 25 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 and 7c), 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 = st tot (19) where st is the number of particles of species i that are included in the stacks and 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 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. 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 = st tot (19) where st is the number of particles of species i that are included in the stacks and 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 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 tot = N st p /N tot p , in which N st p is the total number of particles (including rods and disks) in the stacks and N tot p is the total number of all the particles in the mixture. The average normalized shear stress and the average value of P 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 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 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). For a binary mixture, Equation (17) can be used to calculate the overall fraction of stacking particles = st tot , in which st is the total number of particles (including rods and disks) in the stacks and tot is the total number of all the particles in the mixture. The average normalized shear stress and the average value of 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 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 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). 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 rod and order parameter of disks disk , a volume-averaged order parameter Save can be defined to quantify the degree of particle alignment for the whole mixture The average normalized shear stress and the average order parameter 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 increases monotonically with increasing C. This indicates that the stress decreases as the order parameter increases. The correlation between the system stresses and order parameter is further discussed in the following sections. In a binary mixture, given the order parameter of rods S rod and order parameter of disks S disk , a volume-averaged order parameter S ave can be defined to quantify the degree of particle alignment for the whole mixture The average normalized shear stress and the average order parameter S 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 ave increases monotonically with increasing C. This indicates that the stress decreases as the order parameter S ave increases. The correlation between the system stresses and order parameter is further discussed in the following sections. In a binary mixture, given the order parameter of rods rod and order parameter of disks disk , a volume-averaged order parameter Save can be defined to quantify the degree of particle alignment for the whole mixture The average normalized shear stress and the average order parameter 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 increases monotonically with increasing C. This indicates that the stress decreases as the order parameter increases. The correlation between the system stresses and order parameter is further discussed in the following sections.

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 ν 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 (ν < 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 (ν > 0.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 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 < < 1) are bounded by those of the two monodisperse flows (with = 0 and = 1). When the disks become flatter by reducing the AR, as shown in Figures 12b-d,  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 ≤ 0.4 and the stress curve exhibits a hump shape when = 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 ≤ 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 > 0.4, and the stress difference between the mixtures and monodisperse systems increases as increases. For the mixtures of AR = 0.1 and AR = 6 (Figure 13d), the stress increases monotonically with increasing C when ≤ 0.4, because the stress of monodisperse rods (C 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 ν ≤ 0.4 and the stress curve exhibits a hump shape when ν = 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 ν ≤ 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 ν > 0.4, and the stress difference between the mixtures and monodisperse systems increases  (Figure 13d), the stress increases monotonically with increasing C when ν ≤ 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 ν = 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 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 , 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 (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 are observed compared to those of the monodisperse flows, as shown in Figure 15. The reduction in the order parameters , i.e., the degree of particle alignment, is consistent with the enhanced particle rotation (Figure 14b). 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 ave are observed compared to those of the monodisperse flows, as shown in Figure 15. The reduction in the order parameters S ave , i.e., the degree of particle alignment, is consistent with the enhanced particle rotation (Figure 14b).

Correlation between Stresses and Order Parameter
As discussed in the previous sections, the particle-phase stresses exhibit a strong dependence on the order parameter in dense granular flows ( ≥ 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 is plotted in Figure  16 for various binary mixtures at a solid volume fraction of = 0.5. The pressure is the average value of three normal stress components, i.e., = + + 3 ⁄ . As shown in Figures 16a  and 16b, the data for shear stresses and pressures collapse on master curves, especially for > 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 > 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

Correlation between Stresses and Order Parameter
As discussed in the previous sections, the particle-phase stresses exhibit a strong dependence on the order parameter in dense granular flows ( ≥ 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 is plotted in Figure  16 for various binary mixtures at a solid volume fraction of = 0.5. The pressure is the average value of three normal stress components, i.e., = + + 3 ⁄ . As shown in Figures 16a  and 16b, the data for shear stresses and pressures collapse on master curves, especially for > 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 > 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 Figure 15. The average normalized shear stress and the average order parameter S ave as a function of rods fraction C for the binary mixtures of AR = 0.15 and AR = 6. The solid volume fraction is ν = 0.55.

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 (ν ≥ 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 ave is plotted in Figure 16 for various binary mixtures at a solid volume fraction of ν = 0.5. The pressure P tot is the average value of three normal stress components, i.e., P tot = σ tot xx + σ tot yy + σ tot zz /3. As shown in Figure 16a,b, the data for shear stresses and pressures collapse on master curves, especially for S 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 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 σ col determines the total stress [30]. The collisional stress σ col depends on the particle-particle contact forces F 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 F 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 ave , their ratio, which is known as the bulk friction coefficient, also depends on S ave , as shown in Figure 16c.
The normalized shear stress and pressure are assumed to follow the exponential correlations, 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. Dividing Equation (21) by Equation (22) gives, Energies 2020, 13, x FOR PEER REVIEW 17 of 26 determining the stresses of a randomly oriented particle bed with a small order parameter. In dense flows, the collisional stress component determines the total stress [30]. The collisional stress depends on the particle-particle contact forces (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 may lead to the order parameter-dependence of the particle phase stresses. Since the shear stresses and pressures are functions of the order parameter , their ratio, which is known as the bulk friction coefficient, also depends on , as shown in Figure 16c. The normalized shear stress and pressure are assumed to follow the exponential correlations, = exp ( ⋅ + ) The coefficients and can be determined by fitting to the data in Figure 16a, and and can be determined by fitting to the data in Figure 16b. The values of , , , and are provided in Table 2. Dividing Equation (21) by Equation (22) gives, The curve on Figure 16c represents Equation (23).   (21) and (22)  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 ( , , , and ) can be determined for the different solid volume fractions . Thus, the four coefficients,  ,  , , and , can be expressed as functions of solid volume fraction , and polynomial functions are assumed as follows, By fitting to the obtained data, the coefficients in Equation (24) can be determined for , , , and , as presented in Table 3.   Figure 16. Variation of (a) normalized shear stress, (b) normalized pressure, and (c) ratio of shear stress to pressure with average order parameter S ave for various binary mixtures. The solid volume fraction is ν = 0.5. (21) and (22) for binary flows at the solid volume fraction of ν = 0.5.

Coefficients
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 ν, and polynomial functions are assumed as follows, Energies 2020, 13, x FOR PEER REVIEW 18 of 26 (c) Figure 16. Variation of (a) normalized shear stress, (b) normalized pressure, and (c) ratio of shear stress to pressure with average order parameter for various binary mixtures. The solid volume fraction is = 0.5. (21) and (22)  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 ( , , , and ) can be determined for the different solid volume fractions. Thus, the four coefficients, , , , and , can be expressed as functions of solid volume fraction , and polynomial functions are assumed as follows,

Table 2. Coefficients in Equations
By fitting to the obtained data, the coefficients in Equation (24) can be determined for , , , and , as presented in Table 3.
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 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  Figure 17. Variation of (a) normalized shear stress, (b) normalized pressure, and (c) ratio of shear stress to pressure with average order parameter S ave for the binary mixtures (AR = 0.15 and AR = 6) at various solid volume fractions ν. It is noted that the vertical coordinate is the logarithmic coordinate. Thus, the stresses have the exponential dependence on S ave for various ν.
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.
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 ave in Figure 18. The data Energies 2020, 13, 1841 20 of 25 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 ν 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 Figures 16-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.
Energies 2020, 13, x FOR PEER REVIEW 20 of 26 solid volume fractions 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 Figures 16-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.

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 ( < 0.4), while larger stresses are obtained for the mixtures than for the monodisperse systems at the high solid volume fractions ( ≥ 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.  Figure 17a is also applied in this figure.

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 (ν < 0.4), while larger stresses are obtained for the mixtures than for the monodisperse systems at the high solid volume fractions (ν ≥ 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 µ(I)-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.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Determination of overlap size, contact normal vector, and contact point position for the basic contact types between two cylindrical particles. The equations, which can be found in [24], are defined in particle i's body-fixed frame of reference. The contact detection algorithms for several special contact scenarios are referred to [25].

Face-face contact
Energies 2020, 13, x FOR PEER REVIEW 21 of 26  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 ( )-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.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Table A1. Determination of overlap size, contact normal vector, and contact point position for the basic contact types between two cylindrical particles. The equations, which can be found in [24], are defined in particle i's body-fixed frame of reference. The contact detection algorithms for several special contact scenarios are referred to [25].

Contact types Equations
Face-face contact Overlap size: where Li and Lj are the major axial lengths of contacting cylinders i and j, z cyl j i is the zcoordinate of the center of particle j in particle i's body-fixed frame of reference. where L i and L j are the major axial lengths of contacting cylinders i and j, z i cyl j is the z-coordinate of the center of particle j in particle i's body-fixed frame of reference.
Contact normal vector: n i contact = [0, 0, sign (z i cyl j )] in which the function 'sign' returns the sign of the expression. where R i and R j are the radii, respectively, of the particles i and j, and x i cyl j and y i cyl j are the x and y-coordinates, respectively, of the center of the particle j. where and are the radii, respectively, of the particles i and j, and x cyl j i and y cyl j i are the x and y-coordinates, respectively, of the center of the particle j.

Face-band contact
Overlap size: Contact normal vector: in which , is the unit vector along z-axis of the particle i in the particle i's body-fixed frame of reference.
Contact point position:

Face-edge contact
Overlap size: in which z is the z-coordinate of the point E in the particle i's body-fixed frame of reference.
Contact normal vector: where , is the unit vector along z-axis in the particle i's frame of reference. In which x E i is the position of the point E in the particle i's body-fixed frame of reference.

Parallel band-band contact
Overlap size: Overlap size: δ n = L i/ 2 + R j − z i cyl j .
Contact normal vector: n i contact = sign (z i cyl j )e i cyl i,z in which e i cyl i,z is the unit vector along z-axis of the particle i in the particle i's body-fixed frame of reference.

Contact point position:
x i contact = 1 2 x i A +x i B − sign z i cyl j ( 1 2 δ n )(0,0,1) i where x i A and x i B are the position vectors of the point A and B, respectively, in the particle i's body-fixed frame of reference.

Face-edge contact
where and are the radii, respectively, of the particles i and j, and x cyl j i and y cyl j i are the x and y-coordinates, respectively, of the center of the particle j.

Contact normal vector:
in which , is the unit vector along z-axis of the particle i in the particle i's body-fixed frame of reference. In which x E i is the position of the point E in the particle i's body-fixed frame of reference.

Parallel band-band contact
Overlap size: Overlap size: δ n = L i/ 2 − z i E in which z i E is the z-coordinate of the point E in the particle i's body-fixed frame of reference.

Contact normal vector:
n i contact = sign (z i cyl j )ê i cyl i,z whereê i cyl i,z is the unit vector along z-axis in the particle i's frame of reference.

Contact point position:
x i contact = x i E + sign(z i cyl j )( 1 2 δ n )(0,0,1) i In which x i E is the position of the point E in the particle i's body-fixed frame of reference.  Table A1. Cont. Overlap size: δ n = (R i + R j )-d where d is the shortest distance between the two axes of contacting cylindrical particles.

Contact Types Equations
Contact normal vector: where x i P * and x i Q * are the positions of the intersection points of the shortest distance line with the major axes of the particle i and the particle j, respectively.
in which x i E and y i E are the x and y coordinates of the point E, which is on the edge of the particle j and has the shortest distance to the major axis of the particle i. The overlap size δ n is equal to the shortest distance between AB and CD, δ n = r i ij + µ * ê i i − λ * ê i j in which, r i ij = r i i − r i j , and r i i is the position of the center of CD, and r i j is the position of the center of AB. x i contact = 1 2 P i j + P i i