Effects of volume fraction and particle shape on the rheological properties of oblate spheroid suspensions

The rheology of non-Brownian suspensions of spherical particles has been the topic of extensive studies and reviewed elsewhere.^1–3 The bulk rheological properties of these suspensions are the manifestation of physical processes occurring at the particle scale, such as the fluid–particle interactions and inter-particle interactions. One of the key parameters is the volume fraction of particles, $ϕ$⁠, which has a direct impact on the viscosity of the suspension.^4,5 With a few exceptions,^6–10 most studies in the literature have been focused on spherical particle suspensions. In the present study, we focus instead on the ellipsoidal particle suspensions and investigate the effects of $ϕ$ as well as the particle shape on the rheological properties of these suspensions.

In a simple shear flow, the relative viscosity η_r of the suspension is defined as

where $Σ12$ is a component of the suspension stress tensor $Σ$⁠, η_f is the viscosity of the suspending fluid, and $γ̇$ is the imposed bulk shearing rate. Here, the 1, 2, and 3 directions correspond to the flow, velocity gradient, and vorticity directions of the shear flow, respectively. The first and second normal stress difference coefficients are defined as

respectively, which measure the loss of isotropy in the normal stresses. The rheology of the suspension can largely be described by η_r, $N1*$⁠, and $N2*$⁠, which are the focus of this paper.

The correlation between η_r and $ϕ$ has been well known for sphere suspensions. In the very dilute regime, η_r varies linearly with $ϕ$ according to Einstein.⁴ Other viscosity laws were proposed (see Table I) and have enjoyed considerable success in matching experimental data (see Fig. 4 of Guazzelli and Pouliquen³) for spherical particle suspensions. For the normal stress differences, studies on sphere suspensions suggest that $N2*$ is generally negatively signed, whereas $N1*$ exhibits more uncertainty and can be either positive or negative.

The motion of an ellipsoidal particle in a shearing flow was studied by Jeffery.¹³ Utilizing Jeffrey's solution, Mueller et al.⁸ combined analytical and experimental work to investigate the viscosity of ellipsoidal (both oblate and prolate) particle suspensions. They reported empirical values for the maximum packing fraction, $ϕm$⁠, and the Einstein coefficient, B. Numerical simulations were performed to study the viscosity of plate-like particle¹⁴ and oblate particle suspensions.⁹ The latter work reported values of η_r, $N1*$⁠, and $N2*$ of oblate spheroidal particle suspensions for $ϕ$ up to 0.25. The algorithm used relies on Jeffrey's solution to prescribe the particle motions and introduces equivalent spherical surfaces to approximate the spheroidal surfaces when calculating inter-particle interactions.

In this paper, we employ the lattice Boltzmann (LB) and discrete element (DE) methods to simulate the suspensions of oblate spheroidal particles. The fluid phase, whose motion is governed by Navier–Stokes equations, is solved by the LB method, while the solid phase is solved directly by the DE method. The volume fraction is varied (for $ϕ$ up to 0.40), as well as the aspect ratio (α) of the oblate particles. The aim of the study is to investigate the effects of $ϕ$ and α on the bulk rheological properties (such as η_r, $N1*$⁠, and $N2*$⁠) and quantify parameters (such as B) for oblate particle suspensions.

The coupled LB and DE methods are employed to simulate the fluid-particle flow. The algorithm is implemented through an open-source software package known as OpenLB,^15–17 which is capable of simulating arbitrary shaped particles. Two-way coupling is implemented for the momentum exchange between the particles and fluid and is performed by the procedure developed by Krause et al.,¹⁵ enforcing the “bounce-back” condition^18,19 at the particle-fluid boundary nodes. An efficient scheme for ellipsoidal particle contact detection proposed by Lin and Ng^20,21 is employed. The normal contact forces include Hertz's elastic force^22,23 and a small damping force specifically designed for ellipsoidal particles.²⁴ The tangential contact force includes the Coulomb's friction force. The physical parameters pertinent to the particle contact model are summarized in Table II. The particle Reynolds number, $Re≡γ̇(r*)2/νf$⁠, is 0.014, which lies in the Stokes flow regime. The particle-to-fluid density ratio, $ρp/ρf$⁠, is held constant at 1.2. Preliminary work done for $ρp/ρf=2.7$ (not shown) suggests no significant variation of the rheological properties with the density ratio at this range, which is consistent with the observations by Thorimbert et al.²⁵ For the LB calculations, the D3Q15 lattice with the Bhatnagar–Gross–Krook operator (see, e.g., Krüger²⁶) is employed. Following two previous LB-based studies,^15,25 no explicit lubrication force modeling is included in the LB–DE solver. A detailed documentation of the LB and DE methods for this specific configuration can be found in Chap. 3 of Guo.²⁷ 

As shown in Fig. 1(a), simulations are performed for a $150×50×50 μm$³ fluid channel. A plane Couette flow is set up by the top and bottom boundaries, and the flow is periodic in the x and z directions. The top (⁠$y=2h=50 μm$⁠) and bottom (y = 0) walls move in opposite directions with a specific speed to create the imposed shearing rate $γ̇$ (see Table II). The initial velocity field follows a linear profile, i.e., $uf(t=0)=(γ̇(y−h),0,0)$⁠. The initial position and orientation of the particles, which are randomized as much as possible, are generated by a program known as Packing Ellipsoid 3.²⁸ Four sets of simulations are performed, whose input parameters are tabulated in Table III. No effects of gravity are included in the simulations. As shown in Fig. 1(b), the diameter of the oblate particle along the symmetric axis is $2r2$⁠, and the equal diameter along the other axes is $2r1$⁠. The particle aspect ratio, $α≡r2/r1$⁠, is varied from 0.4 to 1.0 (the latter corresponds to a sphere). Ten simulations are run for each α, with $ϕ$ varied from 0.00 to 0.40. Each simulation is run up to a time corresponding to $γ̇t=20$⁠. An initial transient is typically observed in the statistics for $0<γ̇t<5$⁠. Hence, we regard $5<γ̇t<20$ as the quasi-steady state and produce time-averaged statistics for this specific time window. During the quasi-steady state, an approximately linear velocity profile (as in the initial condition) is observed for the mean flow between the two walls. Upon conducting mesh convergence tests, we set the grid spacing at $0.5 μm$⁠, time step at 41.7 ns, and the dimensionless relaxation time at 1.0 for LB calculations.²⁶ The OpenLB program is executed on the Cedar cluster operated by Compute Canada. Each simulation takes approximately 10 days to run using 40 to 100 parallel processors, depending on the number of particles in the system. Detailed description of how the rheological parameters [Eqs. (1) and (2)] are calculated from the results of the LB and DE calculations can be found in Chap. 2 of Guo.²⁷ 

The relative viscosity η_r of the spherical particle suspension [Eq. (1)] is plotted against $ϕ$ in Fig. 2. The results from our α = 1 (sphere) simulation are compared with the viscosity laws (Table I). Our numerical results reproduce the expected trend between η_r and $ϕ$⁠, and are in general agreement with the Krieger–Dougherty (KD) formula.¹² In Fig. 3(a), the η_r vs $ϕ$ relation is shown all four cases of α. For all cases, η_r increases with $ϕ$⁠. In Fig. 3(b), the ratio $ηr/ηr, α=1$ is shown, where $ηr, α=1$ is the relative viscosity of the spherical particle suspension at each $ϕ$⁠. $ηr/ηr, α=1$ is greater than unity for all $α<1$⁠, suggesting that the oblate spheroid suspensions yield a higher viscosity than their spherical counterparts. The ratio also generally increases with $ϕ$⁠. In Fig. 3(c), we focus on $ηrp$⁠, i.e., the component of η_r due to inter-particle interactions, which can be calculated from the contact stresslet outlined by Eq. (2.24) of Gallier et al.²³ As expected, $ηrp$ makes up a larger portion of η_r as $ϕ$ increases up to 0.40, as inter-particle collisions become more likely in a denser suspension.

Increased particle collisions may also result in anisotropy in the normal stresses. The normal stress difference coefficients $N1*$ and $N2*$ [Eq. (2)] are plotted against $ϕ$ in Fig. 4. Similar to the spherical particle suspension data compiled by Guazzelli and Pouliquen,³ the sign of $N1*$ is inconclusive. As the particles become more oblate, the sign of $N1*$ switches from positive to negative for the denser suspensions (⁠$ϕ≥0.20$⁠). $N2*$ appears to be negative for all cases, which is consistent with the sphere cases.³ As compared to $N1*, N2*$ seems less sensitive to α at a given value of $ϕ$⁠. Particle collisions are expected to occur most frequently within the plane of shear,³ i.e., the x–y plane, which leads to a negative $N2*$—this seems to be the case regardless of the particle shape.

To visualize any reorientation of the particles subject to shearing, we perform equal area projection (EAP) (Sec. 3.3 of Fisher et al.²⁹) to display the normal vector $n$ of each particle on the x–z plane (Fig. 5). EAP, which was previously used to visualize particle orientations by Lin and Ng²⁰ (see their Fig. 8), has the advantage of preserving the number density of the projected points over the flat surface. Two typical cases are shown in Fig. 5, both of $ϕ=0.40$⁠. For both cases, the initial distribution of particle orientation is relatively more random with some preference toward 90° and 270° (⁠$T3$⁠) angles, which is a result of the initialization scheme²⁸ at this large value of $ϕ$⁠. As the shear strain (⁠$γ̇t$⁠) increases, the particles in the $α=0.4$ case notably migrate toward the center of the circle, which corresponds to a particle with normal vector aligned with the y-axis (⁠$T2$⁠). Particles are observed to reorient such that their normal vectors are more aligned with the direction of the velocity gradient, which is consistent with the findings by numerical simulations^9,14 and laboratory measurement.³⁰ For comparison, in the case with $α=0.8$⁠, no systematic reorientation of the particles is observed. Having inspected all simulations, we found no such reorientation for cases of $α=1.0$ or 0.8, where the particles are spherical or too close to a sphere. Presumably, the reorientation is a consequence of inter-particle collisions. Therefore, one may expect the reorientation to be more likely for denser suspensions. Indeed, for $α=0.6$⁠, systematic reorientation occurs for $ϕ≥0.30$⁠; for $α=0.4$⁠, it occurs for $ϕ≥0.15$⁠.

The behavior of $N1*$⁠, which is presented in Fig. 4(a), may be linked to the particle reorientation shown in Figs. 5(a)–5(c). According to the results from particle pair interactions,³ $N1*$ is expected to be negative. It was also found recently that wall confinement^31,32 could be the main contributor to the normal stress in the wall-normal direction, leading to positive $N1*$⁠. A consequence of the particle reorientation could be that the oblate particle suspension is not as effective in transmitting the normal stress in the wall-normal y direction as the spherical case because the longer axes of the particles become more parallel to the walls, i.e., closer to the $T2$ orientation shown in Fig. 5(b). As a result, the wall confinement effects are expected to be weaker for the more oblate particle suspensions, which could be the reason for the sign change in $N1*$ observed in Fig. 4(a) as α decreases from 1.0 to 0.4. In particular, $N1*$ for the $α=0.4$ case [Fig. 4(a)] becomes more negative at larger volume fraction $ϕ$⁠, which is presumably due to increased particle pair interactions³ and weakened wall confinement effects.^31,32

Some microstructural statistics for particle contacts are presented in Fig. 6. First, the dimensionless contact distance, R_d [see Fig. 1(b)], averaged for all contact pairs in each simulation, is plotted in Fig. 6(a) as a function of $ϕ$⁠. For the sphere case (α = 1), R_d is unity by definition. For all oblate cases (⁠$α<1$⁠), R_d is observed to be greater than unity, suggesting that the average contact distance has been increased from the sphere case of the same equivalent radius $r*$ (Fig. 1). Moreover, at a given $ϕ$⁠, R_d increases as the particle becomes more oblate (i.e., as α decreases). In Fig. 6(b), the contact ratio R_c, i.e., the average number of contacts per particle, can be seen to increase with $ϕ$⁠, which is as expected. In general, R_c also increases as α decreases. The latter trend was also reported by Donev et al.³³ in their study of maximum packing fraction for non-spherical particles, which was attributed to the increase in the degrees of freedom as the particle shape becomes less isotropic. The same mechanism could take place in the shear flow although the volume fraction $ϕ$ is smaller than those studied by Donev et al.³³ The increase in R_c implies a higher likelihood for momentum transfer to occur due to inter-particle collisions and thus a higher viscosity, which is in agreement with the trend observed in Fig. 3 between viscosity and the control parameters, i.e., $ϕ$ and α. In short, inter-particle collisions are found to be more frequent on a per-particle basis in the more oblate particle suspensions, which is consistent with the observed increase in viscosity with decreasing α (Fig. 3).

In order to compile the viscosity data in Fig. 3, we formulate the best fits directly to the viscosity laws. Specifically for the KD model, one could find the best fit for both $ϕm$ and B. Instead, we use the experimental results³³ to obtain the maximum packing fraction $ϕm$ for each α and fit the KD model to our data [Fig. 3(a)] to obtain the Einstein coefficients, B. Such results are tabulated in Table IV and can complement an earlier result⁸ for $α=0.13$ (see their Table III). It can be seen that $ϕm$ reaches the maximum at $α≈0.60$³³ and B decreases monotonically with increasing α.

In conclusion, we have presented the rheological properties of oblate particle suspensions for a considerable range of volume fraction and aspect ratio. Other than these two parameters, there could be numerous other factors in determining the rheology which have been studied recently^34–37 for spherical particle suspensions. Investigating these effects for non-spherical particle suspensions as well as extending the simulations to even denser suspensions could be avenues for future research.