Mona Rahmani
email: mrahmani@math.ubc.ca

Abdelkader Hammouti

Anthony Wachs

Momentum balance and stresses in a suspension of spherical particles in a plane Couette flow

Keywords: suspensions, particle-resolved simulations, shear flows, momentum balance

tions ranging between φ = 0.05 and 0.38 and particle Reynolds numbers ranging between Re p = 0.002 and 20, in plane Couette shear flows is investigated using three-dimensional particle resolved numerical simulations. We examine the effects of volume fraction and particle Reynolds number on the macroscopic and microscopic stresses in the fluid phase. The effective viscosity of the suspension is in a good agreement with the previous empirical and experimental studies. At Re p = 20, however, the effective viscosity increases significantly compared to the lower particle Reynolds number simulations in the Stokes flow regime. Examining the stresses over the depth of the Couette gap reveals that this increase in wall shear stresses at high particle Reynolds numbers is mainly due to the significantly higher particle phase stress contributions. Next, we examine the momentum balance in the fluid and particle phase for different regimes to assess the significance of particle/particle interaction and fluid and particle inertia. At the highest particle Reynolds number and volume fraction, the particle inertia plays a dominant role in the momentum balance, the fluid inertia is non-negligible, while the short-lived contact forces are negligible compared to these effects. For all other regimes the fluid inertia is negligible, but the particle inertia and contact forces are important in the momentum balance. Reynolds stresses originated from velocity fluctuations do not contribute significantly to the suspension stresses in any of the regimes we have studied, while the reduction in the shear-induced particle rotation can be a reason for higher wall shear stress at Re p = 20. Finally, we study the kinematics of particles, including their velocity fluctuations, rotation and diffusion over the depth of the Couette gap.The particle diffusion coefficients in the cross flow direction exhibit an abrupt increase at Re p = 20.

I. INTRODUCTION

Suspensions of particles in a fluid are ubiquitous in many natural flows such as mud and gravity currents, as well as in numerous industrial flows, e.g. cosmetics, food processing, slurries, and fresh concrete. In such flows, either the rheological type properties for the entire suspension, or the dynamical behaviour of particles are of interest. Investigations of rheological properties have been particularly focused on finding correlations between the effective viscosity and the volume fraction of the suspension for spheres, e.g. the derivation by Einstein 1 , and the correlations by Krieger-Dougherty 2 and Kruif 3 for dilute to intermediate particle concentrations, and the correlations suggested by Boyer et al. [START_REF] Boyer | Unifying suspension and granular rheology[END_REF] for dense particle concentrations. In addition, attempts have been made to theoretically calculate the effective viscosity of suspensions to second order in volume fraction, e.g. the work by Batchelor & Green 5 and Brady et al. [START_REF] Brady | On the bulk viscosity of suspensions[END_REF] . On the aspect of the dynamical behaviour of particles, the diffusion of particles has been attributed to the anisotropy of local particle stresses in a large number of previous studies [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Morris | Curvilinear flows of noncolloidal suspensions: The role of normal stresses[END_REF][START_REF] Lhuillier | Migration of rigid particles in non-Brownian viscous suspensions[END_REF][START_REF] Dbouk | Normal stresses in concentrated non-Brownian suspensions[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF] . The measured dispersion coefficients are highly dependent on the shear rate gradients and particle/particle interactions [START_REF] Leighton | Measurement of shear-induced self-diffusion in concentrated[END_REF][START_REF] Marchioro | Shear-induced particle diffusivities from numerical simulations[END_REF][START_REF] Ovarlez | Local determination of the constitutive law of a dense suspension of noncolloidal particles through magnetic resonance imaging[END_REF] .

It is well known from previous studies that the microscopic structure of the stresses plays a significant role in determining the macroscopic stresses of the suspension. In a pioneering paper, Batchelor [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] established the contribution of micro-scale stresses to the bulk stresses of the suspension, in an integral form over a volume of the suspension. Conventionally, definitions of macroscopic stresses in suspensions involve averaging of the flow for both phases over an intermediate length scale, sometimes referred to as a mesoscale, which is larger than the scale of particles, but smaller than the large scale of the physical domain, e.g. see references [START_REF] Prosperetti | Pressure forces in disperse two-phase flow[END_REF][START_REF] Zhang | Momentum and energy equations for disperse two-phase flows and their closure for dilute suspensions[END_REF][START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid[END_REF] . By applying the Cauchy's stress principle to linear and angular momentum balance equations and volume averaging, Zhang & Prosperetti [START_REF] Zhang | Physics-based analysis of the hydrodynamic stress in a fluid-particle system[END_REF] found the hydrodynamic stresses in the suspension at the mesoscale. Based on the earlier description of stresses in suspensions by Batchelor [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] , suspension mass and momentum balance models have been developed by Nott & Brady [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF] and Morris & Brady [START_REF] Morris | Pressure-driven flow of a suspension: Buoyancy effects[END_REF] , and later revisited by Lhuillier [START_REF] Lhuillier | Migration of rigid particles in non-Brownian viscous suspensions[END_REF] and Nott et al. [START_REF] Nott | The suspension balance model revisited[END_REF] . These models describe the role of hydrodynamic and non-hydrodynamic forces in contribution of solid particles to the stresses of the whole suspension and particle migration fluxes in Stokes flow regimes.

The bulk stresses in suspensions and migration of particles have been measured in sev-eral experimental studies [START_REF] Boyer | Unifying suspension and granular rheology[END_REF][START_REF] Dbouk | Normal stresses in concentrated non-Brownian suspensions[END_REF][START_REF] Ovarlez | Local determination of the constitutive law of a dense suspension of noncolloidal particles through magnetic resonance imaging[END_REF][START_REF] Fall | Shear thickening and migration in granular suspensions[END_REF] and compared to the volume averaged models mentioned above [START_REF] Dbouk | Normal stresses in concentrated non-Brownian suspensions[END_REF][START_REF] Snook | Dynamics of shear-induced migration of spherical particles in oscillatory pipe flow[END_REF] . However, in laboratory experiments it is difficult to have access to the microstructure of the suspension. To gain more insight into the microscale stresses and their link to mesoscale stresses, numerical simulations that resolve the flow features at the particle length scale are better suited tools. However, in particle-resolved simulations (PRS) the ratio of the largest length scale that can be modelled to the size of particles is quite limited, placing practical constraints on mesoscales of interest. Some previous PRS, although not very numerous, have been dedicated to examining the properties of a suspension of neutrally buoyant rigid spherical particles in a plane Couette flow. Yeo & Maxey 23 used a forced-coupling method to simulate the suspension of spherical particles in a Couette flow for particle volume fractions 0.2 < φ < 0.4 and a particle Reynolds number of Re p = 2. They examined the effects of domain height, which varied between 5 to 15 times the particle diameter, on suspension microstructures and identified wall, buffer and core regions, with the height of each region depending on the volume fraction and total height of the numerical domain. The effects of domain confinement on layering of particles in a plane Couette flow can lead to a significant drop in the effective viscosity of the suspension. This was demonstrated in the simulations of Fornari et al. [START_REF] Fornari | Rheology of confined non-Brownian suspensions[END_REF] , using an immersed boundary method and a maximum height of the Couette gap of 6 times the particle diameter. Using the same method, and in a domain 5 times the particle diameter, Picano et al. [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF] , examined the effects of volume fraction and particle inertia on the effective viscosity of a suspension in a Couette flow. Haddadi & Morris [START_REF] Haddadi | Microstructure and rheology of finite inertia neutrally buoyant suspensions[END_REF] , following an earlier study from the same group 27 , used the lattice-Boltzmann method to examine the microstructure of the suspension in a Couette flow for 0.1 < φ < 0.35 and 0.005 < Re p < 5, in a gap 10 times larger than the particle diameter. They examined the effects of volume fraction and particle Reynolds number on different parts of the particle stresses: the stresslet originated from the surface traction, particle acceleration and particle Reynolds stresses.

Particle-resolved numerical simulations have also been used to study the suspension of spherical particles in a channel flow configuration. Lashgari et al. [START_REF] Lashgari | Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions[END_REF] and Lashgari et al. [START_REF] Lashgari | Channel flow of rigid sphere suspensions: particle dynamics in the inertial regime[END_REF] identified different flow regimes: "laminarlike", "inertial shear-thickening" and "turbulentlike" depending on the particle volume fraction and bulk Reynolds number in a channel flow. Their classification was based on the ratio of the contribution of viscous, particle and Reynolds stresses to the total suspension shear stress. In both studies an immersed bound-and 0.3 the ratio of the channel width to particle diameter was 10.

While volume-averaged suspension balance models [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF] and most numerical simulations of suspensions [START_REF] Marchioro | Shear-induced particle diffusivities from numerical simulations[END_REF][START_REF] Brady | Computer simulation of viscous suspensions[END_REF][START_REF] Drazer | Deterministic and sotchastic behaviour of non-Brownian spheres in sheared suspensions[END_REF][START_REF] Abbas | Dynamics of bidisperse suspensions under Stokes flows: Linear shear flow and sedimentation[END_REF] are focused on the Stokes flow, more recent particle-resolved numerical simulations have explored the inertial effects on stresses in suspensions and the dynamics of the particles. The inertia of particles have been shown to increase the effective viscosity of the suspension [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF][START_REF] Kulkarni | Suspension properties at finite Reynolds number from simulated shear flow[END_REF][START_REF] Lashgari | Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions[END_REF][START_REF] Lashgari | Channel flow of rigid sphere suspensions: particle dynamics in the inertial regime[END_REF][START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF] . Picano et al. [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF] attributed the shear-thickening behaviour of suspensions at higher Reynolds numbers to the creation of shadow regions in the suspension at high fluid inertia that contributes to higher effective volume fractions. For channel flows, Lashgari et al. [START_REF] Lashgari | Channel flow of rigid sphere suspensions: particle dynamics in the inertial regime[END_REF] show that the effective viscosity of the suspension increases at higher Reynolds numbers due to increased Reynolds stresses. They argue that this increase in the effective viscosity at higher Reynolds number is similar to the higher effective viscosity at higher particle fractions as they both can be characterized as the Bagnoldian inertial regime [START_REF] Bagnold | Experiments on a gravity-free dispersion of large solid spheres in a Newtonian fluid under shear[END_REF] . In the latter flow, however, the reason for higher shear stresses on the walls was higher contributions from particle stresses. For a Couette flow configuration, Haddai & Morris [START_REF] Haddadi | Microstructure and rheology of finite inertia neutrally buoyant suspensions[END_REF] showed that increasing the particle Reynolds number increases the particle normal stresses due to high contributions from the Reynolds stresses , while the particle acceleration stresses remain negligible. However, in their study particle and wall shear stresses were not significantly affected by the increase in Re p . Based on laboratory measurements in a coaxial-cylinder rheometer, Linares et al. [START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF] discussed that at sufficiently large Reynolds numbers a transition to turbulence leads to an increase in the effective viscosity. Even before this transition, they observed an increase in the effective viscosity with increasing Reynolds number for volume fractions φ ≤ 0.3.

Our goal in the present study is to use PRS to study the suspension of neutrally buoyant spherical particles in a plane Couette flow, with a gap height 10 times the particle diameter, for volume fractions 0.05 < φ < 0.38 and 0.002 < Re p < 20. We examine the effect of increasing φ and Re p on the bulk suspension shear stress and the microscopic stresses in the fluid phase. This part of our study reveals that at Re p = 20 the shear stress on the walls increase significantly, mainly due to the increase in particle phase shear stress, compared to lower Reynolds number suspensions. To gain more insight into the significance of inertial and particle/particle interaction effects, in different regimes and particularly at high particle Reynolds numbers, we focus on the volume-averaged momentum balance for each phase. We also investigate the potential role Reynolds stresses and particle rotation can play in increasing the suspension stresses. Finally, we briefly examine the diffusion of the particles in the cross flow of the Couette gap for different flow regimes.

Our particle resolved simulations are performed using a distributed Lagrange multiplier/ fictitious domain (DLM/FD) method, e. g. see reference [START_REF] Wachs | PeliGRIFF, a parallel DEM-DL/FD direct numerical simulation tool for 3d particulate flows[END_REF] , and using 24 points to resolve the diameter of each particle, see reference [START_REF] Wachs | Accuracy of Finite Volume/Staggered Grid Distributed Lagrange Multiplier/Fictitious Domain simulations of particulate flows[END_REF] for a sensitivity study of PRS to the number of points per particle diameter. This method has previously been used by Gallier [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF] to study the suspension of spherical particles, however for a fixed configuration of particles.

While most previous simulations of suspension flows have used either the lattice-Boltzmann method [START_REF] Haddadi | Microstructure and rheology of finite inertia neutrally buoyant suspensions[END_REF][START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. part 2. Numerical results[END_REF][START_REF] Ladd | Lattice-Boltzmann simulations of particle-fluid suspensions[END_REF][START_REF] Shardt | Direct simulations of dense suspensions of non-spherical particles[END_REF] , the immersed boundary method [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF][START_REF] Lashgari | Channel flow of rigid sphere suspensions: particle dynamics in the inertial regime[END_REF] , or a forced-coupling method [START_REF] Yeo | Dynamics of concentrated suspensions of non-colloidal particles in Couette flow[END_REF] , the FD method, similarly, has the advantages of using a fixed mesh and simple force calculations for the purpose of simulations of suspension flows [START_REF] Glowinski | A distributed Lagrange multiplier/fictitious domain method for particulate flow[END_REF] .

The paper is organized as follows. In section 2 we describe the physical problem we will study. The numerical method is briefly reviewed in section 3. The macroscopic and microscopic stresses are presented in section 4, and the momentum balance in section 5.

The kinematics of the particles, i.e. their velocity fluctuations, rotation and diffusion are discussed in section 6, with the conclusions stated in section 7.

II. DESPRIPTION OF THE SUSPENSION PROBLEM

We consider a suspension of neutrally buoyant particles in a Newtonian plane Couette flow in a Cartesian system of coordinates (x, y, z), with x, y and z denoting the streamwise, vertical and spanwise directions, respectively; see figure 1. The shear flow is generated by the horizontal movement of planes, vertically at a distance H apart, in opposite directions at velocity U . This gives a bulk shear rate of γ = 2U/H. Spherical particles, with diameter D, are initially randomly positioned in this flow so that the volume fraction of solid particles to the fluid is φ. The particle Reynolds number is defined as Re p = ρ γD 2 /µ, with ρ being the fluid density (and here also the particle density), and µ the fluid viscosity. For the fluid phase, the equations of balance of mass and momentum are stated as [START_REF] Nott | The suspension balance model revisited[END_REF][START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] ∇.u = 0, and where u is the velocity field, σ f the stress tensor for the fluid phase, and b the forces due to the presence of solid particles or the non-hydrodynamic forces per unit volume, with D/Dt denoting the material derivative. In the fluid phase, the stress tensor is

ρ Du Dt = ∇.σ f + b, (1)
σ f = -P I + µ(∇u + ∇u T ), (3)
with I being the identity matrix. The rate of the strain tensor is e = (∇u + ∇u T )/2. The motion of each particle is described by the motion of a solid body in a fixed cartesian frame of reference:

m p du p dt = F + F c , (4)
J dω p dt + ω p × J ω p = T + T c , (5)
with the fluid velocity at the boundary of the particle, u b , found from the rigid body motion constraint:

u b = u p + ω p × r, (6)
where u p and ω p denote the translational and rotational velocity of the centre of mass of the particle, respectively, m p is the mass of each particle, J the inertia tensor, and r is the distance from the centre of mass of particle. It should be noted that ω p × J ω p = 0 for spheres. The hydrodynamic force and torque are denoted by F , and T , respectively, while the contact forces and torques are denoted by F c , and T c , respectively. The hydrodynamic force and torque are calculated as: F = S σ f .ndS, and T = S r × σ.ndS, with n being the vector normal to the boundary of body, and S the boundary enclosing the body. The contact forces between particles (and consequently contact torques) are modelled by by a soft-sphere contact model that includes a normal Hookean elastic restoring force, a normal viscous dissipative force, and a tangential friction force, see Wachs [START_REF] Wachs | GRAIN3D, a flexible DEM approach for particles of arbitrary convex shape -Part I: Numerical model and validations[END_REF] for more details. The description of each of these forces is as follows.

FIG. 2. Illustration of a soft-sphere contact between two spherical particles i and j: G i and G j are the gravity centers of the particles, δ is the overlap between the two particles, k n the stiffness of the equivalent spring, γ n the damping coefficient and n ij and t ij are the normal and tangential unit vectors at the point of contact, respectively.

The normal Hookean elastic restoring force, F el , is defined as

F el = k n δn ij , (7)
where k n is the stiffness of an equivalent spring, δ is the overlap distance between particles i and j (note that δ is negative during contact) and n ij is the unit normal vector along the centerline of the two particles, as illustrated in Figure 2. The viscous dissipative force, F dn , is

F dn = -2γ n M ij U n , (8)
where γ n is the damping coefficient, M ij = (m i m j )/(m i +m j ) is the reduced mass of particles i and j, with m i and m j being the mass of each particle, and U n denotes the normal relative velocity between particles i and j. The tangential friction force, F t , is found from

F t = -min{µ c |F n + F dn |, |F dt |}t ij , (9)
where F n is the total normal force and F dt is the dissipative friction force defined as

F dt = -2γ t M ij U t -k s t ts U t dt, (10)
with µ c being the Coulomb dynamic friction coefficient, γ t the tangential dissipative friction coefficient, U t the tangential relative velocity between particles i, k s the static tangential friction coefficient, t s the time at which contact started and t ij the tangential unit vector along the tangential relative velocity. The term -k s t ts U t dt matters in static or quasistatic particle configurations only and thereby k s has been set to 0 from now on [START_REF] Wachs | GRAIN3D, a flexible DEM approach for particles of arbitrary convex shape -Part I: Numerical model and validations[END_REF] . The total contact force exerted on particle i is the sum of all the contact forces described above originating from contacts with the neighbouring particles j:

F c i = j F c ij = j (F el + F dn + F t ) ij . (11)
In this contact model, we set the stiffness k n such that the maximum ratio of the overlap distance to the radius, 2δ max /D, is less than 0.1% 42 . We also set the normal restitution coefficient e n and the Coulomb tangential friction coefficient µ c to 0.9 and 0.4, respectively.

Based on the value of e n , the damping coefficient, γ n , is computed, while γ t is set to a value with the same order of magnitude as k n , see Wachs [START_REF] Wachs | GRAIN3D, a flexible DEM approach for particles of arbitrary convex shape -Part I: Numerical model and validations[END_REF] for details. We found the effect of the friction coefficient or the normal restitution coefficient to be negligible for volume fractions less than about 40%, in line with Gallier et al. [START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF] who reported that for φ = 0.38 the increase in the Coulomb friction factor has a small increasing effect only on the effective viscosity of the suspension. Moreover, based on a test that compared simulations with and without lubrication forces, we found the effect of lubrication forces to be small for our range of volume fractions. Therefore, our simulations compute lubrication forces up to the limit of the grid size, and no explicit lubrication correction (e.g. see Nguyen & Ladd [START_REF] Nguyen | Lubrication corrections for lattice-Boltzmann simulations of particle suspensions[END_REF] ) is added to the force balance to capture the unresolved sub-grid effects.

III. NUMERICAL METHODS

To solve the equations of motion of the fluid and particle phase, we use a finite volume version of PeliGRIFF, a numerical simulation tool that makes use of distributed Lagrange multiplier and fictitious domain techniques, as described in detail by Wachs [START_REF] Wachs | PeliGRIFF, a parallel DEM-DL/FD direct numerical simulation tool for 3d particulate flows[END_REF] . The numerical aspects and validations of this code have been discussed in many other publications [START_REF] Wachs | PeliGRIFF, a parallel DEM-DL/FD direct numerical simulation tool for 3d particulate flows[END_REF][START_REF] Wachs | Accuracy of Finite Volume/Staggered Grid Distributed Lagrange Multiplier/Fictitious Domain simulations of particulate flows[END_REF][START_REF] Yu | A fictitious domain method for particulate flow with tea transfer[END_REF][START_REF] Yu | A fictitious domain method for dynamic simulation of particle sedimentation in Bingham fluids[END_REF][START_REF] Wachs | A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions[END_REF] .

To capture the Stokes behaviour of the suspensions, as well as the transition to inertial regimes, we perform a series of simulations at Re p = 0.002, 0.02, 0.2, 2, and 20. For each of these particle Reynolds numbers, volume fractions of φ = 0.05, 0.1, 0.2, 0.3 and 0.38 are considered. These simulations are performed in a box with

L x = L y = L z = H = 10D.
According to the simulations of Yeo & Maxey 23 and Fornari et al. [START_REF] Fornari | Rheology of confined non-Brownian suspensions[END_REF] this domain height is sufficient for particle layering near walls not to influence the bulk stresses. To resolve the diameter of each particle, 24 points are used, with the total domain mesh size being 240 3 .

For high Reynolds numbers and volume fractions, we have tested the sensitivity of the results to the mesh size. For Re p = 20 and φ = 0.38, using 32 points per diameter, the average shear stress on the walls was within 7% of its value using 24 points. turbulence [START_REF] Daviaud | Subcritical transition to turbulence in plane couette flow[END_REF][START_REF] Bottin | Discontinuous transition to spatiotemporal intermittency in plane couette flow[END_REF][START_REF] Philip | From temporal to spatiotemporal dynamics in transitional plane couette flow[END_REF][START_REF] Duguet | Oblique laminar-turbulent interfaces in plane shear flows[END_REF][START_REF] Couliou | Growth dynamics of turbulent spots in plane couette flow[END_REF] . Therefore, we do not expect any large-scale velocity fluctuations to develop due to hydrodynamic instabilities. Based on experimental measurements in pipe flows, Matas et al. [START_REF] Matas | Transition to turbulence in particulate pipe flow[END_REF] showed that due to the presence of particles, the critical Reynolds number for transitionto turbulence can either increase or decrease depending on the ratio of pipe diameter to particle size and the volume fraction of solid particles. However, here we did not find any particle induced transitional behaviour as our large-scale Reynolds numbers are small. The range of large-scale Reynolds numbers we consider here are well below the critical Reynolds numbers for the transition to turbulence found by Linares et al. [START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF]

in coaxial

Couette flow experiments.

The dimensional time step of the fluid solver in all simulations is less than 0.0005 γ-1 .

This time step is small enough to ensure the inertial and viscous time scales of the flow are well resolved. The time step for solving the motion of particles is at least 15 times smaller than the physical contact time of the spheres.

IV. STRESSES IN THE SUSPENSION A. macroscopic shear stresses

One of the main objectives of studying suspension flows is to characterize the macroscopic and microscopic stresses in the flow to describe the rheological behaviour of the suspensions.

A measure of the macroscopic shear stresses of the suspension is the effective viscosity obtained by non-dimensionalizing the mean shear stress on the walls by µ γ, the viscous shear stress on walls without the particle effects. The time evolution of the effective viscosity of the suspension for Re p = 20 and different volume fractions, is presented in figure 3. This figure shows that in our simulations the wall shear stress starts from zero, in an initial state in which the fluid and particle velocities are set to 2(z -H/2)U , and evolves to a statistically steady state, where the micro-structure features of the flow are fully developed but still evolving in time. In this steady state, the mean shear stresses fluctuate about a mean value as the interaction of particles through collisions and hydrodynamic forces continuously modifies the microstructure of the suspension. These fluctuations are much more pronounced at higher volume fractions due to the stronger interaction of particles.

The dimensionless time required to reach a steady state is less than 100, as reported by Drazer et al. [START_REF] Drazer | Microstructure and velocity fluctuations in sheared suspensions[END_REF] . The horizontally averaged streamwise fluid velocity and volume fraction, also time averaged over the statistically steady state period, are presented in figure 4. The averaged velocity profiles are close to the initial shear profile for all cases, with a small deviation observed at φ = 0.38. The averaged volume fraction profile however exhibits strong fluctuations, especially at φ = 0.38. These fluctuations are caused by the random initial positioning of the particles and the fact that their net deviation from these initial positions is small over time. At φ = 0.38, the concentration of particles close to the walls is on average higher than that at the center of the Couette gap, which is an indication of layering of particle close to the solid walls. The particle resolved simulations of Gallier et al. [START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF] revealed a similar profile of the local concentration of particles at φ = 0.4.

The dimensionless shear stress on the walls (or the effective viscosity), time-averaged in the statistically steady state, as a function of the volume fraction is shown in figure 5 for different particle Reynolds numbers. Our numerical effective viscosity are compared to the correlations between the effective viscosity and the volume fraction proposed by e.g. show a good agreement with the empirical correlations of Eilers, Krieger-Dougherty, and number of Re p = 20, the wall shear stress exhibits a significant increase above the assumed empirical correlations. The abrupt increase in τ w xy at Re p = 20 follows a monotonic increase of the wall shear stress with the Reynolds number. Figure 5 suggests that the bulk shear stress goes through a transition as the Reynolds number increases to Re p = 20.

The effect of particle Reynolds number on the effective viscosity of suspensions has been previously examined in Couette flows. Kulkarni & Morris 27 observed a particle Reynolds number effect only at φ = 0.3, the highest volume fraction they simulated. Their computed effective viscosity increased at higher particle Reynolds numbers (Re p = 1 was their maximum particle Reynolds number), while the effect was non-monotonic for lower particle Reynolds numbers. Picano et al. [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF] obtained a monotonic increase in the effective viscosity when increasing Re p from 1 to 5, with this effect being more pronounced at higher particle volume fractions. Our results are in agreement with these finding of Picano et a. [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF] , for a larger range of particle Reynolds numbers. For a suspension of particles in a coaxial Couette flow and over a wide wide range of 3 < Re p < 100, Linares et al. [START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF] measured a monotonic increase of the effective viscosity with particle Reynolds number, but only for φ ≤ 0.3. For higher volume fractions, the effective viscosity revealed either a non-monotonic behaviour or a decreasing trend with increasing Re p . Our simulations with φ > 0.3 clearly do not agree with this finding of Lineras et al. [START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF] , and for φ < 0.3 our simulations show a slower increase in the effective viscosity with increasing Re p compared to their measurement. The former finding is in agreement with the simulations of Picano et al. [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF] The increase in shear stress in suspensions at higher Reynolds numbers has been attributed to the development of shadow region in case of inertial particles in viscous flow regimes by Picano et al. [START_REF] Picano | Shear thickening in non-Brownian suspensions: an excluded volume effect[END_REF] and increase in Reynolds stresses in turbulent flow regimes with large velocity fluctuations by Lashgari et al. [START_REF] Lashgari | Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions[END_REF][START_REF] Lashgari | Channel flow of rigid sphere suspensions: particle dynamics in the inertial regime[END_REF] . Numerical simulations of Mikulencak & Morris [START_REF] Mikulencak | Stationary shear flow around fixed and free bodies at finite Reynolds number[END_REF] for a stationary sphere in a shear flow indicated that inertial effects enhance the contributions of the particle to the total shear stress of the suspension. Patankar & Hu 57 , using numerical simulations of two-dimensional circular particles in a shear flow, pointed out that the reduction in the angular velocity of particles, that occurs evidently for Re p > 1 contributes to higher velocity gradients between particles and wall and therefore larger shear stresses. In the limit of high Reynolds numbers, Re p > 100, and relatively low volume fractions, φ < 0.3, Linares at al. [START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF] discuss that high fluid velocity fluctuations and a transition to turbulence result in higher effective viscosity. We examine the importance of particle and fluid inertia in momentum balance at higher Reynolds numbers in section V, and the significance of Reynolds stresses and reduction in particle rotation at higher Reynolds numbers in section VI.

B. microscopic shear and normal stresses

It is established in the literature that the macroscopic behaviour of the suspension is highly dependent on its microscopic features, e.g. see [START_REF] Zhang | Physics-based analysis of the hydrodynamic stress in a fluid-particle system[END_REF][START_REF] Haddadi | Microstructure and rheology of finite inertia neutrally buoyant suspensions[END_REF][START_REF] Esteghamatian | Micro/meso simulation of a fluidized bed in a homogeneous bubbling regime[END_REF] among other studies. At the microscopic level, to distinguish between fluid and particle contributions to the stresses, conventionally the whole suspension stresses are split into fluid and solid parts as

7 σ = (1 -χ)σ f + χσ p , ( 12
)
where χ is the phase indicator, with χ = 0 in the fluid phase and χ = 1 in the solid phase.

The fluid phase stress tensor is σ f = -P I + 2µe, and σ p denotes the stress tensor in the particle phase.

The microstructure of the fluid phase stresses, for an instantaneous snapshot, is shown in figure 6 for the shear stress, τ f xy , and in figure 7 for the deviatoric part of the normal stress, σ f yy . These figures show that the local shear and normal stresses can be significantly larger than the bulk shear stress, µ γ. Particularly in the gaps between the particles and particles and walls, where the hydrodynamic interactions between the particles and particles and walls strain and pressurize the fluid and the stresses exhibit high local values. These high local shear and normal stresses are larger in magnitude and size of regions of occurrence when the volume fraction increases from φ = 0.05 to φ = 0.38. By increasing the particle Reynolds number from Re p = 0.02 to Re p = 20, the regions of high concentration of stresses increase slightly in area in the fluid phase. However, this change is much less significant compared to the effects of increasing the volume fraction of particles.

In figure 6, in regions of the fluid phase between the particles, where the distance between two particles is of order of one particle size, the fluid experiences zero or smaller than µ γ shear stresses. This effect is particularly pronounced in the wakes behind the particles. Local τ f xy is marginally higher than the bulk background shear stress. Near the walls, where the solid volume fraction goes to zero, the fluid shear stress abruptly increases to values close to the computed time averaged wall shear stress. This rapid change in τ f xy occurs over a one particle diameter length scale and is similar to the sharp increase in the viscous stresses near the wall in simulations of Lashgari et al. [START_REF] Lashgari | Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions[END_REF] higher than the average shear stress from the fluid phase contribution in the middle of the domain. The difference is the contribution of the solid phase, τ p xy , that includes the stresslet, particle acceleration stresses and particle Reynolds stresses [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Haddadi | Microstructure and rheology of finite inertia neutrally buoyant suspensions[END_REF] . Similar to this finding, Shakib-Manesh et al. [START_REF] Shakib-Manesh | Shear stress in a Couette flow of liquid-particle suspensions[END_REF] showed that in their PRS of circular particles in a Couette flow, the increased shear stress at higher Reynolds numbers was due to enhanced particle shear stress.

Figure 9 delineates that τ p xy (the difference between the shear stress on the walls and the fluid shear stress in the middle of the domain shown by the vertical dashed line) increases with increasing φ and Re p . While at higher φ the stronger particle/particle and particle/wall interactions lead to larger contributions from the stresslet to the shear stress 7,61 , at higher Re p it is expected that the acceleration and Reynolds stresses also contribute significantly to the particle shear stress [START_REF] Lashgari | Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions[END_REF][START_REF] Shakib-Manesh | Shear stress in a Couette flow of liquid-particle suspensions[END_REF] . To gain more insight into the significance of particle and fluid inertia in comparison to particle/particle interactions, we will examine the momentum balance for particle and fluid phase in section 5. The other potential reasons for higher particle shear stress at higher Reynolds numbers, larger Reynolds stresses and reduced particle rotation, will be investigated in section 6.

V. MOMENTUM BALANCE ANALYSIS

The momentum balance of the suspension is described by Cauchy's equation of motion

∇.σ + b = ρ Du Dt . (13)
A continuum volume-averaged model of equation 13 for Stokes regimes is given by the suspension balance model (SBM) [START_REF] Nott | Pressure-driven flow of suspensions: simulation and theory[END_REF][START_REF] Nott | The suspension balance model revisited[END_REF] . SBM describes closures for the balance of mass and momentum for solid and fluid phases, and also the whole suspension. It should also be noted that SBM's continuity ensemble equation, stated as ∇.u = 0, is inherent in our fictitious domain formulation, e.g. see [START_REF] Wachs | PeliGRIFF, a parallel DEM-DL/FD direct numerical simulation tool for 3d particulate flows[END_REF] , which makes the discussion of mass conservation trivial here.

By volume averaging equation 13, in this section we focus on the significance of each term in the force balance in different regimes.

By multiplying equation 13 by the phase indicator χ and averaging over the entire domain, the momentum balance for the particle phase is obtained 7

∇. χσ p + n p F p + b p = n p ρ Du Dt p . (14)
In the equation above, the first term on the right hand side is the divergence of volume averaged stresses in the particle phase, n p F p is the volume averaged sum of hydrodynamic particle forces, and b p represents the contributions from particle contact forces per volume. The term on the right hand side is the volume averaged sum of particles acceleration.

Similarly, the fluid phase momentum balance can be written as

∇. (1 -χ)σ f -n p F p + b f = ρ Du Dt f , (15)
where the first term is the divergence of the volume averaged stresses in the fluid phase, b f is the contribution of the particle contacts to the forces in the fluid phase, and ρ Du Dt f is the inertial force due to the acceleration of the fluid phase. The volume averaged sum of hydrodynamic forces on particles represents the momentum exchange between the particle and fluid phase.

First, we consider the instantaneous fluid phase momentum balance as stated in equation (15). In the Stokes flow regime, the fluid phase inertia can be neglected. In our particle resolved simulations the forces in the fluid phase due to particle contacts are already captured as part of the divergence of the fluid phase stresses. Moreover, an exact evaluation of this term is not feasible. Therefore, we did not include b f explicitly. This reduces equation 15, in the limit of low Re p , to simply

∇. (1 -χ)σ f = n p F p . ( 16
)
This equation states that the average of divergence of the stresses in the fluid phase is equal to the average of the interphase forces on the particles surface, and has also been previously derived in a similar form by Lhuilier 9 and Dbouk et al. [START_REF] Dbouk | Shear-induced particle migration: Predictions from experimental evaluation of the particle stress tensor[END_REF] . Upon applying the divergence theorem, making use of the periodicity of the boundaries in x and z directions and assuming that at the top and bottom walls χ = 0 everywhere (particle contact area with the wall is negligible), equation 16 simply gives the following force balance in three directions

(n p F x p , n p F y p , n p F z p ) = 1 H (∆τ w xy , ∆σ w yy , ∆τ w zy ). (17)
In the equation above ∆ denotes the difference between stresses on the top and bottom walls.

The force balance in equation ( 17) is examined in figure 10 Figure 10 shows that at φ = 0.05 and Re p = 0.02 the net particle hydrodynamic forces are very small compared to the shear force as the flow around each particle is close to being perfectly symmetric at this low Reynolds number and particle concentration. At φ = 0.38, however, the net particle hydrodynamic forces fluctuate with larger amplitudes, particularly in the y direction due to the interactions between particles and particles and walls. Intermittent spikes in the forces occur due to the particle/particle and particle/wall contacts. At φ = 0.05 and Re p = 20, the sum of particle hydrodynamic forces exhibits high frequency fluctuations while its amplitude remains relatively small. For all these three cases, we get a close match between the instantaneous n F p and ∇. Now, we turn our attention to the particle phase momentum balance, as stated in equation (14). Since χ = 0 on the top and bottom wall, the application of the divergence theorem to equation ( 14) makes the term ∇. χσ p vanish. Therefore, the momentum balance is between the other three terms in equation ( 14). Since computing the exact contributions of the particle/particle and particle/wall contacts to the stresses of the particle phase, i.e. b p , is not straightforward, we examine the balance between terms n p F p and n p ρ Du Dt In the streamwise (x) direction, there is a close match between the PDFs of n p F p and n p ρ Du Dt p for the Re p = 20 and φ = 0.38 case. This indicates that for this high particle volume fraction and Reynolds number the sum of the streamwise hydrodynamic forces on the particles captures the entire forces on particles and is balanced by the sum of the particle acceleration forces. Therefore, the contributions from the contacts are already accounted for in the hydrodynamic force part. This is because the enhanced particle/particle interactions lead to short-lived contact forces acting on particles that directly transfer into higher instantaneous particle accelerations. So the term b p is negligible in the particle momentum

Re p = 0.02 ∇. (1 -χ)σ f -n p F p = 0 (1) (1)
n p F p + b p = n p ρ Du Dt p (2) (2)
Re p = 20 (1) ∇. (1 -χ)σ f -n p F p = ρ Du Dt f (3)
(2)

n p F p = n p ρ Du Dt p (4)
TABLE I. Summary of approximations for momentum balance equations in the fluid and particle phase for different flow regimes. Equation ( 4) is valid only in the streamwise and spanwise directions, in the vertical direction equation ( 2) should be used instead.

the wall effects on particles, the variation in PDFs of n p F p increase dramatically, while To summarize the findings in figures 11, 12 and 13, we present the approximate fluid and particle momentum balance equations in table I. This table shows that in both low volume fraction and low Reynolds number regimes the fluid inertia can be neglected in the balance, while the inter-particle contact forces are non-negligible, even at φ = 0.05. At the high Reynolds number and volume fraction regime both fluid and particle inertia terms are important and their sum balances the gradients in the fluid stresses (except in the vertical direction). The inter-particle contact forces transfer into particle inertia rapidly and can be neglected as a distinct term (except in the vertical direction).

VI. KINEMATICS OF PARTICLES A. Particle velocity fluctuations

We examine the magnitude of time averaged particle velocity fluctuations in all three directions, x, y, and z, as a function of y in figure 14. The magnitude of particle velocity fluctuation in x, y and z directions are the root mean square of the instantaneous velocity fluctuations, defined, respectively, as

u p = u p -2 z p - H 2 U 2 t 1/2 , v p = (v p ) 2 t 1/2 , w p = (w p ) 2 t 1/2 , ( 18
)
where t denotes temporal averaging for each particle. In equations ( 18), each component of particle velocity fluctuation is defined based on the difference between the particle velocity and the local undisturbed fluid velocity. These velocity fluctuation measures bear similarity with the granular temperature, T , defined as T = (u pu p t ).(u pu p t ) t /3, e.g. see Wylie et al. [START_REF] Wylie | Rheology of suspensions with high particle inertia and moderate fluid inertia[END_REF] , which is a measure of the kinetic energy of particles due to their velocity fluctuations. However, here we are more interested in the magnitude of each component of velocity fluctuations rather than the rms of all components as defined in the granular temperature.

Figure 14 shows that the amplitude of particle velocity fluctuations in the streamwise (shear) direction is larger than those in the other two directions. Velocity fluctuations in the vertical direction also exhibit larger amplitudes compared to the spanwise direction. These indicate that the particle velocity fluctuations are not isotropic in three directions as shearing and wall effects favour velocity fluctuations. As φ increases and the particle interactions are enhanced, the particle velocity fluctuations increase in amplitude in all three directions. Additionally, by increasing φ, the vertical direction dependence of u p significantly changes and its maximum shifts toward the center at φ = 0.38. The reason is that at higher φ the role of interaction of particles with each other in inducing streamwise direction particle velocity fluctuations becomes comparable to the effects of the interaction of particles with walls. By increasing Re p , u p , v p and w p become more uniform over the depth of the domain as the effects of inertia in inducing velocity fluctuations become comparable to particle/particle and particle/wall interaction effects. However, the dependence of the amplitude of the velocity fluctuations on Re p is not monotonic and much weaker than the dependence on φ. The velocity fluctuations v p and w p are less than 0.02U over the entire depth of the gap. Therefore, the contribution of Reynolds stresses to the bulk shear stress is negligible in all our simulations. These velocity fluctuation amplitudes fall into the "laminarlike" (low volume fraction and low inertia) or the "inertial shear-thickening" (high volume fraction and presence of inertia at particle scale) regimes of velocity fluctuations classified by Lashgari et al. [START_REF] Lashgari | Laminar, turbulent, and inertial shear-thickening regimes in channel flow of neutrally buoyant particle suspensions[END_REF] in a channel flow. They also concluded that in these regimes the contribution of Reynolds stresses to total wall shear stress is small. While particle accelerations are significantly higher at Re p = 20 compared to those at Re p = 0.02 (see figures 11 to 13 ), time averaged particle velocity fluctuations are not necessarily higher at Re p = 20. This is due to the highly intermittent nature of particle accelerations and forces at high particle Reynolds numbers (see figure 10 ) that can generate high velocity fluctuations over a short period of time, but the time averaged velocity fluctuations are smoothed out to values comparable to those at low Reynolds numbers. However, since the Stokes number of particles are small in all our simulations (with the Stokes number being St = Re p /18) we do not expect any transition to highly collisional regimes, in a sense of "agitated regimes" as described by Tsao & Koch 64 , Sangani et al. [START_REF] Sangani | Simple shear flows of dense gas-solid suspensions at finite stokes numbers[END_REF] and Parmentier & Simonin 66 , or a significant increase in velocity fluctuations by increasing the Reynolds number as found by Abbas et al. [START_REF] Abbas | Flow of particles suspended in a sheared viscous fluid: Effects of finite inertia and inelastic collisions[END_REF] .

It is noteworthy that we found the fluid phase velocity fluctuations to have a very similar behaviour to the particle velocity fluctuations. We did not find any transition to turbulence, as described by Linares et al. [START_REF] Linares-Guerrero | Effects of inertia and turbulence on rheological measurements of neutrally buoyant suspensions[END_REF] in their laboratory experiments, in any of our simulations.

Therefore, at high Reynolds numbers, while our particles become inertial through experiencing high accelerations, the fluid remains laminar.

B. Particle rotation

The rotational kinetic energy of the particles serves as a measure of the interaction of particles with each other and with the fluid phase. The reduction in particle rotation at higher Reynolds numbers in shear direction, Ω z , has been linked to the increase in the shear stress of the suspension [START_REF] Patankar | Finite Reynolds number effect on the rheology of a dilute suspension of neutrally buoyant circular particles in a Newtonian fluid[END_REF] . Figure 15 shows the rotational velocities in all directions (root of time-averaged squares of rotational velocities). This figure shows that Ω x , and Ω y , rotational velocities caused by the particle collisions and particle hydrodynamic interactions, increase by increasing φ, as was also observed by Drazer et al. [START_REF] Drazer | Microstructure and velocity fluctuations in sheared suspensions[END_REF] . Increasing Re p from 0.002 to 2 had a slight increasing effect on the average Ω x and Ω y (data not shown here). At

Re p = 20, however, the average particle rotations in x and y directions are suppressed at high volume fractions where the particle interactions are strong. Increasing the Reynolds number results in a more uniform distribution of particle rotational velocities across the depth of the gap. The particle rotation in z direction, mainly caused by the background shear and rolling of particles, is about 0.5 γ, in agreement with previous studies [START_REF] Drazer | Microstructure and velocity fluctuations in sheared suspensions[END_REF][START_REF] Mikulencak | Stationary shear flow around fixed and free bodies at finite Reynolds number[END_REF] . For all particle Reynolds numbers, following a monotonic increase with increasing φ, the average Ω z diminished at φ = 0.38 due to the suppression of free particle rolling at this high volume fraction. Similarly, the increase in Re p has an enhancing effect on Ω z up to Re p = 2 (not all data shown here). At Re p = 20, the shear induced particle rotations are suppressed by Hu's 57 analysis, this reduction Ω z can contribute to higher wall shear stresses at higher Reynolds numbers.

C. Particle diffusion

Shear and particle interaction induced diffusion of particles in suspensions have been studied extensively previously [START_REF] Leighton | Measurement of shear-induced self-diffusion in concentrated[END_REF][START_REF] Marchioro | Shear-induced particle diffusivities from numerical simulations[END_REF][START_REF] Drazer | Deterministic and sotchastic behaviour of non-Brownian spheres in sheared suspensions[END_REF][START_REF] Eckstein | Self-diffusion of particles in shear flow of a suspension[END_REF] and attributed to the asymmetric and anisotropic stresses around the particles. Since in our Couette flow simulations the shear rate is constant, the diffusion of the particles is solely due to the interaction of particles. To make our computation of the diffusion coefficients comparable with these previous studies, we have normalized the dimensional diffusion coefficients by γ(D/2) 2 , i.e. now we are using the particle radius instead of particle diameter.

In general D yy and D zz have an increasing trend with increasing φ. However, at higher volume fractions this trend becomes occasionally non-monotonic, as also observed in previous studies [START_REF] Breedveld | Measuring shear-induced self-diffusion in a counterrotating geometry[END_REF][START_REF] Sierou | Shear-induced self-diffusion in non-colloidal suspensions[END_REF] for averaging in our and other numerical simulations. For most cases, particles are more diffusive in the vertical direction compared to the streamwise direction, in agreement with the previous estimations of the particle diffusion coefficients in both directions.

At Re p = 20, the diffusion coefficients exhibit significantly higher values compared to lower Reynolds numbers and also compared to previous quantifications of diffusion coefficients at zero of finite Reynolds numbers. We attribute this transition at Re p = 20 to high spikes in particle acceleration (see figures 11, 12 and 13) that result in high intermittent velocity fluctuations and consequently large particle displacements in one direction rather than the oscillatory motion of particles at lower Reynolds numbers. It is noteworthy that the maximum transverse displacement of the particles in the duration of our simulations is 4D.

This is less than half height of the domain and so our diffusion coefficients do not reflect the effects of confined domain height on averaged particle diffusion. Simulations should be run for a longer time to gain more insight into confinement effects. However, this investigation is beyond the scope of the present study.

VII. CONCLUSIONS

We have studied the suspension of spherical particles in a plane Couette flow using a Lagrange multiplier/ fictitious domain method, for particle volume fractions 0.05 < φ < 0.38 and particle Reynolds numbers 0.002 < Re p < 20.

At the macroscopic level, our wall shear stress, presented in terms of the effective viscosity, was in a good agreement with previous empirical correlations for the effective viscosity as a function of particle volume fraction. At Re p = 20, the effective viscosity increased significantly compared to the lower particle Reynolds numbers due to an abrupt increase in the particle phase shear stress. At the microscopic level, the local shear and normal stresses in the fluid phase exhibited a much wider variation at higher particle volume fractions owing to the local shearing or compressing of the fluid in the small gaps between the particles or the low shear rate and low pressure regions in the wakes of particles. The increase in the particle Reynolds number caused only a small increase in the variation of local shear and normal stresses in the fluid phase. To gain more insight into the nature of the transitional behaviour at Re p = 20, we examined the instantaneous Cauchy momentum balance equation for each phase for the bulk of the suspension for four different regimes with low to high particle volume fractions and particle Reynolds numbers. Our analysis revealed that in the limit of large φ and Re p the momentum balance is distinctly different form other regimes due to particle inertia and enhanced particle interactions by inertial effects. In this regime, the contact forces are short-lived compared to the inertial effects and they instantaneously contribute to the particle acceleration. Therefore, the momentum balance for the particle phase at Re p = 20 and φ = 0.38 is simplified to the equality of the particle hydrodynamic forces and particle Re p , the fluid inertia was negligible and the balance was simply between the gradient of fluid phase stresses and the particle hydrodynamic forces. In these regimes, the particle inertia term was however important in the particle phase momentum balance. The momentum balance for the particle phase was between the particle hydrodynamic forces, the contact forces and the particle acceleration forces.

At Re p = 20 and for all φ, the sum of hydrodynamic forces on particles, and the sum of particle acceleration forces revealed high intermittency with high amplitudes of fluctuations.

While these spikes in particle acceleration resulted in intermittently high velocity fluctuations, the time-averaged amplitude of velocity fluctuations were not higher than those at lower particle Reynolds numbers. Therefore, similar to lower Reynolds numbers, at Re p = 20 we did not find the effect of Reynolds stresses on the bulk stress to be significant. High amplitude particle accelerations at Re p = 20 however enhanced the particle diffusion in the cross flow direction of the Couette gap. The effect of particle Reynolds number on particle diffusion coefficients in the vertical and spanwise directions was small and non-monotonic for the range 0.002 < Re p < 2. The diffusion coefficients in general increased at higher φ as the particle interactions were stronger. The shear-induced particle rotational velocities diminished significantly at Re p = 20, which can contribute to higher wall shear stress.

Finally, our simulations revealed that the particle diffusion, velocity fluctuations and hydrodynamic forces are anisotropic in the streamwise, spanwise and vertical directions.

At Re p = 20, the particle velocity fluctuations and rotational velocities exhibit a more uniform and isotropic distribution across the Couette gap as the inertial effects become more dominant. However, addressing the isotropy of the suspension in PRS certainly requires larger domain sizes and more computational resources, directions to be pursued in future studies. While we have explored some of the distinct properties of the suspension at Re p = 20, many other aspects of the inertial effects on stresses and particle migration in suspensions at high Reynolds numbers and different flow configurations are yet to be understood.

FIG. 1 .

 1 FIG. 1. Configuration of the suspension problem considered.

3 .

 3 figure 3. Our large-scale Reynolds numbers (defined based on the flow height), are smaller than the critical Reynolds number for the occurrence of secondary flows or transition to

38 FIG. 3 .

 383 FIG. 3. Time evolution of the mean shear stress on the wall, τ w xy , non-dimensionalized by µ γ (which gives the effective viscosity) for different volume fractions: φ = 0.05, 0.1, 0.2, 0.3, and 0.38 at Re p = 20.

  . Our computed values of effective viscosity are in general lower than those obtained in numerical simulations by Kulkarni & Morris 27 and Picano et al. 25 , while they are closer to values computed by Yeo & and Maxey 23 .

FIG. 5 .

 5 FIG. 5. Variation of the effective viscosity of the suspension, time average of the mean shear stress on walls, < τ w xy > t , divided by µ γ, with the volume fraction for different particle Reynolds numbers. The dashed line delineates the empirical correlation of Eilers 55 , the dashed-dotted line the correlation of Krieger-Dougherty 2 , and the dotted line, the correlation of Boyer et al. 4 .

  FIG. 6. Snapshots of the instantaneous micro-structure of the shear stress in the fluid phase, τ f xy /(µ γ), at a time in the statistically steady state. Each snapshots shows a slice at z = L z /2, for Re p = 0.02, panels (a) and (b), and Re p = 20, panels (c) and (d). In each row the left panels (a) and (c) corresponds to a low volume fraction: φ = 0.05, and the right panel (b) and (d) correspond to a higher volume fraction: φ = 0.38. The stresses inside the solid particles are shown as zero.

  in a channel flow with high Reynolds numbers or high particle volume fractions. In a fluid layer very close to the walls the shear stress in the fluid phase represents the entire suspension's shear stress as χ = 0 and τ xy = τ f xy . The wall shear stress in all cases, and especially for high volume fractions, is considerably

FIG. 8 .

 8 FIG. 8. Probability density functions of the local shear and normal stresses in the fluid phase: (a) τ f xy /(µ γ) and (b) σ f yy /(µ γ).

FIG. 9 .

 9 FIG. 9. The fluid phase horizontally averaged shear stress, τ f xy , for Re p = 0.02, panels (a) and (b), and Re p = 20, panels (c) and (d), and φ = 0.05, left panels (a) and (c), and φ = 0.38, right panels (b) and (d). The dashed lines show the average shear stresses within the suspension, i.e. averaged from y = D to y = L z -D. The dots indicate the time averaged shear stress on the walls.

  for the lowest and highest volume fraction simulations, i.e. φ = 0.05 and 0.38, at Reynolds numbers Re p = 2, and 20. In this figure, the instantaneous volume averaged sums of the hydrodynamic forces on all particles, n p F p , are presented in all three directions. The volume average of the divergence of the fluid phase stresses are shown for comparison, however only at a few times due to the noisy nature of the data. All terms have been non-dimensionalized by the shearing force scale: F shear = µ γA, with A being the surface area of the top (or bottom) wall.

( 1 -

 1 FIG. 10. Instantaneous fluid phase momentum balances over a 4 dimensionless time unit span for Re p = 0.02, the top panels (a) and (b), and Re p = 20, the bottom panels (c) and (d). The left panels (a) and (c) show φ = 0.05 and the right panels (b) and (d) show φ = 0.38. The lines indicate the sum of hydrodynamic forces on particles, and the dots the divergence of stresses in the fluid phase, non-dimensionalized by the bulk shear force.

  difference to b p . Figures 11, 12 and 13 compare the balance between these two terms in directions x, y and z, respectively, by examining the probability density function (PDF) of their instantaneous values over a sufficiently large time span in the statistically steady state. As both terms n p F p and n p ρ Du Dt p fluctuate rapidly in time, their instantaneous comparison cannot be interpreted in a meaningful way. On the other hand, comparing their PDFs provides a means to examine their balance statistically.

FIG. 11 .

 11 FIG. 11. The probability density function of the instantaneous volume averaged sum of hydrodynamic forces on particles, n p F p , and particle phase inertia n p ρ Du Dt p

FIG. 13 .

 13 FIG. 13. The probability density function of the instantaneous volume averaged sum of hydrodynamic forces on particles, n p F p , and particle phase inertia n p ρ Du Dt p

  the variation in PDFs of n p ρ Du Dt p increase to a lesser extent. Particularly at Re p = 0.02 and φ = 0.38, the contributions of the contact forces in b p are significant in closing the momentum balance in the y direction. At Re p = 20 and φ = 0.38, also the term b p is necessary to close the balance in the y direction. In the spanwise (z) direction, similar to the other directions the variance in PDFs of n p F p and n p ρ Du Dt p increases by increasing either φ or Re p , while increasing φ and Re p simultaneously has a much more pronounced effect. In this direction a good balance holds between the terms n p F p and n p ρ Du Dt p and therefore, b p is negligible.

FIG. 14 .

 14 FIG. 14. Dependence of the time-averaged magnitude of particle velocity fluctuations, u p , v p and w p , normalized by U , on the time-averaged vertical particle location. The colors show the particle velocity fluctuations in: the streamwise direction (blue), vertical direction (red), and spanwise direction (green). The top panels (a), (b) and (c) show Re p = 0.02 and the bottom panels (d), (e) and (f) show Re p = 20. The left panels (a) and (d) correspond to φ = 0.05, the middle panels (b) and (e) show φ = 0.2 and the right panels (c) and (f) show φ = 0.38.

FIG. 15 .

 15 FIG. 15. Dependence of the time-averaged particle rotational velocity on the time-averaged vertical particle location. The colors show the particle velocity fluctuations in: the streamwise direction Ω x (blue), vertical direction Ω y (red), and spanwise direction Ω z (green). The top row shows Re p = 0.02 and the bottom row shows Re p = 20. The left panels (a) and (d) correspond to φ = 0.05, the middle panels (b) and (e) show φ = 0.2 and the right panels (c) and (f) show φ = 0.38.

Figure 16 presents

 16 Figure 16 presents the trajectories of the transverse motion of all particles with respect to their initial position in time. This figure delineates the diffusive nature of the particle motion in the transverse direction. The envelope of particle trajectories suggest that at the low volume fraction, φ = 0.05, particles are almost similarly diffusive in z and y directions. At the high volume fraction, φ = 0.38, some particles are clearly more diffusive in the y directions, where wall confinement induces higher particle/particle interactions. By comparing the envelope of particle trajectories, it can be estimated that the particles are approximately 2 to 3 times more diffusive at φ = 0.38 as opposed to at φ = 0.05, with diffusivity being measured as D y = ∆y 2 /(tD 2 γ) and D z = ∆z 2 /(tD 2 γ). The insets in figure16show that the diffusion coefficients D z and D y reach a plateau after some time, where their values are statistically steay. This is the stage where particle motions reach a diffusive regime.The diffusion coefficients, averaged over all particles initially positioned at 2 < y/D < 8 in the domain and averaged over the time interval in the diffusive mode, denoted by D yy and D zz , are shown as a function φ for different Reynolds numbers in figure17. This averaging excludes the particles that are initially located close to walls and therefore have lower diffusion coefficients due to wall effects. Also, in this figure we show comparisons to previous experimental measurements and numerical simulations of the diffusion coefficients in the vertical or shear direction (D yy ), and in the spanwise or vorticity direction (D zz ). To

FIG. 16 .

 16 FIG. 16. Trajectories of particle motion compared to their initial position for Re p = 20. The left panels (a) and (c) delineate φ = 0.05 and the right panels (b) and (d) correspond to 0.38. The insets show the time evolution of the diffusion coefficients, (∆y) 2 /(D 2 t γ) and (∆z) 2 /(D 2 t γ), averaged over all particles.

FIG. 17 .

 17 FIG. 17. Particle diffusion coefficients from the present study and their comparisons to previous studies: (a) D yy and (b) D zz (diffusion coefficients normalized by γ(D/2) 2 ), averaged over all particles initially located at 2 < y/D < 8 and on the statistically steady state diffusive time interval, are shown for different φ and Re p . See figure 5 for particle Reynolds number color coding.

  [START_REF] Mikulencak | Stationary shear flow around fixed and free bodies at finite Reynolds number[END_REF] ,

		8
				Re p = 0.002
		7		Re p = 0.02
				Re p = 0.2
		6		Re p = 2
				Re p = 20
	τ w xy > t /(µ γ)	4 5		Eilers Krieger-Dougherty Boyer et al.
	<	3
		2
		1
		0	0	0.05	0.1	0.15	0.2	0.25	0.3	0.35	0.4
							φ

  . The particle Reynolds number does not have a distinctive effect on the averaged diffusion coefficients for 0.002 < Re p < 2. At lower Reynolds numbers (i.e. Re p ≤ 2), our diffusion coefficients are closer to the experimental measurements of Leighton & Acrivos 12and Breedveld et al.[START_REF] Breedveld | Measuring shear-induced self-diffusion in a counterrotating geometry[END_REF][START_REF] Breedveld | The measurement of the shear-induced particle and fluid tracer diffusivities in concentrated suspensions by a novel method[END_REF] and numerical simulations of Abbas et al.[START_REF] Abbas | Dynamics of bidisperse suspensions under Stokes flows: Linear shear flow and sedimentation[END_REF] , while they are signifi-We conjecture that this difference might either fall within the error bars in the previous experimental measurements, or be caused by the low number of particles available

cantly higher than the Stokesian dynamics simulations of Marchioto & Acrivos 13 and Sierou & Brady

[START_REF] Sierou | Shear-induced self-diffusion in non-colloidal suspensions[END_REF]

. The Stokesian dynamics simulations, that make use of very small Reynolds number assumption, are known to give lower diffusion coefficients

[START_REF] Breedveld | Measuring shear-induced self-diffusion in a counterrotating geometry[END_REF]

. At lower volume fractions, in general our diffusion coefficients are higher than the quantifications in previous studies.