The dynamics of suspensions of solid particles in a fluid medium is inherently complex, yet such complex fluids are ubiquitous in nature, e.g. lavas, slurries and † Email address for correspondence: hormozi@ohio.edu D. Alghalibi, I. Lashgari, L . Brandt and S. Hormozi debris, and are important for a wide variety of industrial processes, such as paints, pastes, concrete casting, drilling muds, waste disposal, food processing, crude oil flows with rocks and medicine. This wealth clarifies why the behaviour of these suspensions has been extensively studied both experimentally, theoretically and, more recently, via numerical simulations. These studies have uncovered numerous complex features of suspension flows, which are, however, not yet fully understood (see Stickel & Powell 2005). The complexities and challenges are attributed to the large variety of interactions among particles (hydrodynamic, contact, interparticle forces), the physical properties of the particles (shape, size, deformability, volume fraction) and the properties of the suspending fluid (Newtonian or non-Newtonian). In this work we numerically study suspensions of neutrally buoyant rigid spheres in both Newtonian and inelastic non-Newtonian fluids.
The rheology of neutrally buoyant non-Brownian particles suspended in a Newtonian fluid has been largely investigated, see Batchelor (1970), Brady & Bossis (1988), Larson (1999), Stickel & Powell (2005). This is determined by two dimensionless numbers: the solid volume fraction and the particle Reynolds number Re p D f Pa 2 = (where a is the particle radius, f the fluid density, the fluid viscosity and P the flow shear rate). Many studies focused on the limiting case of Stokesian suspensions when inertia is negligible (i.e. Re p ! 0) and the effective viscosity of the suspension depends only on . Theoretical works mainly address the limiting cases of ! 0 and ! m a x , where max is the maximum packing fraction. When a suspension is dilute, its effective viscosity follows the linear behaviour derived by Einstein (1906Einstein ( , 1911) ) eff D .1 C 2:5/ (particle interactions are neglected) or the quadratic formulation of Batchelor (1977) eff D .1 C 2:5 C 6:95 / (with mutual particle interactions included). The recent theoretical works of Wyart and co-workers (see e.g. DeGiuli et al. 2015) make use of perturbations around the jamming state to show that the effective viscosity diverges as ! max .
At moderate to large values of , however, the rheology becomes more complex due to multi-body and short-range interactions. The suspension shear viscosity increases with before diverging at max . In addition, normal stresses appear when the suspension is subject to shear. A number of studies have been performed to measure the suspension effective shear viscosity and the normal stresses, see the experiments of Krieger & Dougherty (1959), Zarraga, Hill & Leighton (2000), Singh & Nott (2003), Ovarlez, Bertrand & Rodts (2006), Deboeuf et al. (2009), Bonnoit et al. (2010), Boyer, Guazzelli & Pouliquen (2011), Couturier et al. (2011), Dbouk, Lobry & Lemaire (2013b) and numerical simulations of Sierou & Brady (2002), Yurkovetsky & Morris (2008), Yeo & Maxey (2010), Dbouk et al. (2013a).
In the presence of weak inertia (Re p D 0) the rheological measurements start to differ from those in the Stokesian regime. For suspensions with 0:02 6 Re p 6 10 and 0:1 6 6 0:3, the recent numerical studies of Kulkarni & Morris (2008), Picano et al. (2013), Yeo & Maxey (2013) show an increase in the suspension stresses as Re p increases, although discrepancies exist in the reported values of the stresses. Due to the improvement of computational methods (see e.g. Kulkarni & Morris 2008;Yeo & Maxey 2011, 2013;Lashgari et al. 2014;Picano, Breugem & Brandt 2015;Fornari et al. 2016), interface-resolved simulations of solid particles in Newtonian fluids can now reveal details of the suspension microstructure, and shed light on the role of inertia in the overall dynamics (Prosperetti 2015). More precisely, the study of Picano et al. (2013) shows that the inertia affects the suspension microstructure, resulting in an enhancement of effective shear viscosity. As Re p increases, the pair distribution function becomes more anisotropic, almost zero at the rear of the particles, so-called excluded volumes. This increases the effective solid volume fraction, and consequently, the effective shear viscosity. Taking into account this excluded volume effect (which depends on both and Re p ), Picano et al. (2013) scaled the effective shear viscosity in the presence of inertia that is small compared to that of Stokesian suspensions. While there are hardly any studies addressing the rheology of suspensions for Re p > 1 0 and 0:1 6 6 0:45 (Bagnold 1954;Lashgari et al. 2014Lashgari et al. , 2016;;Linares-Guerrero, Hunt & Zenit 2017), there is a considerable body of research addressing the rheology of dry granular materials ( > 0:45), where viscous effects are negligible and friction, collision and particle phase momentum govern the flow dynamics (Trulsson, Andreotti & Claudin 2012;Andreotti, Forterre & Pouliquen 2013;DeGiuli et al. 2015;Amarsid et al. 2017).
The behaviour of suspensions is even more complex when the carrier fluid is non-Newtonian, such as generalised Newtonian or viscoelastic fluids. Only a few studies have been devoted to non-colloidal particles suspended in non-Newtonian fluids, attempting to address the bulk rheology from a continuum-level closure perspective. These studies mainly focus on non-inertial suspensions with few exceptions e.g. Hormozi & Frigaard (2017).
On the theoretical front, the homogenisation approach is adopted by Chateau, Ovarlez & Trung (2008) to derive constitutive laws for suspensions of non-colloidal particles in yield stress fluids. The authors consider a Herschel-Bulkley suspending fluid and show that the bulk rheology of suspensions also follows the Herschel-Bulkley model, with an identical power-law index, but with a yield stress and consistency that increase with the solid volume fraction. Generally, adding large particles to a fluid enhances both the effective viscosity of the bulk and the local shear rate of the fluid phase. While the latter has no influence on the viscosity of a Newtonian suspending fluid, it strongly influences the local apparent viscosity in the case of a non-Newtonian suspending fluid. In the homogenisation theory of Chateau et al. (2008) the mean value of the local shear rate is estimated via an energy argument and used to derive the suspension constitutive laws. A number of experimental works have been carried out (see Ovarlez et al. 2006;Chateau et al. 2008;Mahaut et al. 2008;Coussot et al. 2009;Vu, Ovarlez & Chateau 2010;Ovarlez et al. 2012;Dagois-Bohy et al. 2015) showing a general agreement with the homogenisation theory.
It has also been shown that adding large particles to a shear thinning fluid not only enhances the effective shear viscosity, but also promotes the onset of the non-Newtonian behaviour to a smaller shear rate (Poslinski et al. 1988;Liard et al. 2014). These features are also observed for shear thickening suspending fluids in both continuous shear thickening (CST) and discontinuous shear thickening (DST) scenarios (Cwalina & Wagner 2014;Liard et al. 2014;Madraki et al. 2017). The advancement in the onset of shear thinning, C S T or DST, is explained by the increase of the local shear rate in the suspending fluid due to the presence of larger particles (Chateau et al. 2008;DeGiuli et al. 2015). Recent studies show that the homogenisation theory underestimates the occurrence of DST (Madraki et al. 2017) and overestimates the values of the yield stress, which depends on the shear history and, consequently, on the suspension microstructure, see Ovarlez et al. (2015). While the former can be attributed to the fact that the mean value of the local shear rate is not sufficient to predict DST, instead controlled by the extreme values of the local shear rate distribution, the latter implies that constitutive laws derived from homogenisation theory need to be refined by taking into account the microstructure and shear history. Here, we present a numerical study of rigid-particle suspensions in pseudoplastic (via the Carreau model) and dilatant fluids (via the power-law model) employing an immersed boundary method. The method was originally proposed by Breugem (2012) to enable us to resolve the dynamics of finite-size neutrally buoyant particles in flows.
The simulations are performed in a planar Couette configuration where the rheological parameters of the suspending fluids and the particle properties are kept constant while the volume fraction of the solid phase varies in the range of 10 % 6 6 40 % and the bulk shear rate varies to provide particle Reynolds numbers in the range of 0:1 6 Re p 6 10. The well-resolved simulations benefit the study via providing access to the local values of particle and fluid phase velocities, shear rate and particle volume fraction. We explore the confinement effects, microstructure and rheology of these complex suspensions when the particle Reynolds number is non-zero. We compare our results with the recently proposed constitutive laws based on the homogenisation theory (see Chateau et al. 2008) and refine these laws for the case when inertia is present.
The paper is organised as follows. The governing equations and numerical method are discussed in § 2. We present our results in § 3, first considering the local distribution of flow profiles and the Stokesian rheology of non-colloidal suspensions with generalised Newtonian suspending fluids. We then discuss the role of inertia and propose new rheological laws including inertial effects. A summary of the main conclusions and some final remarks are presented in § 4. where u D .u; v; w/ is the velocity vector containing the spanwise, streamwise and wall-normal components corresponding to the .x; y; z/ coordinate directions respectively (see figure 1). The pressure is denoted by P while the density of both the fluid and particles is indicated by as we consider neutrally buoyant particles. The fluid viscosity O varies as a function of the local shear rate P .u/ following the rheological Carreau-law or power-law models defined below. Finally the body force f is added to the right-hand side of the equation to indicate how the no-slip condition at the particle surface is effectively implemented and it indicates the forcing from the dispersed phase on the carrier fluid.
The motion of the rigid spherical particles is described by the Newton-Euler equations, where U p and p are the translational and angular velocity of the particle p, while m p and I p are the mass and moment of inertia, 2m p a 2 =5, of a sphere with radius a. F p and T p are the net force and momentum resulting from hydrodynamic and particle-particle interactions,
In these equations @V p represents the surface of the particles with outwards normal vector n and I the identity tensor. The radial distance from the centre to the surface of each particle is indicated by r. The force and torque, F c and T c , act on the particle as a result of particle-particle or particle-wall contacts. The no-slip and no-penetration boundary conditions on the surface of the particles are imposed by forcing the fluid velocity at each point on the surface of the particle, X , to be equal to particle velocity at that point, u.X / D U p .X/ D U p C p r. This condition is not imposed directly in the immersed boundary method used in the current study, but instead included via the body force f on the right-hand side of (2.1).
Several models have been developed to capture the inelastic behaviour of some non-Newtonian fluids such as polymeric solutions. In the current work, we employ the Carreau law to describe the behaviour of shear thinning (pseudoplastic) fluids. In the expression above O is the non-dimensional viscosity O D = 0 , where 0 is the zero shear rate viscosity. In this work the non-dimensional viscosity takes the value O D 1 at zero shear rate and O D 1 = 0 D 0:001 in the limit of infinite shear rate, as shown in figure 2(a). The second invariant of the strain-rate tensor P is determined by the dyadic product of the strain tensor P D 2s ij V s ij , where s D . r u C r u T /=2 (see Bird et al. 1987). The power index n indicates the non-Newtonian fluid behaviour. For n < 1 the fluid is shear thinning where the fluid viscosity decreases monotonically with the shear rate. The constant is a dimensionless time, scaled by the flow time scale, and represents the degree of shear thinning. In the present study the power index and the time constant are fixed to n D 0:3 and D 10. We report in figure 2(a), using markers, the range of shear rate and the corresponding viscosity of the carrier fluid used for the different simulations of particle-laden flows discussed later.
For shear thickening (dilatant For a more detailed description of the parameters appearing in the Carreau-and powerlaw models we refer the readers to the book by Morrison (2001).
The particle motion in the flow is simulated by means of an efficient immersed boundary method (IBM) coupled with a flow solver for the generalised Navier-Stokes equations. We follow the formulation by Breugem (2012) developed to increase the numerical accuracy and stability for simulations of neutrally buoyant particles.
The governing equations for the fluid phase are discretised using a second-order central difference scheme; the Crank-Nicholson scheme is used for the time integration of the viscous term while the nonlinear term is treated explicitly using the three-step Runge-Kutta scheme. A fixed and staggered Eulerian grid is used for the fluid phase whereas a Lagrangian grid is attached to the surface of each particle. These two grid points communicate via the IBM forcing to satisfy the no-slip and no-penetration boundary conditions on the surface of the particles.
The interactions between the particles and/or wall are taken into account using the lubrication correction and soft collision model described in detail in Costa et al. (2015). When the gap distance between the particle-particle or particle-wall becomes smaller than a certain threshold, a mesh dependent lubrication correction based on the asymptotic solution by Brenner (1961) is employed to reproduce correctly the interaction between the particles. At smaller gaps the lubrication correction is kept constant to account for the surface roughness. Finally, a soft-sphere collision model is activated based on the relative velocity and the overlap between the two particles (particle-wall), where both the normal and tangential components of the contact force (2015), among others. One issue with the IBM scheme is the fluid trapped inside of the particles. Although, the effect of the fluid acceleration inside the particles is accounted for and subtracted from the forces acting on the particle, it can still result in non-realistic values of the local shear rate for the grid points close to the interface. To fix this issue, a volume of fluid scheme is used to create a velocity field that replaces the velocity of the fluid inside the particles with their rigid body motion.
At the interface (the surface of the particles), a weighted average of solid and fluid velocities is considered based on the local solid volume fraction of the Eulerian grid in this region.
Another issue concerns the use of the lubrication correction as suggested in the original IBM method of Breugem (2012) for Newtonian fluids. To extend this model to non-Newtonian fluids and approximate the lubrication correction force when particles are less than one grid cell away from each other, we use in the asymptotic solution by Brenner (1961) the local viscosity at the Eulerian point closest to the mid-point of the line connecting the centres of two particles in interaction. It is noteworthy to mention that the viscosity is calculated explicitly from the local shear rate, using the velocities from the previous substep of the Runge-Kutta time D. Alghalibi, I. Lashgari, L . Brandt and S. Hormozi Figure 3 shows a typical example of this validation.
In this study, we perform interface-resolved simulations of suspensions of neutrally buoyant spheres in shear thinning and shear thickening fluids. The flow is driven by the motion of the upper and lower walls in a plane Couette configuration. Periodic boundary conditions are imposed in the streamwise and spanwise directions. Similar to Picano et al. (2013), we use a box size of 2h 3:2h 3:2h with h the half-channel width; the number of uniform Eulerian grid points is 80 128 128 in the wall-normal, streamwise and spanwise directions. The particles have all the same radius, a D h=5, which corresponds to 8 Eulerian grid points per particle radius, whereas 746 Lagrangian grid points are used on the surface of each particle to resolve the fluidparticle interactions. The fluid is sheared in the y-z plane by imposing a constant streamwise velocity of opposite sign V w D Ph at the two horizontal walls.
In this work we fix the rheological parameters and vary the wall velocity (the shear rate P) and the particle volume fraction . We explore a wide range of shear rates, 0:1 6 P 6 10, corresponding to the particle Reynolds numbers 0:1 6 Re p D P a 2 =6 10 for the Newtonian fluid, 0:5 6 P 6 10 corresponding to 0:0624 6 Re p 6 9:8113 for the shear thinning fluid and 1 6 P 6 10 000 corresponding to 0:1 6 Re p 6 10 for the shear thickening fluid, see figure 2. Note that in the definition of the particle Reynolds number is, for each case, the non-Newtonian fluid viscosity in the absence of particles. For each suspending fluid, (Newtonian, shear thinning or shear thickening fluid), the range of shear rate as well as the values of zero shear rate viscosities are chosen such that we cover the same range of particle Reynolds number, i.e. 0:1 6 Re p 6 10. summarised in table 1. Results are collected after the flow reaches a statistically steady state. We ensure the convergence by repeating the analysis using half the number of samples and comparing the statistics with those from the entire number of samples.
In the present study, we investigate the flow of rigid particles in a simple shear where the carrier fluid is Newtonian, shear thinning or shear thickening. We focus on the bulk properties of the suspension as well as its local behaviour. We present the distribution of particle and fluid phase velocity, particle volume fraction and local shear rate. Then we study the rheology of these suspensions and compare our results with predictions from the homogenisation theory presented by Chateau et al. (2008), valid for Stokesian suspensions. We therefore focus on how inertia affects the suspension behaviour.
The wall-normal profiles of the local particle volume fraction .y/ are computed for all the simulations listed in table 1 by averaging the local solid volume fraction over time and in the spanwise direction. A typical example of the results is shown in figure 4. Here the wall-normal distribution of .y/ is displayed across half of the gap for four particle volume fractions, D T0:11; 0:21; 0:315; 0:4U, and two particle Reynolds numbers Re p D T0:1; 6U for Newtonian, shear thinning and shear thickening suspending fluids. Three features are evident here. First, particles tend to form layers due to the wall confinement; the number of particles at the wall is larger than in the bulk. Second, the particle layering increases as the bulk solid volume fraction increases. Third, the distribution of .y/ changes slightly with the type of suspending fluid over the range of the particle Reynolds number studied here (i.e. 0 < Re p 6 10), suggesting that the local particle volume fraction is mainly controlled by geometry and confinement. A detailed study of confinement is beyond the scope
0:5 0:5 0:5 0:5 0:5 0:5 0:5 0:5 0:5 of this manuscript. However, more recent results to be published elsewhere show that decreasing confinement by 50 % when dealing with oblate particles has a less than 1 % effect on the bulk rheology.
We report the normalised mean fluid V f =V w and particle velocity V p =V w in figure 5 for the same values of the particle Reynolds number and bulk solid volume fraction in figure 4, using the same symbol and colour scheme throughout the manuscript. The statistics of the fluid phase velocity have been computed neglecting the points occupied by the solid phase in each field (phase-ensemble average). Generally, independent of the type of suspending fluid, the normalised mean fluid velocity decreases in the intermediate region between the wall and the centreline, 0:02 . y=h . 0:4 as the particle volume fraction and shear rate or Re p increase. The larger differences between the different cases are found close to the wall. The deviation of the normalised mean fluid velocity profile from linearity is more pronounced for the shear thinning suspending fluid and less evident for the shear thickening fluid. This is due to the fact that the local viscosity seen by the particles in the case of generalised Newtonian fluids depends on the local shear rate. The statistics pertaining to the solid phase, depicted in figure 5(c,d), are calculated using quantities related to each individual particle, and taking the phase-ensemble average over time and space. As for the carrier fluid, the normalised mean particle velocity decreases when increasing the volume fraction and inertial effects, Re p . Again this is more significant for the case of shear thinning suspending fluid due to the smaller apparent viscosity at the particle scale (or equivalently the larger local particle Reynolds number). Moreover, comparing figures 5(a) and 5(c) (and similarly figures 5b and 5d), we note that the slip velocity between the solid phase and fluid phase increases close to the wall. In particular, particles move faster close to the walls, y=h . 0:2, something explained by the different boundary conditions (particles can roll and slide on the wall). This slip is more evident for the case of a shear thinning suspending fluid, suggesting that as we increase the local particle Reynolds number, any modelling would need to take into account a boundary layer close to the wall governed by a two-phase equation of motion instead of mixture equations (Dontsov & Peirce 2014;Costa et al. 2016). Indeed, continuum models are more prone to failure when the slip velocity between fluid and particles is not negligible, i.e. the particle Stokes number is large.
In this section, we briefly explain the theoretical prediction for the rheology of an inertialess suspension of rigid spherical particles in generalised Newtonian fluids. This theoretical approach is based on homogenisation theory and was first developed by Chateau et al. (2008) In order to predict the rheology of the suspensions we need to know the value of the local shear rate P local .x; y; z/ for a bulk shear rate P. In homogenisation theory, it is assumed that the bulk rheology is determined by the mean value of the local shear rate, i.e. P local . Following Chateau et al. (2008), we assume that viscous dissipation is responsible for the entire energy loss in the suspension (Chateau et al. 2008;DeGiuli et al. 2015). This gives the following estimate for the mean local shear rate The apparent viscosity of the suspending fluid can be estimated via a linearisation of the fluid behaviour at each imposed shear rate, similar to (2.4) and (2.5) for powerlaw and Carreau suspending fluids. We assume that the same relative viscosity as in Newtonian suspensions can be adopted when the suspending fluid is non-Newtonian. Therefore, the shear stress of the suspension can be written as ij D G./ 0 O P ij ;
( It is noteworthy to mention that the above framework is developed for Stokesian suspensions. Also, G./ includes all the information about the microstructure of the suspensions which might not be independent of the shear rate for the case of generalised Newtonian suspending fluids. We therefore investigate these assumptions by means of numerical simulations. One of our goals is to shed light on the recent experimental results in the Stokesian regime by Madraki et al. ( 2017) that show discrepancies with the homogenisation predictions. Moreover, we will show how inertia affects the rheology, G./, and provide a stress closure for non-colloidal suspensions in the presence of inertia. (2017). It is also noteworthy to mention that the variation of the normalised standard deviation for the case of a Newtonian suspending fluid is smaller than that of a shear thinning and larger than that of a shear thickening suspending fluid. This difference is mainly due to the inertial effects (i.e. local particle Reynolds number, see table 1) and it is not related to the particle layering observed in figure 4, as this is similar for all suspending fluids. We show the normalised mean local shear rate for all the simulations in figure 8(a-d). The dashed lines show the prediction of the homogenisation theory given by (3.1). It can be seen that the homogenisation theory gives a good estimate for the mean local shear rate when 6 0:3. However, the theory does not provide accurate predictions of the mean local shear rate as we approach the dense regime, as also observed experimentally in the work of Dagois-Bohy et al. (2015). In particular, it overestimates the values of the mean local shear rate when the solid volume fraction is larger than 0:3. The discrepancy between the mean local shear rate from our simulations and the homogenisation theory is not due to the confinement effects, and consequently, particle layering. Indeed figure 5(a,b) shows a very slight variation of the fluid phase velocity across the gap even when strong particle layering exists, e.g. when D 0:4 for a shear thickening suspending fluid. We can also see that, unlike the bulk solid volume fraction, inertia has a secondary effect on the values of mean local shear rate. Figure 8(d) finally reports all the data in one graph to emphasise that the normalised mean local shear rate is almost independent of the type of suspending fluid.
It is noteworthy to emphasise that a number of experimental works have been carried out showing a general agreement with the homogenisation theory (see e.g. Ovarlez et al. 2006;Chateau et al. 2008;Mahaut et al. 2008;Coussot et al. 2009;Vu et al. 2010;Ovarlez et al. 2012;Dagois-Bohy et al. 2015;Ovarlez et al. 2015;Madraki et al. 2017). Hence, the comparison of our results with the homogenisation theory (see figure 8) is an indirect comparison with the available experimental results. Unlike the previous experimental results, in our simulations we have access to the local velocity and solid volume fraction of the suspension. Therefore, we can provide a better comparison with the theory. Here, we show that F ./ numerical is proportional to ./ homogenisation , such that ./ numerical D ./ homogenisation . Our simulations match the theory with varying from 0.8 (for > 20 %) to 1 (for 6 20 %) almost independently of the suspending fluid. In summary, we show that the general form of the lever function F . / is the same when obtained from homogenisation theory, experiments or simulations.
The effective viscosity of the suspension is defined by the space-and time-averaged wall shear stress divided by the bulk shear rate. Figure 9 shows the non-dimensional effective viscosity for the suspensions considered here (see table 1) versus the particle volume fraction for 0:1 6 Re p 6 10. Note that the non-dimensional effective viscosity is calculated after the transient stage of the simulation and convergence tests have been performed by comparing the statistical results with those obtained using half the number of samples. We also display the non-dimensional effective viscosity predicted by (3.8) in the graphs of figure 9. The results show that, independent of the type of suspending fluid, the non-dimensional effective viscosity increases with the solid volume fraction. In the limit of small Reynolds number, Re p 0:1 (when the inertia is negligible), the results of our simulations follow the predictions by the homogenisation theory. However, as we increase the particle Reynolds number, we deviate from the homogenisation theory, which has been developed for Stokesian suspensions. The increase of the non-dimensional effective viscosity with Re p can be explained as an inertial shear thickening (see e. g. Picano et al. 2013) and will be discussed more in § 3.5.
Here, we provide a more detailed understanding of the role of inertia in the rheology of non-colloidal suspensions. To this end, we first calculate the total suspension stress () as well as the contributions of viscous stress ( v ), Reynolds stress ( R ) and particle where, R includes both particle and fluid Reynolds stresses and it is given as follows:
R D .1 /hv f w f i C ./hv p w p i: (3.12)
Our calculation shows that, for all the simulations performed in this study, the Reynolds stress R is almost zero as the bulk Reynolds number (Re 100Re p ) is in the range of laminar flow. The particle stress includes contributions from the hydrodynamic stresslet, particle acceleration, short-range lubrication correction and collision forces (Lashgari et al. 2014, 2016). It is noteworthy to mention that the total stress is essentially equal to the suspension stress at the wall w as we are not in the dense regime and particle-wall collisions are negligible.
Figures 101112show the stress budget for all of our simulations. The following trend appears immediately clear: independent of the type of suspending fluid, the particle stress mainly contributes to the total stress for > 0:21. This contribution magnifies D. Alghalibi, I. Lashgari, L . Brandt and S. Hormozi To shed light on the role of inertia, we provide scalings for the different components of the stress. First, we show that the viscous stress v indeed scales viscously. To this end, we compute the ratio of the values of the viscous stresses in two successive simulations at constant solid volume fraction, cases .i/ and .i C 1/, where we change the bulk shear rate from P .i/ to P .i C 1/. In the case of viscous scaling, we expect v .i C 1/= v .i/ D TP.i C 1/=P.i/U n , where n D 1 for the Newtonian suspending fluid, n D 0:3 for the shear thinning suspending fluid and n D 1:5 for the shear thickening suspending fluid studied in this work. Figure 14(a-c) shows that our results follow the viscous scaling for all the cases. This also confirms the validity of our computational results.
Second, we investigate the scalings for the particle stress. Unlike the viscous stresses, the scale of particle stresses is not known a priori due to non-zero Re p in our simulations. We expect a scaling of the order of a 2 P 2 when inertia is the main mechanism to transport momentum in a suspension. We therefore divide the values of the particle stresses in two successive simulations at constant solid volume fraction .i/ and .i C 1/, differing in the bulk shear rate P .i/ and P .i C 1/.
Figure 15(a-c) shows that for all types of suspending fluid p .i C 1/= p .i/ is slightly larger than TP.i C 1/=P.i/U indicating that the viscous scaling is still a reasonably good approximation and the inertial contribution to the transport of momentum is a secondary effect. To further confirm this, we report in figure 15 viscous stresses. Assuming an assembly of particles suspended in a viscous fluid that is subject to steady shear under a confining particle pressure P, there exists only one dimensionless control parameter, i.e. the ratio of the microscopic time scale for the particle rearrangement to the macroscopic time scale of the flow P 1 . The microscopic time scale can be determined by balancing the forces applied to a particle, i.e. the imposed pressure P and the drag force F d (Cassar et al. 2005;Andreotti et al. 2013). This implies the existence of the three following microscopic time scales (3:13a c) Here, p and d p are the particle density and particle diameter, whereas the particle terminal velocity and its drag coefficient are denoted by V and C d . When the drag force is negligible and the rearrangement of particles is solely governed by the imposed pressure, the relevant microscopic time scale is t 1 , and consequently, the controlling dimensionless parameter I D t 1 P. When the im osed pressure balances the viscous drag force the microscopic time scale t 2 D d p = P=. f Cd St /, leading to the controlling dimensionless parameter J D t 2 P. Finally, the third scenario, when the imposed pressure balances the turbulent drag force on the particle, is associated with the microscopic time scale t 3 and the controlling dimensionless parameter is J I D t 3 P. Note that the multiplier 4 appears in front of Re p;local since our local particle Reynolds number is based on the particle radius, not the particle diameter. The above closure for the drag coefficient allows us to calculate the ratio of the microscopic time scales: t 1 =t 2 D . p = f / 4 Re p;local =24 and t 1 =t 3 D . p = f / 1=2 .4 Re p;local / 5=16 . Figure 16 shows the phase diagram for the three regimes of non-colloidal suspensions in the plane .t 1 =t 2 ; t 1 =t 3 /. This map is styled after figure 7.13 momentum transport is that associated with particle contacts and collisions, i.e. in this regime particle inertia dominates. The second regime, when t 2 & t 1 and t 2 & t 3 , has t 2 as the dominant microscopic time scale and J f .P local /P =P as the primary controlling dimensionless number. In this regime, the viscous stresses are the main mechanism for momentum transport, the viscous drag regime. The third regime, when t is the dominant microscopic time scale, with relevant the dimensionless number J I f d 2 P 2 =P, is dominated by the Reynolds stresses of the fluid phase, i.e. the turbulent drag at the particle scale (inertial drag regime).
Note that this frictional framework is relevant for dense suspensions in which the particle phase pressure is significant. Moreover, the boundaries between the three regimes discussed above and shown in figure 16 are unknown a priori and finding the exact locations would require an extensive experimental or computational study.
Here, we wish to determine the regime that our simulations for suspensions with 0:4 belong to. It reasonable to assume that a 40 % solid volume fraction is dense enough that the particle phase pressure P is significant, and consequently, the frictional view just introduced can be used to characterise the suspension regimes.
Figure 16 shows that all of the dense simulations considered here (i.e. D 0:4) fall into the regime where the viscous stresses play the main roles in the momentum transport, and hence, the viscous number J is the controlling dimensionless parameter. Therefore, viscous forces are mainly responsible for the dissipation of the power or energy given to the bulk of these suspensions, although Re p is non-zero.
The above discussion raises the question of the role of inertia for the suspensions under investigation here and its influence on the effective viscosity. As shown in § 3.4, increasing Re p results in an enhancement of the suspension effective viscosity. This is attributed to excluded volume effects in the work by Picano et al. (2013): the effective volume fraction of the suspension is higher than the nominal one because of the 'shadow' region (a region with statistically vanishing relative particle flux) around the particles due to inertia in the relative particle motion. Therefore, finite inertia affects the suspension microstructure by increasing the effective , which, in turn, enhances the viscous dissipation. This is why the enhancement of the effective viscosity with inertia is secondary with respect to increases of the nominal solid volume fraction and the suspension viscosity can be scaled back to the Eilers fit, obtained for Stokesian suspensions (see Picano et al. 2013).
Our study shows that the boundary between the inertial I-dominated and the viscous J-dominated regimes, as estimated by scaling arguments with no prefactor, is a good approximation to the real one at least for the suspensions examined here, see figure 16. Despite non-zero values of Re p in our simulations, we are not in the regime where the particle inertia dominates, and therefore, closures for the suspension stresses should D. Alghalibi, I. Lashgari, L . Brandt and S. Hormozi Note that the form of (3.16) is chosen to satisfy the following conditions. First, at small , we expect to be proportional to 2 , since at least two-particle interactions are needed to give rise to the excluded volume effect. Second, the excluded volume must vanish as ! max . In summary, the closure (3.15) provides the suspension shear stress for Newtonian and generalised Newtonian suspending fluids where the microstructure modifications induced by inertia are coded into the G function. We note again here that this closure is valid for suspension in which viscous effects are the still main mechanism for momentum transport, while inertial effects are mainly responsible for changing the microstructure. Note that the local Reynolds number seen by the particles, Re p;local , is the appropriate dimensionless number to model the effects of the microstructure.
We study suspensions of neutrally buoyant spheres in both Newtonian and inelastic non-Newtonian fluids, using interface-resolved direct numerical simulations. The Carreau and power-law models are employed to describe the rheological behaviour of shear thinning and shear thickening carrier fluids, where the fluid viscosity varies instantaneously with the local flow shear rate, P local .x; y; z/. The simulations are performed in a plane Couette configuration. The rheological parameters of the suspending fluids and the particle phase properties are kept constant while we change the volume fraction of the solid phase, 0 6 6 40 %, and the applied shear rate. This allows us to study inertial suspensions with 0 < Re p 6 10. We focus on the bulk properties of the suspension as well as the local behaviour of the particles and the carrier fluid. The main findings of our work are summarised as follows.
The local profile of solid volume fractions show particle layering for nominal volume fractions > 30 %, with the peak of the layers located close to the wall. The distribution of particles across the Couette cell is mainly controlled by geometry and confinement effects, with a weaker dependency on the type of suspending fluid and inertial effects. However, the latter has a strong influence on the velocity profiles of the particle and fluid phases. In the cases with non-Newtonian suspending fluids, viscosity at the particle scale depends on the local shear rate. This results in a larger and smaller local particle Reynolds number as the shear thinning and thickening degree increase, respectively. Inertia also induces a slip velocity between solid and fluid phases with a maximum value close to the wall and zero at the gap centre. This suggests that a suspension flow should be formulated with a two-phase continuum framework close to the walls, where inertial effects are first apparent, although a mixture continuum framework can still provide reasonable predictions of the flow away from the walls, i.e. around the central region. Note that this is also observed for turbulent suspensions and has been successfully used in Costa et al. (2016) to predict turbulent drag in channel flows.
We present the probability density function of the local shear rate, revealing the existence of a wide spectrum of local shear rates; this depends primarily on and secondarily on Re p , which implies the deficiency of theoretical approaches based on mean field values in explaining the mechanics of suspensions. Mean field theories (e.g. Chateau et al. 2008) should be refined to include higher moments, including possibly minima and maxima, of the local shear rate distribution. This requires further investigation and improvement of the available model frameworks.
We demonstrate that the non-dimensional relative viscosity of the suspensions with non-Newtonian carrier fluids can be well predicted by the homogenisation theory of Chateau et al. (2008) in the limit of Re p ! 0, and more accurately for lower . However, adding inertia to the system alters the microstructure and results in a deviation of the relative viscosity of the suspensions from the Stokesian prediction, while the main dissipation mechanism is still viscous. In fact, we show that for the parameter range explored here both particle stresses and fluid stresses are clustered about viscous scalings. We therefore adopt the frictional view of Cassar et al. (2005), Andreotti et al. (2013) to show that the main dimensionless number controlling the mechanics of suspensions is the so-called viscous number, J f .P local /P =P, confirming that viscous stresses are responsible for the momentum transport even when the particle Reynolds number is finite. However, due to inertial effects, the microstructure becomes anisotropic and so-called shadow regions form around particles (these are regions with zero probability of finding another particle, see Picano et al. 2013). This enhances the effective particle volume fraction, and consequently, the viscous dissipation and relative viscosity. We have estimated the volume of these shadow regions from our simulations and included this microstructural effect into a functional form for the relative viscosity. In this way, we provide a prediction for D. Alghalibi, I. Lashgari, L . Brandt and S. Hormozi the added excluded volume due to inertia and a closure for the suspension stress in the case of both Newtonian and generalised Newtonian suspending fluids valid for O.Re p / 10. Note also that once the effective volume fraction is considered, Eilers fit is able to predict the suspension viscosity, confirming that we are still in the viscous regime.
Concerning the frictional framework of Cassar et al. (2005) and Andreotti et al. ( 2013), we note that this view provides scaling for the suspension stresses but no information about the microstructural effects. Therefore, refinements of existing rheological laws would need to include microstructural effects; these local details of the suspension flow can be obtained by well-resolved experimental or computational data, as done here. Our study shows that the so-called inertial shear thickening mode introduced by Picano et al. (2013) belongs to the viscous class of suspensions, yet the inertial microstructure needs to be taken into account to accurately predict the suspension rheology.
This study raises a series of research questions requiring further investigation. First, we show that as we change the type of suspending fluid the contribution of particle stress to the total suspension stress changes. In our system, this is most likely mainly due to the change in the stresslet. A systematic study may therefore be devoted to the investigation of the dependency of the stresslet on the non-Newtonian behaviour of the suspending fluid and its influence on the bulk rheology. Second, the simulations presented here fall into the viscous-dominated suspension regime, although we show how microstructural effects due to inertia are not at all negligible. One may therefore further explore the frictional map by performing simulations/experiments for a wider range of and Re p , one goal being to determine the boundaries between the viscous-and inertia-dominated transport mechanisms. Such an analysis would also suggest how to include microstructural effects in order to refine rheological laws. Third, we consider here the simplest shear flow, with homogeneous bulk shear rate. It is not obvious how inhomogeneous bulk shear rates change the results and how non-local rheology comes into play (see e.g. Pouliquen & Forterre 2009;Bouzid et al. 2013;Henann & Kamrin 2014;Kamrin & Henann 2015). Finally, theoretical approaches based on mean field theory should be improved to include the consequences of the wide spectrum of the local shear rate. To this end, it is essential to study complex features of suspensions with non-Newtonian suspending fluids.
https://doi.org/10.1017/jfm.2018.532 Published online by Cambridge University Press Particle suspensions in shear thinning and shear thickening carrier uids
https://doi.org/10.1017/jfm.2018.532 Published online by Cambridge University Press