2.1. Structure solver and 3-D rigid body rotation equations

The surface of any arbitrarily shaped particle, ${\it\Gamma}(t)$ , is discretized with an unstructured triangular mesh, forming a list of $K$ material points. These points can be tracked over time by their global Lagrangian position vector $\boldsymbol{r}^{k}(t)$ ( $k=1,K$ ):

where $\boldsymbol{r}_{0}^{k}$ is the initial location of the $k\text{th}$ material point on ${\it\Gamma}(0)$ . The particle rotates in the background fluid due to the torque exerted by flow. The rotation matrix which describes such rotation in three-dimensional space is denoted as $\unicode[STIX]{x1D64D}$ . First, we define the shape of the rigid body in terms of a coordinate reference system attached to the body at its centre of mass (this coordinate system rotates with the same angular velocity as the rigid body and is referred to as the body frame). Therefore, the position vector of the material point $k$ remains unchanged over time in this reference frame. Having the geometric description of the rigid body in the body frame, we use $\unicode[STIX]{x1D64D}(t)$ to transform the body-frame description into the inertial frame and find the instantaneous position vector of point $k$ in the inertial coordinate system as:

where $\boldsymbol{r}_{c}$ is the position vector of the particle’s centre of mass and $\boldsymbol{r}_{0}^{k}$ is the position vector of point $k$ expressed in the body space, i.e. the position vector of the same point in the global frame at the beginning of the motion. Equation (2.2) reveals that we need to compute the rotation matrix at each specific time in order to define the orientation of the body. To do so, we should solve the dynamic equation of motion for the particle while it moves in the fluid. Motion of the immersed body in the fluid results in a hydrodynamic reacting force (torque) exerted on the surface of the body. The total external torque is defined as:

where $\boldsymbol{F}_{k}(t)$ denotes the total force from the fluid acting on the $k\text{th}$ point at time  $t$ . Knowing the total hydrodynamic torque exerted on the body, we can compute the angular momentum of the body by integrating the following equation:

where $\boldsymbol{L}$ is the angular momentum of the body defined as:

where $\unicode[STIX]{x1D644}$ and ${\bf\omega}$ are the instantaneous moment of inertia matrix and the angular velocity vector, respectively. In addition, $\unicode[STIX]{x1D644}(t)$ depends on the orientation of the body, and we can find the moment of inertia tensor in terms of the constant-value moment of inertia computed in the body reference frame and the rotation matrix as:

where $\unicode[STIX]{x1D644}_{b}$ is calculated with respect to the body-frame coordinate system and remains constant over the simulation. Considering (2.4)–(2.6), it is clear that we need one more equation to solve for all the unknown variables, including $\boldsymbol{L}(t)$ , $\unicode[STIX]{x1D644}(t)$ , ${\bf\omega}(t)$ and $\unicode[STIX]{x1D64D}(t)$ . The final equation to close the system of differential equations and fully determine the dynamics of rigid body rotation is the time evolution of the rotation matrix due to the angular velocity of the body. Instead of the conventional method of solving for $\unicode[STIX]{x1D64D}(t)$ , it is preferable to use quaternions, $\boldsymbol{q}(t)$ , to represent the orientation of the body in three-dimensional space (Shoemake Reference Shoemake1985; Baraff Reference Baraff1997; Kuipers Reference Kuipers1999) to avoid the numerical drift. Since we still need the rotation matrix to calculate the moment of inertia, i.e. (2.6), $\unicode[STIX]{x1D64D}(t)$ will be computed as an auxiliary variable from $\boldsymbol{q}(t)$ as follows:

We also need the following equation to calculate the quaternion vector during the simulation:

Now we have enough equations ((2.4)–(2.8)) to close the system of equations and solve for all the unknown variables of the particle dynamics $(\boldsymbol{L},\unicode[STIX]{x1D644},{\bf\omega},\unicode[STIX]{x1D64D},\boldsymbol{q})$ .

2.2. Flow solver

The equations governing the motion of an incompressible Newtonian viscous fluid are the unsteady 3-D Navier–Stokes equations. Non-dimensionalized by a characteristic velocity $\dot{{\it\gamma}}a$ and the particle radius  $a$ , the governing equations in curvilinear coordinates read as follows:

where $J$ is the determinant of the Jacobian of the transformation, ${\it\xi}_{q}^{r}=\partial {\it\xi}^{r}/\partial x_{q}$ , $g^{rm}={\it\xi}_{q}^{r}{\it\xi}_{q}^{m}$ are the contravariant metrics of the transformation, $\mathfrak{U}^{r}=u_{m}{\it\xi}_{m}^{r}/J$ is the contravariant volume flow rate, $u_{q}$ is the non-dimensional Cartesian velocity of the fluid, $p$ is the pressure non-dimensionalized by ${\it\rho}\dot{{\it\gamma}}^{2}a^{2}$ and $a$ is one-half the largest diameter of the particle.

Equations (2.9) are solved over curvilinear grids using a hybrid staggered/non-staggered discretization with the QUICK (Quadratic Upstream Interpolation for Convective Kinematics) scheme for the convective terms and the central difference for the viscous terms. The equations are integrated in time using an efficient, second-order-accurate fractional step methodology with a Newton–Krylov solver for the momentum equations, and a Generalized Minimal Residual Method (GMRES) solver enhanced with multigrid as a preconditioner for the Poisson equation (Ge & Sotiropoulos Reference Ge and Sotiropoulos2007; Borazjani et al. Reference Borazjani, Ge and Sotiropoulos2008). The code is parallelized using Message Passing Interface (MPI) and the parallel libraries of PETSc (Portable, Extensible Toolkit for Scientific Computation) (Balay et al. Reference Balay, Abhyankar, Adams, Brown, Brune, Buschelman, Dalcin, Eijkhout, Gropp, Kaushik, Knepley, McInnes, Rupp, Smith, Zampini and Zhang2001) to efficiently utilize the supercomputing facilities available to the authors.

To handle the moving boundaries, a sharp-interface immersed-boundary method is used (Gilmanov & Sotiropoulos Reference Gilmanov and Sotiropoulos2005; Ge & Sotiropoulos Reference Ge and Sotiropoulos2007; Borazjani et al. Reference Borazjani, Ge and Sotiropoulos2008). In this method, the fluid grid is fixed and does not move with the boundary, making it more suitable for large deformations because it does not create highly skewed grids. The sharp interface is maintained by reconstructing the boundary conditions at the nodes that are exterior to, but adjacent to the immersed-boundary surface using interpolation along the local normal to the boundary, as has been shown to be second-order accurate (Gilmanov & Sotiropoulos Reference Gilmanov and Sotiropoulos2005). The background nodes at each time step are classified into fluid, wall and immersed-boundary nodes using an efficient ray-tracing algorithm described in detail in Borazjani et al. (Reference Borazjani, Ge and Sotiropoulos2008). The method has been fully validated and applied to a wide range of problems, including cardiovascular flows (Borazjani, Ge & Sotiropoulos Reference Borazjani, Ge and Sotiropoulos2010a ; Borazjani & Sotiropoulos Reference Borazjani and Sotiropoulos2010a ; Borazjani Reference Borazjani2013a ), aquatic swimming (Borazjani et al. Reference Borazjani, Sotiropoulos, Malkiel and Katz2010c ; Borazjani & Sotiropoulos Reference Borazjani and Sotiropoulos2010b ; Borazjani Reference Borazjani2013b ; Borazjani & Daghooghi Reference Borazjani and Daghooghi2013), and vortex-induced vibrations (Borazjani & Sotiropoulos Reference Borazjani and Sotiropoulos2009).

2.3. Coupling methods for partitioned FSI problems

The structure domain is coupled to the fluid domain through the dynamic boundary condition (2.4), i.e. the torque on the right-hand side of the equation is a function of fluid variables (velocity and pressure). In addition, the fluid domain is coupled to the structure domain through the no-slip boundary condition on the immersed boundaries. Using the strong-coupling partitioned approach with under-relaxation (Borazjani et al. Reference Borazjani, Ge and Sotiropoulos2008), we solve the fluid and structured domain together. In this method, the fluid and structure domains are separately advanced in time, but the simultaneous (implicit) solution of both equations is achieved through fluid/structure sub-iterations. The strong-coupling scheme is computationally expensive since each sub-iteration requires the solution of both structural and fluid equations. However, due to its stabilizing capability, implementing strong coupling enables us to choose time steps up to 20 times larger than with loose coupling.

2.4. Formulation of rheological parameters

The last step in our numerical framework is the calculation of the particle stress from the simulated motion of fluid and particles. We implemented the approach developed by Batchelor (Reference Batchelor1970) to compute the bulk stress generated by a flowing suspension in terms of volume averages. The volume-averaged stress (non-dimensionalized by ${\it\rho}\dot{{\it\gamma}}^{2}a^{2}$ ) in a statistically homogenous suspension is calculated as:

where $V_{f}$ is the volume occupied by the fluid and $\partial \mathscr{U}_{i}/\partial x_{j}$ is the volume-averaged velocity gradient. The quantity ${\it\Sigma}_{ij}^{p}$ , named the ‘particle stress’ by Batchelor (Reference Batchelor1970), represents the contribution to the bulk stress due to the presence of particles. For rigid particles, this contribution reads as follows:

where $V_{p}$ is the particle volume, $N$ is the total number of particles in the suspension, ${\it\sigma}_{ij}$ is the Newtonian stress, and $u_{i}^{\prime }$ are velocity fluctuations. As is evident in (2.11), particle stress is generated by two mechanisms: the actual stress of the particles and the momentum flux caused by velocity fluctuations due to the presence of particles (Batchelor Reference Batchelor1970; Haddadi & Morris Reference Haddadi and Morris2014). Using the divergence theorem, the volume integral appearing in the first term of (2.11) can be broken into two integrals to clarify the origin of the particle stress:

where $x_{i}$ is the displacement vector measured from the centre of mass of the particle (choice of reference point is arbitrary) and $a_{i}$ is the particle acceleration (which comes from $\boldsymbol{{\rm\nabla}}{\bf\sigma}={\it\rho}\boldsymbol{a}$ ) calculated at $\boldsymbol{r}$ (the position vector of a point inside the particle relative to its centre of mass) as $\boldsymbol{a}={\bf\alpha}\times \boldsymbol{r}+{\bf\omega}\times {\bf\omega}\times \boldsymbol{r}$ , with ${\bf\alpha}$ being the angular acceleration of the particle. In the absence of external torque on the particles, the symmetric part of the particle stress is given by:

where $\unicode[STIX]{x1D61A}_{ij}$ is the symmetric part of the first moment of surface stress called the ‘stresslet’ (Batchelor Reference Batchelor1970), $\unicode[STIX]{x1D617}_{ij}$ is the acceleration stress, and $\unicode[STIX]{x1D619}_{ij}$ is the Reynolds stress. Neglecting Brownian effects for a suspension consisting of buoyant rigid particles of arbitrary shape, these terms are as follows:

where $\unicode[STIX]{x1D64E}$ is non-dimensionalized with respect to ${\it\rho}\dot{{\it\gamma}}^{2}a^{5}$ , $\unicode[STIX]{x1D64B}$ and $\unicode[STIX]{x1D64D}$ are non-dimensionalized with respect to ${\it\rho}\dot{{\it\gamma}}^{2}a^{2}$ , and all variables on the right-hand side of (2.14) are non-dimensional variables from the fluid and structure solvers. Finally, the normal stress differences and the relative viscosity are defined by: where $\overline{\bullet }$ denotes the cycle-average value of $\bullet$ , and $(x,y,z)$ correspond to vorticity, flow and velocity gradient directions, respectively, as shown in figure 1. For a given suspension, these functions can completely describe the rheological behaviour that can be obtained from bulk stress components (Dealy & Wang Reference Dealy and Wang2013).

Figure 1.

2-D schematic of the computational set-up. The simulations are in 3-D.

2.5. Computational set-up and validations

A schematic of the computational set-up is shown in figure 1. Simulations start by placing the centre of mass of the particle in the centre of the computational domain, as depicted in the figure. To accelerate the steady state, we initialized the velocity field of the domain with $u_{z}=\dot{{\it\gamma}}y$ , which results in a shear field as follows:

where $\unicode[STIX]{x1D640}^{\infty }$ and ${\it\bf\Omega}^{\infty }$ are the strain-field tensor and rotation vector of undisturbed flow. The boundary condition in the velocity gradient direction ( $y$ -direction in figure 1) is a moving wall with velocity $\dot{{\it\gamma}}L_{y}/2$ in opposite directions, while the other directions are set to be periodic. The implementation of periodic boundary conditions creates a planar array of particles that rotate in phase, i.e. a periodic suspension. Such simulations can give a good approximation of a homogeneous suspension in low concentrations in which particles are sufficiently far from each other and close-distance interactions are less likely to happen. Note that, in a real suspension, particles are located randomly and rotate with a phase difference under shear, but the particle–particle interactions are low in the dilute and semi-dilute regimes. We found that the maximum fluctuations of the velocity field at the periodic boundaries are less than 5 % of the characteristic velocity $\dot{{\it\gamma}}a$ for the highest volume fraction ${\it\phi}=0.112$ and $\mathit{Re}=0.01$ . This shows that hydrodynamic interactions of particles at the periodic boundaries do not require a near-contact strategy and are not that significant. In addition, this boundary condition is a better approximation than the far-field boundary condition that is considered in the analytical Stokes solution for dilute regimes. The dimensions of the computational cuboid are $[L_{x}\times L_{y}\times L_{z}]$ , which vary from $[0.46a\times 2.5a\times 2.5a]$ to $[5a\times 5a\times 5a]$ in order to change volume fractions of suspensions for different aspect ratios.

The volume fraction of particles in suspension studies can range from very dilute to highly concentrated; for each, different assumptions and strategies are needed. At very low concentrations (dilute assumption) particles are sufficiently well separated such that any interactions between them can be assumed negligible. At slightly higher volume fractions (semi-dilute regime) particles are close enough such that they can interact hydrodynamically, but still far enough away to squeeze the flow in the small gaps between approaching particles – which can cause local viscous heating and lubrication or the formation of chains and networks of touching particles (Mueller et al. Reference Mueller, Llewellin and Mader2010). When the concentration increases, close-distance particle interactions and particle collisions need to be taken into account, which requires collision detection and modelling of the lubrication forces (Nguyen & Ladd Reference Nguyen and Ladd2002). To avoid such difficulties, this study is aimed at dilute and semi-dilute suspensions, where particle–particle contact is not likely to happen. In addition, the Stokes number which governs the collision dynamics (Davis, Serayssol & Hinch Reference Davis, Serayssol and Hinch1986) in our simulations is always less than or equal to 10, making viscous forces dominant and the particle inertia too weak to influence the collision (Joseph et al. Reference Joseph, Zenit, Hunt and Rosenwinkel2001). Our simulations are performed for $0.008\leqslant {\it\phi}\leqslant 0.112$ , for prolate ellipsoids with varying aspect ratio $2\leqslant \mathit{Ar}\leqslant 5$ . Following Doi & Edwards (Reference Doi and Edwards1978), for a given $N$ particles of length $L$ and diameter $D$ in unit volume, the dilute and semi-dilute regimes are defined as $N\ll 1/L^{3}$ and $1/L^{3}<N<1/L^{2}D$ , respectively. For an ellipsoidal particle these limits in non-dimensional form read as ${\it\phi}\ll ({\rm\pi}/6)\mathit{Ar}^{-2}$ and $({\rm\pi}/6)\mathit{Ar}^{-2}<{\it\phi}({\rm\pi}/6)<\mathit{Ar}^{-1}$ for the dilute and semi-dilute regimes, respectively. For instance, a suspension consisting of ellipsoids with aspect ratio $\mathit{Ar}=3$ is considered dilute if ${\it\phi}\ll 0.058$ , and semi-dilute if $0.058<{\it\phi}<0.175$ . All of the performed simulations in our study fall into the dilute and semi-dilute regimes.

The effect of Brownian motion is ignored in this work because we are focusing on the effect of inertia on the suspension under similar conditions. Furthermore, Brownian motion mainly affects the orientation distribution function of a suspension, which is irrelevant in a periodic suspension. The impact of Brownian motion is governed by the Péclet number $\mathit{Pe}=6{\rm\pi}{\it\mu}a^{3}\dot{{\it\gamma}}/kT$ , where $k=1.38\times 10^{-23}~\text{J}~\text{K}^{-1}$ is the Boltzmann constant, ${\it\mu}$ the dynamic viscosity of suspending liquid, and $T$ the absolute temperature. The Péclet number measures the strength of shearing to Brownian motion, and experiments show that the rheology of a suspension becomes independent of this parameter for $\mathit{Pe}>10^{3}$ (Stickel & Powell Reference Stickel and Powell2005). However, even at high $\mathit{Pe}$ in the limiting case of very high aspect ratios (slender particles), the dependence of the orientation of a particle on its initial orientation is much reduced (Hinch & Leal Reference Hinch and Leal1972, Reference Hinch and Leal1976; Petrie Reference Petrie1999). In this study the Peclet number is very high ( $\mathit{Pe}\approx O(10^{5})$ , even at $\mathit{Re}=0.01$ and a particle dimension of a few microns) and the particle aspect ratio is $\mathit{Ar}\leqslant 5$ . Therefore, the effect of Brownian motion can be neglected.

In the absence of Brownian motion and particle interactions, particles with an axis of symmetry move in orbits (as will be explained in detail in § 3) which are independent of time and determined just by the initial orientations. Therefore, the initial orientation of the ellipsoids plays a very important role in the particle’s motion and stress, and consequently on the relative viscosity of the suspension. For instance, if the axis of symmetry is aligned with the vorticity vector, the particle rolls along in the flow with a steady motion (Mueller et al. Reference Mueller, Llewellin and Mader2010), i.e. similar to a sphere. Therefore, the angular velocity and particle stress remain constant during the shear motion. In contrast, if the axis of revolution of the ellipsoid is perpendicular to the vorticity vector ( $x$ -axis), the aspect ratio of the ellipsoid has the maximum effect on the motion and particle stress. We used this orientation for our ellipsoids to observe the maximum variation in the rheological properties of the suspension. For any other initial orientation, the effect is intermediate between the spherical and presented ellipsoidal particles, depending on the angle between the axis of revolution of the ellipsoid and the vorticity vector.

We carried out extensive validation and grid sensitivity studies to verify the reliability of our solutions – see appendix A for validation of calculations using spherical particles. For ellipsoidal particles, we used a fine surface mesh with approximately 5000–15 000 triangular elements, depending on surface area, to satisfy the required accuracy for boundary condition reconstruction as mentioned in appendix A. We also tested the grid sensitivity of the numerical results by changing the grid size from ${\it\Delta}=0.04a$ to ${\it\Delta}=0.01a$ and compared the shear stresslet of particles. Figure 2 shows the results on the three grids at $\mathit{Re}=1.0$ and ${\it\phi}=0.118$ . We found the angular velocities computed on the three grids are almost identical, and the shear stresslet becomes insensitive to the grid size for ${\it\Delta}=0.02a$ and finer. In fact, the average values of the stresslet for the fine grid ${\it\Delta}=0.01a$ and the medium grid ${\it\Delta}=0.02a$ are very close, while for the coarse grid it is approximately 4 % higher than for the other grids. This indicates that the results are not sensitive to the grid size on the medium grid with resolution ${\it\Delta}=0.02a$ which is used for all the simulations reported in this work.

Figure 2.

Grid sensitivity: the shear stresslet for a prolate ellipsoid of aspect ratio $\mathit{Ar}=2$ in a suspension with volume fraction ${\it\phi}=0.112$ at $\mathit{Re}=1.0$ plotted during one cycle for different grid sizes.

4.1. Stresslet

The dominant part of a hydrodynamic stress system in the suspension due to the presence of particles is the symmetric part of the first moment of stress calculated on the surface of particles, called the stresslet. Based on (2.13) and (2.15), the contribution of the stresslet to the normal stress differences and the relative viscosity, $N_{1}$ , $N_{2}$ and ${\it\eta}_{r}$ are proportional to $\unicode[STIX]{x1D61A}_{zz}-\unicode[STIX]{x1D61A}_{yy}$ , $\unicode[STIX]{x1D61A}_{yy}-\unicode[STIX]{x1D61A}_{xx}$ and $\unicode[STIX]{x1D61A}_{yz}$ , respectively. To identify the role of inertia and verify our results, we compare our results against the analytical solution of a single ellipsoid at $(\mathit{Re},{\it\phi})\rightarrow 0$ (Kim & Karrila Reference Kim and Karrila2013) – see appendix B.

In figure 5, the Reynolds-scaled stresslet for a ellipsoid with aspect ratio $\mathit{Ar}=2$ in a suspension of volume fraction ${\it\phi}=0.058$ at various Reynolds numbers is compared with the analytical solution during one complete rotation about its centre of mass. The numerical stresslet calculation is verified by comparing it to the analytical calculation in figure 5, as it follows the same trend, has similar phases and amplitudes as the analytical values, and the deviations are the results of inertia, which increase as Re increases. The influence of inertia on the stresslet can easily be seen in figure 5. It can be observed that, as inertia increases, the shear stresslet $\mathit{Re}\unicode[STIX]{x1D61A}_{yz}$ and the second normal stresslet difference $\mathit{Re}(\unicode[STIX]{x1D61A}_{yy}-\unicode[STIX]{x1D61A}_{xx})$ increase, whereas the first normal stresslet difference $\mathit{Re}(\unicode[STIX]{x1D61A}_{zz}-\unicode[STIX]{x1D61A}_{yy})$ decreases during the cycle. However, the inertial effects on the these components of the stresslets are not uniform. For instance, deviation from the analytical solution in the shear stresslet is higher in some phases, or in the normal stresslet differences, inertia increases or decreases the stresslet depending on the orientation of particles. This asymmetrical growth is due to the fact that inertia acts differently on each component of the stresslet: the strain-field stresslet $\unicode[STIX]{x1D64E}^{\unicode[STIX]{x1D640}}$ , and the rotational stresslet $\unicode[STIX]{x1D64E}^{{\it\bf\Omega}}$ .

Figure 5.

Instantaneous shear stresslet (a) and normal stresslet differences (b) and (c), as a function of orientation for a prolate ellipsoid of aspect ratio $\mathit{Ar}=2$ in a suspension with volume fraction ${\it\phi}=0.058$ at various Reynolds numbers are compared with the exact solution of Kim & Karrila (Reference Kim and Karrila2013).

The particle’s stresslet is generated by two mechanisms based on the analytical solution of a single prolate ellipsoid in the linear flow field, i.e. the superposition of a rotation ${\it\bf\Omega}^{\infty }$ and a strain field $\unicode[STIX]{x1D640}^{\infty }$ (Kim & Karrila Reference Kim and Karrila2013) – see appendix B and (B 1):

where superscripts $\unicode[STIX]{x1D640}$ and ${\it\bf\Omega}$ denote the strain and the rotational parts of the stresslet, respectively. The right-hand side of (4.1c ) for a prolate ellipsoid rotating in the $C=\infty$ orbit (shown in figure 4) is given by (Kim & Karrila Reference Kim and Karrila2013): where $F_{i}$ are functions of the aspect ratio $\mathit{Ar}$ , and ${\it\theta}$ is measured as shown in figure 1. We will focus on the influence of inertia on the rotation and strain-field stresslet mechanisms to explain the variation of shear and normal stresslets with Reynolds number.

4.2. Acceleration stress

We discussed in § 4.1 how inertia affects the total particle stress of the suspension (by increasing the shear stresslet $\unicode[STIX]{x1D61A}_{yz}$ and the second normal stresslet difference $\unicode[STIX]{x1D61A}_{yy}-\unicode[STIX]{x1D61A}_{xx}$ and decreasing the first normal stresslet difference $\unicode[STIX]{x1D61A}_{zz}-\unicode[STIX]{x1D61A}_{yy}$ of the suspension). However, this is not the sole effect of inertia on the rheology of suspension. Based on (2.13), inertia plays a direct role in the total particle stress by generating acceleration and Reynolds stresses, which will be discussed in more detail in this subsection and § 4.3. Acceleration stress, by definition, is generated by the acceleration of particles moving within the suspension ((2.12) and (2.14b )). Similar to § 4.1, we have derived the analytical expression based on the Stokes solution of the problem (appendix C) to verify our numerical results, have a better understanding of the nature of this term, and establish the effect of inertia by comparing it against our solutions. It is worth noting that the particle acceleration is calculated based on the rotation of the particle, and translational acceleration has been ignored due to the fact that the centre of mass of the particle is located on the zero surface of shear, and the net of the force exerted on the particle is negligible. In appendix C, we obtain a formula for the acceleration stress of any arbitrarily shaped particle based on the Stokesian solution of Jeffery (Reference Jeffery1922) (assuming a particle rotating in the $C=\infty$ orbit). For our ellipsoidal particles, this formulation can be simplified and expressed in terms of the volume fraction, aspect ratio and angular kinematics of the particle as:

where ${\it\alpha}_{x}$ is the angular acceleration in the vorticity direction and ( $\unicode[STIX]{x1D617}_{yz}$ , $\unicode[STIX]{x1D617}_{zz}-\unicode[STIX]{x1D617}_{yy}$ , $\unicode[STIX]{x1D617}_{yy}-\unicode[STIX]{x1D617}_{xx}$ ) will be referred to as the shear acceleration stress, first normal acceleration stress difference, and second normal acceleration stress difference, respectively. In the remainder of this subsection we discuss the effect of the volume fraction, aspect ratio and inertia on the shear and normal components of the instantaneous acceleration stress.

Figure 6.

Shear acceleration stress as a function of orientation for a suspension of (a) prolate ellipsoids with aspect ratio $\mathit{Ar}=2$ at different Reynolds numbers and volume fraction ${\it\phi}=0.058$ , (b) ellipsoids with various aspect ratios at $\mathit{Re}=1$ and volume fraction ${\it\phi}=0.058$ , and (c) ellipsoids with aspect ratio $\mathit{Ar}=2$ and $\mathit{Re}=1$ at different volume fractions.

4.3. Reynolds stress

The last term on the right-hand side of (2.11) is the stress generated due to disturbances caused by the motion of particles in the flow. Mathematically, the Reynolds stress originates from the nonlinear (inertial) term of the Navier–Stokes equations. Physically speaking, interactions between the particles themselves and the flow cause velocity fluctuations for both the particles and the flow. Reynolds stress, as defined in (2.14c ), is the volume average of the velocity fluctuation products for the entire suspension volume ( $V=V_{f}+V_{p}$ ). Yeo & Maxey (Reference Yeo and Maxey2013) considered just the fluid part in the formulation of particle stress and ignored this term for the range $\mathit{Re}\leqslant 2.0$ . Kulkarni & Morris (Reference Kulkarni and Morris2008) considered both parts for the Reynolds stress and reported a linear increase of this term with respect to $\mathit{Re}$ . Both studies concluded that, for spherical particles, the Reynolds stress is negligible for flow regimes in the range $\mathit{Re}\sim O(1)$ . Recently, Haddadi & Morris (Reference Haddadi and Morris2014) investigated a dense suspension of spheres and found that the Reynolds stress influences the total second normal stress considerably when $\mathit{Re}\sim O(1)$ . Here, we calculate the Reynolds stress considering both fluid and particle fluctuations, i.e.  $\unicode[STIX]{x1D64D}=\unicode[STIX]{x1D64D}_{p}+\unicode[STIX]{x1D64D}_{f}$ . In order to calculate the velocity fluctuations ( $\boldsymbol{u}^{\prime }=\boldsymbol{u}-\mathbf{\mathscr{U}}$ ), the average velocity $\mathbf{\mathscr{U}}$ is considered as the shear velocity of the undisturbed flow, as described in (2.16a,b ) and shown in figure 1. The particle’s instantaneous velocity at each spatial point is calculated based on the vector position and angular velocity vector of the particle. In this section, $\unicode[STIX]{x1D619}_{yz}$ , $\unicode[STIX]{x1D619}_{zz}-\unicode[STIX]{x1D619}_{yy}$ and $\unicode[STIX]{x1D619}_{yy}-\unicode[STIX]{x1D619}_{xx}$ will be referred to as the shear Reynolds stress, first normal Reynolds stress difference and second normal Reynolds stress difference, respectively.

5.1. Intrinsic viscosity

As clearly seen in table 1, and previously shown in figures 6 and 9, the cycle averages of the shear acceleration $\overline{\unicode[STIX]{x1D617}_{yz}}$ and the shear Reynolds stress $\overline{\unicode[STIX]{x1D619}_{yz}}$ are zero, i.e. the acceleration stress and the Reynolds stress are symmetric parameters about zero (because of the symmetrical particle shape and symmetrical set-up). Therefore, the cycle average of the shear particle stress $\overline{{\it\Sigma}_{yz}}$ , which determines the relative viscosity of a suspension, is equal to the cycle average of the shear stresslet per unit volume $\overline{\unicode[STIX]{x1D61A}_{yz}}/V$ . Comparing the average value of the shear stresslet $\overline{\unicode[STIX]{x1D61A}_{yz}}/V$ with its standard deviation ${\it\sigma}_{\unicode[STIX]{x1D61A}_{yz}}/V$ , one can see that the relative average variation of the shear stresslet (ratio of standard deviation with respect to the average value) increases as the aspect ratio increases. For $\mathit{Ar}=2$ , ${\it\sigma}_{\unicode[STIX]{x1D61A}_{yz}}/V$ is approximately 10 % of the average value, while this ratio is approximately 20 % for $\mathit{Ar}=3$ , and more than 30 % when $\mathit{Ar}=5$ . Therefore, the higher the aspect ratio of the particles in a suspension, the higher these oscillations. On the other hand, for the semi-dilute regime ${\it\phi}>0.05$ and $\mathit{Re}\sim O(1)$ and higher, the variation of the shear acceleration and Reynolds stress terms ( ${\it\sigma}_{\unicode[STIX]{x1D617}_{yz}}$ and ${\it\sigma}_{\unicode[STIX]{x1D619}_{yz}}$ ) become important when compared to ${\it\sigma}_{\unicode[STIX]{x1D61A}_{yz}}/V$ . For instance, at $\mathit{Re}=10$ and $\mathit{Ar}=2$ , the average variation of the acceleration stress is approximately 50 % of the variation of the stresslet. Therefore, we can conclude that for suspensions with particle Reynolds number up to $\mathit{Re}\sim O(0.1)$ , regardless of the aspect ratio of the particle, the acceleration stress and Reynolds stress can be neglected. However, for higher Reynolds numbers ( $\mathit{Re}\geqslant 1$ ), the variation of the acceleration stress becomes important and is no longer negligible. This effect decreases slightly on increasing the aspect ratio of the particle, i.e. ellipsoids with $\mathit{Ar}=3$ and $\mathit{Ar}=5$ have a lower ratio of ${\it\sigma}_{\unicode[STIX]{x1D617}_{yz}}$ with respect to ${\it\sigma}_{\unicode[STIX]{x1D61A}_{yz}}/V$ , when compared to $\mathit{Ar}=2$ .

Table 1.

Cycle-average contribution of stresslet, acceleration stress and Reynolds stress in total shear particle stress for various aspect ratio $\mathit{Ar}$ , particle Reynolds number $\mathit{Re}$ and suspension volume fraction ${\it\phi}$ .

Since our simulations are at the limit of the dilute regime, we can ignore the second term of (1.1) and focus on the intrinsic viscosity $[{\it\eta}]$ . As mentioned in § 1, this constant for non-Brownian suspensions depends on the shape, size and orientation of the particles, and can be determined either analytically for simple shapes (Jeffery Reference Jeffery1922; Guth & Gold Reference Guth and Gold1938; Einstein Reference Einstein1956; Leal & Hinch Reference Leal and Hinch1971; Burgers Reference Burgers1995) or experimentally for more general geometries and concentrations (Jeffrey & Acrivos Reference Jeffrey and Acrivos1976; Mueller et al. Reference Mueller, Llewellin and Mader2010). Based on the cycle-average values of components of the shear particle stress presented in table 1, the intrinsic viscosity can be calculated as follows:

Figure 12.

Intrinsic viscosity of a suspension as a function of (a) Reynolds number for different aspect ratios and volume fractions, (b) volume fraction for prolate ellipsoids with aspect ratio $\mathit{Ar}=2$ and (c) aspect ratio for constant volume fraction ${\it\phi}=0.058$ .

Using the cycle-average shear stresslet, we present the intrinsic viscosity of the suspension $[{\it\eta}]$ in the last column of table 1, and plot the results in figure 12. The effect of $\mathit{Re}$ on the intrinsic viscosity is shown in figure 12(a), in which $[{\it\eta}]$ is plotted as a function of $\mathit{Re}$ for different particle shapes and volume fractions. As can be observed in the figure, the intrinsic viscosity of the suspension consistently increases as the inertia increases. The highest increase occurs for a very dilute suspension ( ${\it\phi}=0.008$ ) of ellipsoids with $\mathit{Ar}=2$ , which is approximately 23 % when inertia increases from $\mathit{Re}=0.01$ to $\mathit{Re}=10$ . Other cases experience an 11–14 % increase of $[{\it\eta}]$ for the same variation of Reynolds number. This increase is a consequence of increasing the Reynolds-scale shear stresslet of the particles due to the strain field $\mathit{Re}\unicode[STIX]{x1D61A}_{yz}^{\unicode[STIX]{x1D640}}$ , which was discussed in § 4.1.1, and is in agreement with the numerical results of Mikulencak & Morris (Reference Mikulencak and Morris2004) and Yeo & Maxey (Reference Yeo and Maxey2013) and analytical prediction of Lin et al. (Reference Lin, Peery and Schowalter1970) and Subramanian et al. (Reference Subramanian, Koch, Zhang and Yang2011) for spherical particles.

The effect of suspension concentration on the intrinsic viscosity can be seen in figure 12(b), which shows the intrinsic viscosity of suspensions of ellipsoids with constant aspect ratio $\mathit{Ar}=2$ as a function of volume fraction at various Reynolds numbers. The analytical value of the intrinsic viscosity (which ignores inertia) for this aspect ratio, i.e.  $[{\it\eta}]=2.819$ (Jeffery Reference Jeffery1922), is also plotted for comparison with the numerical results. It can be observed that $[{\it\eta}]$ increases uniformly with volume fraction in the range $0.01\leqslant \mathit{Re}\leqslant 1.0$ , while at $\mathit{Re}\sim O(10)$ the slope of variation becomes nonlinear and small. Comparing figures 12(a) and 12(b), one can observe the percentage changes for the intrinsic viscosity with respect to inertia and volume fraction are very close and comparable. The combined effects of inertia and concentration are also clear in figure 12(b), where the higher the ${\it\phi}$ and $\mathit{Re}$ , the higher the deviation from the Stokesian solution.

The final investigated parameter which influences the intrinsic viscosity is the aspect ratio of the ellipsoids. In figure 12(c), $[{\it\eta}]$ is shown as a function of aspect ratio $\mathit{Ar}$ at constant volume fraction ${\it\phi}=0.058$ for various Reynolds numbers and compared with the Jeffery’s Stokes solution at $(\mathit{Re},{\it\phi}\rightarrow 0)$ and the empirical fit of Mueller et al. (Reference Mueller, Llewellin and Mader2010). Based on the Stokesian solution (as clearly shown in the figure), suspensions consisting of ellipsoids with a higher aspect ratio have higher relative viscosity. However, the rate of increase of $[{\it\eta}]$ with respect to $\mathit{Ar}$ is higher than the analytical prediction. In fact, even in the very low inertia regime ( $\mathit{Re}=0.01$ ) and dilute concentration ( ${\it\phi}=0.058$ ), Jeffery’s solution underestimates the intrinsic viscosity by approximately 12 % when the aspect ratio of the ellipsoid is $\mathit{Ar}=5$ . Another interesting fact about this figure is that the influence of the aspect ratio on the intrinsic viscosity of a suspension is significant. For instance, on changing the aspect ratio from $\mathit{Ar}=2$ to $\mathit{Ar}=5$ , $[{\it\eta}]$ increases by approximately 30 % for the whole range of Reynolds numbers.

The most successful experimental data fitting for ${\it\eta}_{r}={\it\eta}_{r}({\it\phi})$ , considering the entire range of volume fraction (from very dilute to highly concentrated), is the model proposed by Krieger & Dougherty (Reference Krieger and Dougherty1959), which considers the maximum packing volume fraction ${\it\phi}_{m}$ as a fitting parameter:

In contrast to the polynomial correlation suggested for the dilute and semi-dilute regimes, i.e. (1.1), (5.2) is able to predict the behaviour of the relative viscosity for the entire range of volume fraction, especially close to the packing regime where ${\it\eta}_{r}\rightarrow \infty$ as ${\it\phi}\rightarrow {\it\phi}_{m}$ . A simpler relation suggested by Maron & Pierce (Reference Maron and Pierce1956) can be used when setting $[{\it\eta}]{\it\phi}_{m}=2$ . Using this correlation between the intrinsic viscosity and the packing volume fraction, the binomial expansion of both relations would be equivalent to (1.1) for the dilute and semi-dilute regimes when ${\it\phi}\rightarrow 0$ . Mueller et al. (Reference Mueller, Llewellin and Mader2010) found the Maron & Pierce (Reference Maron and Pierce1956) fit to be very accurate for particles with aspect ratios in the range $\mathit{Ar}\in (1,9)$ and suggested the empirical relation $[{\it\eta}]=3.02+0.321\mathit{Ar}$ . This experimental data fit is shown and compared with our results for intrinsic viscosity in figure 12(c). Our results are generally in good agreement with the empirical fit, and a higher value of the experimental fit is somehow expected. Pabst et al. (Reference Pabst, Gregorová and Berthold2006) explain this overestimate of the experimental intrinsic viscosity as the difference between assumptions of intrinsic viscosity (perfectly homogeneous and stable) in the dilute regime, and systematic and statistical errors associated with non-dilute measurements. Note that the empirical fits of Mueller et al. (Reference Mueller, Llewellin and Mader2010) based on the Krieger & Dougherty (Reference Krieger and Dougherty1959) and Maron & Pierce (Reference Maron and Pierce1956) models are obtained by considering the entire range of ${\it\phi}\in (0,{\it\phi}_{m})$ , not just the small range in volume fraction that corresponds to dilute and semi-dilute regimes, where ${\it\phi}\ll {\it\phi}_{m}$ .

Finally, considering all the parts of figure 12, we can conclude that the inertia, volume fraction and aspect ratio of the particles are all influential parameters on the intrinsic viscosity of suspensions consisting of ellipsoids. Even though an increase in the viscosity of the suspension with these variables could have been expected separately, our numerical results show the synergy of these influential parameters for various cases. We found that the inertia and volume fraction reinforce the effects of each other to increase the viscosity of the suspension, and this increase becomes more important at higher aspect ratios. To show the synergistic effect of $\mathit{Re}$ , ${\it\phi}$ and $\mathit{Ar}$ on $[{\it\eta}]$ , for example, we can see that the intrinsic viscosity of a suspension can increase by 55 % when $(\mathit{Re},{\it\phi},\mathit{Ar})$ change from $(0.01,0.008,2)$ to $(10,0.058,5)$ , while the Stokesian solution predicts the increase to be just 18 %.

5.2. Normal stress differences

The normal stress differences, $N_{1}$ and $N_{2}$ , indicate the non-isotropic nature of the diagonal components of the stress tensor, and play an important role in particle migration. The average values and the corresponding standard deviations of the components of these variables, which are expected to be zero in the absence of non-hydrodynamic forces, are calculated and shown in tables 2 and 3. The major component of these particle-oriented variables, as discussed extensively in § 4.1, are the stresslet parts ( $\unicode[STIX]{x1D61A}_{zz}-\unicode[STIX]{x1D61A}_{yy}$ and $\unicode[STIX]{x1D61A}_{yy}-\unicode[STIX]{x1D61A}_{xx}$ ), which oscillate between positive and negative values. Particles’ stresslets, as can be seen in (B 1) and figure 5, are inversely proportional to the Reynolds number. This dependency is clearer in the average variation of the stresslet contribution to the normal stress differences, i.e.  $(1/V){\it\sigma}_{(\unicode[STIX]{x1D61A}_{zz}-\unicode[STIX]{x1D61A}_{yy})},(1/V){\it\sigma}_{(\unicode[STIX]{x1D61A}_{yy}-\unicode[STIX]{x1D61A}_{xx})}\approx K(\mathit{Ar},{\it\phi})/\mathit{Re}$ – see tables 2 and 3. However, these oscillations around zero are not exactly symmetric, and the asymmetric behaviour increases with inertia (figure 5 b,c). Consequently, the cycle-average value of the normal stresslet differences becomes larger when the inertia decreases. In contrast to the standard deviation of the normal stresslet difference, the slope of variation of $(1/V)\overline{\unicode[STIX]{x1D61A}_{zz}}-\overline{\unicode[STIX]{x1D61A}_{yy}}$ and $(1/V)\overline{\unicode[STIX]{x1D61A}_{yy}}-\overline{\unicode[STIX]{x1D61A}_{xx}}$ with Reynolds number is nonlinear and complicated, because on the one hand the normal stresslet differences become more symmetric (figure 5 b,c), and on the other hand the amplitude of these variables increases rapidly with inertia. The latter effect is obviously more drastic, as can be seen in the cycle-average values of the normal stresslet differences in tables 2 and 3, and the average values decrease monotonically as inertia decreases. The signs of both $(1/V)\overline{\unicode[STIX]{x1D61A}_{zz}}-\overline{\unicode[STIX]{x1D61A}_{yy}}$ and $(1/V)\overline{\unicode[STIX]{x1D61A}_{yy}}-\overline{\unicode[STIX]{x1D61A}_{xx}}$ remain unchanged and consistent for $\mathit{Re}\geqslant 0.1$ , i.e. the first normal stresslet difference is negative and the second normal stresslet difference is positive.

Table 2.

Cycle-average contribution of stresslet, acceleration stress and Reynolds stress in the first normal stress differences for various aspect ratio $\mathit{Ar}$ , particle Reynolds number $\mathit{Re}$ , and suspension volume fraction ${\it\phi}$ .

Table 3.

Cycle-average contribution of stresslet, acceleration stress and Reynolds stress in the second normal stress differences for various aspect ratio $\mathit{Ar}$ , particle Reynolds number $\mathit{Re}$ , and suspension volume fraction ${\it\phi}$ .

The cycle average and standard deviation of the other two components of the normal stress differences, acceleration and Reynolds stress parts change very little with Reynolds number for constant volume fraction and aspect ratio of particles, as can be observed in tables 2 and 3. In fact, what makes these components important at a moderate inertial regime is that their values remain almost constant, while the stresslet part decreases rapidly with Reynolds number.

For the first normal stress difference, the cycle-average acceleration stress, $\overline{\unicode[STIX]{x1D617}_{zz}}-\overline{\unicode[STIX]{x1D617}_{yy}}$ , is zero in all cases and has no contribution in $N_{1}$ , whereas the Reynolds stress becomes influential for $\mathit{Re}\geqslant 0.1$ . Comparing values of $N_{1}$ and $\overline{\unicode[STIX]{x1D619}_{zz}}-\overline{\unicode[STIX]{x1D619}_{yy}}$ in table 2, one can clearly observe the effect of the Reynolds stress in the calculation of the first normal stress difference. In fact, in some cases, the contribution of the Reynolds number stress is equal to or even greater than the stresslet in the first normal stress difference. Since the value of Reynolds stress is positive for all cases (figures 10 and 11), neglecting this term would result in an overestimation of $N_{1}$ , which is negative. Considering the variation of these terms about the cycle-average values, the acceleration stress would become as important as the Reynolds stress components, i.e.  ${\it\sigma}_{(\unicode[STIX]{x1D617}_{zz}-\unicode[STIX]{x1D617}_{yy})}\approx {\it\sigma}_{(\unicode[STIX]{x1D619}_{zz}-\unicode[STIX]{x1D619}_{yy})}$ . The negative values of $N_{1}$ in table 2 are in agreement with the Stokesian results of Kumar & Higdon (Reference Kumar and Higdon2011b ) for dicolloidal particles at comparable volume fractions, even though our results are one order of magnitude larger. This difference in magnitude is probably due to a lower value of dicolloidal particle’s aspect ratio, neglecting inertial effects, and considering different orientations for the particles in their simulation, whereas we only discuss ellipsoids for a specific initial orientation ( $C=\infty$ orbit).

For the second normal stress difference, cycle averages of the acceleration and Reynolds stresses are of comparable magnitude but opposite signs. Therefore, they cancel out the effects of each other on the cycle average of the second normal stress difference. Comparing the averages of the three components in table 3 reveals that the contributions of different mechanisms of particle stress are comparable at $\mathit{Re}\geqslant 1.0$ . However, as the Reynolds stress and acceleration stress have opposite average values, their influence cannot be seen in the second normal stress difference, i.e.  $N_{2}$ and $(1/V)\overline{\unicode[STIX]{x1D61A}_{yy}}-\overline{\unicode[STIX]{x1D61A}_{xx}}$ are almost equal, not because $\overline{\unicode[STIX]{x1D617}_{yy}}-\overline{\unicode[STIX]{x1D617}_{xx}}$ and $\overline{\unicode[STIX]{x1D619}_{yy}}-\overline{\unicode[STIX]{x1D619}_{xx}}$ are negligible, but because they are equal and opposite. Nevertheless, similar to the first normal stress and shear stress, these inertial terms become important at small time scales, where the variation of these components about the average values need to be considered.

Finally, comparing the normal stress differences in tables 2 and 3 (except for very low inertial regimes ( $\mathit{Re}\sim O(0.01)$ )) $N_{1}$ has a negative value, while $N_{2}$ has a positive value with $|N_{1}|>|N_{2}|$ . This trend is in agreement with the analytical predictions of Lin et al. (Reference Lin, Peery and Schowalter1970) and Subramanian et al. (Reference Subramanian, Koch, Zhang and Yang2011) and the numerical simulations of Mikulencak & Morris (Reference Mikulencak and Morris2004), Kulkarni & Morris (Reference Kulkarni and Morris2008) and Haddadi & Morris (Reference Haddadi and Morris2014) for spherical particles at comparable volume fractions and inertia, as well as the Stokesian solution of Rallison (Reference Rallison1978) for ellipsoids with aspect ratios less than 5 in a Brownian dilute suspension. Note that, to the best of our knowledge, no experimental data are available on normal stress differences for direct comparison with our results in the range of volume fractions we have studied in this paper.