Fluid medium effect on stresses in suspensions of high-inertia rod-like particles

The rheology of suspensions of high-inertia (or granular) non-spherical particles characterized by high particle Stokes and Reynolds numbers is rarely investigated. In this study, we investigate the rheology of suspensions of inertial rod-like particles of aspect ratio 4 subjected to shear ﬂow. In particular, the effect of ﬂuid medium (air, water) against dry granular simulations on the developed stresses is assessed. CFD-DEM simulations are performed for a periodic shear box for a range of shear rates and volume frac- tions of particles. The dependence of rheological properties like shear stress, normal stress difference, pressure and relative viscosity on volume fraction, shear rate, granular temperature and the particle ori- entation are discussed. These results provide insight into the macroscopic rheology of suspensions of rods and demonstrate that the effect of particle shape and surrounding ﬂuid cannot be completely ignored. Air as a ﬂuid medium shows similar scaling as compared to dry granular simulations, but the stress values are generally lower. We observe drastic change in both scaling and values for water as ﬂuid medium. In all cases, the rods show strong alignment in the direction of shear. This study can be further extended to develop stress closures for use in Eulerian ﬂow models. (cid:1) 2019 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


Introduction
Suspensions consisting of particles suspended in fluid are present everywhere in nature (milk, blood, rivers, sandstorms, landslides) and industry (cosmetics, detergents, paints and various fluids in the oil and mineral industries). Suspensions, particularly of solids dispersed in a continuous fluid phase (liquid or gas) exhibit a wide range of behaviour, including Newtonian, shear thinning, shear thickening and discontinuous shear thickening (Mari et al., 2014a). The behaviour of these materials under different conditions is still not well understood and therefore has been a topic of research for decades. Better understanding of the fundamental physics of particle suspensions is necessary to control industrial processes and products.
Studying particle suspensions via experiments is a difficult and tricky process (Campbell, 1989). Rheometers are used to study suspensions, but only offer measurement of shear stress and total pressure in the system. Non-intrusive experimental techniques, which can improve our understanding of suspension rheology, are rare. Also they only measure the combined macroscopic result of all the interactions and not the behaviour of individual particles on a microscopic level. Additionally, most suspension experiments reported in literature are performed on small particles in a liquid continuous phase, meaning that the particle inertia is relatively unimportant (very low Stokes number). In the past 40 years, the discrete element method (DEM), originally developed by Cundall and Strack (1979) for soft-sphere interactions, has been widely used to study granular systems in depth. This method offers many advantages over experimentation, like precise measurements of volume fraction, particle orientation and complete access to the stress tensor and residual kinetic energy of the system. DEM coupled with computational fluid dynamics (CFD) has been used extensively to understand a number of granular suspensions, encountered in shear flows, risers and fluidized beds (Campbell, 1989;Campbell and Gong, 1986;Deen et al., 2007;Mahmood and Elektorowicz, 2016). The CFD-DEM method, due to its ability to accurately deal with a variety of shapes, is also a popular method of choice to study suspensions of non-spherical particles (Mahajan et al., 2018;Zhong et al., 2016).
On the other hand, to understand the rheology of dilute and dense granular flows, granular kinetic theory was developed based on the kinetic theory of dense gases by Chapman and Cowling (1970). Jenkins and Savage (1983) extended the kinetic gas theory to predict the collisional stresses in granular flows. Lun et al. (1984) improved on the theory presented by Jenkins and Savage.
Their theory predicts not only collisional stresses but also kinetic stresses produced due to particle motion. Campbell and Gong (1986) performed 2-D numerical simulations, which were found to be in good agreement with the predictions by Lun et al. (1984). The most widely used kinetic theory models are developed for the flow of inelastic, smooth spheres excluding the effect of interstitial fluid (Fan, 1994;Garzó and Dufty, 1999). Kinetic theory models have been extended or modified to account for different particle properties (e.g. friction) and physical phenomena (e.g. long lasting collisions) (Campbell, 2002;Jenkins, 2006Jenkins, , 2007Jenkins and Berzi, 2010;Chialvo and Sundaresan, 2013).
Another popular theory in literature to describe granular flows is l eff (I) rheology (MiDi, 2004;Forterre and Pouliquen, 2008), which was developed in the mid 2000s. In this theory, constitutive equations, which describe the conservation laws of granular flows, are characterized by an inertial number, I ¼ _ cD= ffiffiffiffi P q p q , with _ c the shear rate, D the particle diameter, P the particle pressure, and q p the particle density. An important result of the theory is that the effective friction (l eff ), which is the ratio of shear stress to the pressure, is a function of inertial number I. These functions have been determined by simulations (Jop et al., 2006;Da Cruz et al., 2005;Hatano, 2007) as wells as by experiments (Fall et al., 2015;Jop et al., 2005). Most of the studies discussed above deal with circular disks in 2D or spheres in 3D, primarily because dealing with nonspherical particles in kinetic theory is not trivial. However, there is a recent increase in studies dealing with such particles. Researchers have performed experiments (Baosheng et al., 2010;Li et al., 2004) and simulations with other particle shapes. Simulations performed by Pena et al. (2007) for polygonal disks and by Cleary and Sawley (2002) for blocky particles show that particle shape has strong effects on the granular temperature and volume fraction in the core of the flow, both of which are smaller compared to volume-equivalent spherical particles. Moreover, particles with non-regular shape (higher aspect ratio or more angular geometry) are hard to shear due to interlocking of the particles. One of the main attributes of elongated particles is that they have a preferred alignment towards the main flow stream, as shown by Pena et al. (2007), Reddy et al. (2009Reddy et al. ( , 2010 and Campbell (2011). Campbell (2011 found that for smooth frictionless ellipsoid particles, smaller stresses were observed when compared to volume-equivalent spheres. It was also found that large surface friction can lead to particle rotation which can block the flow leading to higher stress values. Guo et al. (2012aGuo et al. ( ,b, 2013Guo et al. ( , 2015 extensively studied the dry granular shear flow of rod-like particles of various aspect ratio, wet flexible fibers, flat disks and elongated rods in the presence and absence of friction. Nagy et al. (2017) have also performed a 3D simulation with spherocylindrical particles up to aspect ratio of 2.5. Their study demonstrated that l eff I ð Þ rheology can be extended to non-spherical particles at high volume fractions. All of these studies exclude the effect of fluid on the developed stresses. Guazzelli and Pouliquen (2018) explained, in a comprehensive paper on rheology of dense granular suspensions, that hydrodynamic interactions or lubrication forces between the particles are important in the dilute regime. At high concentrations they become less significant compared to direct particle contacts. Although in reality granular particles are usually surrounded by an interstitial fluid (like air) the influence of the latter on the dynamic properties of the solid particles is generally neglected in most theoretical and computational works. Needless to say, the effect of the interstitial fluid on solid particles turns out to be significant in a wide range of practical applications and physical phenomena (Xu et al., 2013;Chamorro et al., 2015). Garzó et al. (2012) presented an Enskog kinetic theory for monodisperse gas-solid flows. Their theory demonstrates that the effect of the fluid phase on the constitutive equations for the solid phase shear viscosity is non-negligible. However, their theory, like most kinetic theories, is limited to spherical particles and therefore is not applicable to granular non-spherical particle suspensions.
While the majority of studies so far have focused on Stokes c is the shear rate, g f the viscosity of the fluid phase, and q p and d p are the density and volume equivalent diameter of the solid particles, respectively) of the order of 1 or smaller, this paper focuses instead on particle suspensions that have Stokes number and Reynolds number much larger than 1. The range of numbers is given in Table 1. Although such high Stokes and Reynolds numbers are encountered in a number of applications (a variety of chemical, pharmaceutical and process industries), the main motivation behind understanding the high Stokes number flows comes from the use of biomass particles as raw material in fluidized bed gasifiers. These biomass particles are often preprocessed into elongated pellets with aspect ratio ranging from 2 to 8 before being fed to gasifier. As a starting point, the shape is simplified to spherocylinders and an intermediate aspect ratio 4 is chosen for this study. However, this work can be extended to any aspect ratio mentioned above. The pellets typically have a size of a few millimeters and are processed in the fluidized bed gasifier with an upward flow of gas. This combination of large particle size and low viscosity gas results in high Stokes numbers. More generally, fluidized beds containing high inertia particles are used in a variety of chemical, pharmaceutical and process industries and often experience the range of Stokes and Reynolds numbers studied in this paper. The values for the water case given in Table 1 are consistent with values encountered in water treatment plants. There, activated carbon particles, nonspherical in shape (often elongated) are fluidized in water. These particles have a volume equivalent diameter of a few millimeters resulting in the Stokes numbers studied in this work.
In summary, the rheology of particulate suspensions is highly dependent on the properties of the particles (size, shape and surface properties) as well as the properties of the continuous fluid (viscosity, relative density). Understanding the interplay of interparticle forces and coupling to the fluid is crucial for understanding the rheological behaviour of granular suspensions, and producing meaningful results. Most existing studies assume mono-disperse spherical particle granular systems to simplify the analysis. However, real systems nearly always differ from this ideal case, being poly-disperse or containing ellipsoids, platelets or even irregular shapes. The relatively few works on the rheology of nonspherical granular particles in the past have dealt mainly with dry frictionless elastic particles. There are very few studies available on the rheology of non-spherical and inelastic granular particle suspensions via CFD-DEM simulations. The present study focuses on exploring this topic, in particular for rod-like granular particles. Not only the rheological properties of rod-like particles are investigated, but also the theories discussed above are used to understand the observed behaviour. O(10 2 À 10 5 ) O(10 0 À 10 3 )

Model framework
The modelling done in this work is performed using a combined CFD-DEM method. A modified version of the CFDEM code is used to perform the simulations (Goniva et al., 2012). The individual particles are tracked and their trajectories are numerically integrated over time. The local contact forces and torques develop when adjacent particles spatially overlap. Fig. 1 shows an example of an overlapping contact between two spherocylinder particles. Details about the contact force, torque and overlap calculations can be found in our previous paper (Mahajan et al., 2018).
Consider a spherocylinder p in a dense environment of similar spherocylinders. The translational motion for spherocylinder p can be calculated by integrating where F ! p;c is the contact force due to collisions and F ! p;f !p is the total interaction force exerted by the fluid phase on the particle. The rotational motion of particle p is calculated by integrating where I ! !
p ; x ! p and T ! p are the moment of inertia tensor, angular velocity and torque for particle p, respectively. There are a number of models in the literature to calculate contact forces between particles. The one which is mostly used, and one of the simplest, is described by Cundall and Strack (1979), where the normal contact force between particle i and j is calculated by a linear springdashpot model as follows: where k n is the normal spring constant, d ij;n is the normal overlap, f n is the normal damping coefficient and v ! ij;n is the normal relative velocity between the particles i and j at the location of the contact point. The reader is referred to the paper of Mahajan et al. (2018) for more details about the contact model. Note that in this paper we limit ourselves to smooth particles without tangential friction. This is closest to the assumptions of (most common forms of) kinetic theory. Additionally, by reducing the number of parameters (in this case, friction) we can independently study the effect of fluid on the rod suspension. In a forthcoming paper, we investigate the additional role of tangential friction.
The fluid phase is described in an Eulerian manner by imposing a mesh of equal sized cells on the simulation domain. The PISO algorithm (Versteeg and Malalasekera, 1995) is used to solve the phase continuity and momentum transport equation for incompressible, Newtonian flow. Eqs. (4) and (5) are known as Model A in the literature as proposed by Zhou et al. (2011).
where is the fluid volume fraction, q f the fluid density, u ! f the fluid velocity, and s ! ! f is the fluid stress tensor. f ! p!f represents the momentum exchange between the fluid and particle phases. The volumetric force acting on the fluid phase by its interaction with the particles is given by where F ! p;f !p is the force acting on an individual particle p and n is the total number of particles in the CFD cell. The details of the transient solver with a PISO loop for pressure velocity coupling used in this work are given in the paper by Goniva et al. (2010). Hölzer et al. (2008) derived an equation describing the drag coefficient for a single non-spherical particle in a gas flow (Eq. (7)). Here / is the particle sphericity, v ! r is the relative velocity between the particle and the gas, This equation incorporates the orientation of the particle in the crosswise (/ ? ) and lengthwise sphericity (/ k ). These are calculated based on the angle between the gas velocity vector and the particle orientation vector, otherwise referred to as the incident angle or angle of attack.
Di Felice (1994) developed a correlation describing the effect of local fluid volume fraction on the drag force (Eqs. (10) and (11)).
It must be noted that the effect of torque and lift on particles due to the fluid are neglected in this work. The lubrication force, which arises from the additional hydrodynamic pressure when the interstitial fluid is squeezed out from the space between two solid surfaces, is also not accounted for. This is acceptable for large granular particles with higher Stokes numbers. The only force acting on the particles due to the surrounding fluid medium is the drag force acting at the center of mass of each particle. Therefore, it is safe to say that the orientation of particles in the simulations are solely a result of particle-particle and particle-wall collisions.
Our prime quantity of interest is the total particle stress tensor r ! ! t . Due to its dynamic nature, momentum transport can occur by simple flux of particles or by particle-particle contact. Therefore, the particle stress is a combination of two independent contributions, i.e. the streaming stress tensor ( r ! ! s ) due to particle momentum flux and the collisional stress tensor ( r ! ! c ) due to particle collisions (Cleary, 2008). These stress tensors are given as follows: where v ! 0 i is the velocity fluctuation for particle i and . . . h i indicates averaging over time (under steady state conditions) and over all particles located in the bulk control volume V cb . The velocity fluctuation is calculated relative to the locally averaged velocity of the To avoid wall-effects, the measurements are performed in a bulk-like control volume V cb which is one particle length away from the walls located at z ¼ 0 and z ¼ h d . The granular temperature (T) is a measure for the residual kinetic energy in the system. It is calculated from the velocity fluctuation of the particles.
Based on the required volume fraction, a number of particles were placed inside the computational domain. Shear in the domain is achieved by moving the wall at z ¼ 0 in the negative x-direction and the wall at z ¼ h d in the positive x-direction (see Fig. 3). The no-slip boundary conditions at the walls also drive the fluid. The boundary conditions are periodic in the x and y directions. The gravitational force acting on the system is set to zero. Note that some amount of wall-induced particle ordering occurs. This influence was found to quickly decay away from the walls. The particles used in this study are smooth (frictionless), hard (large contact stiffness), and dissipative (restitution coefficient of 0.7). Campbell (2002) identified the importance of the elastic properties of particles in determining the overall rheology of dense granular flows. He studied the effect of interparticle contact stiffness on developed stresses. The results showed that the stresses scale with the stiffness saturating at sufficiently high stiffness, indicating that the stresses are generated by the particle elasticity. Therefore, in this study, the particles are assumed to be hard with large contact stiffness, frictionless with no tangential force and partly inelastic. There is however a price to pay for simulating very large contact stiffness: the timestep necessary for accurate integration of the equations of motion becomes very small (see Fig. 2). Moreover, higher values of restitution coefficient demand even smaller timesteps, making simulations extremely expensive.
The simulation domain needs to be large enough to produce sufficient statistics and results which are independent of domain size, while at the same time keeping the computational cost minimal. Therefore, a comparison was done for two domain sizes: 6; 4; 4 ð ÞÂ L p , and 6; 4; 8 ð ÞÂ L p . The results were checked for dependence of the stress on the domain size. The two different-sized domains showed good agreement for all measured parameters: stress components, pressure and granular temperature. Therefore, the smaller domain was employed in all simulations. All simulation parameters are reported in Table 2.
Under a shear flow, the system achieves steady state when an equilibrium is reached between the energy supplied by the driving walls and the energy dissipated by interparticle collisions and fluid friction. For each simulation, the time required to achieve a steady state is different and is dependent on the initial configuration, volume fraction and shear rate. All quantities reported in this paper are obtained after a steady state is achieved, typically for a shear strain of at least 10. The computing times depend primarily on the overall frequency of particle collisions which depends on volume fraction and applied shear rate. For low volume fractions and low shear rates, simulations are faster as compared to high   volume fractions and high shear rates. On an average, each simulation (for each shear rate and volume fraction) took 2 days on 24 processors. The stresses and granular temperature initially increase due to initial particle collisions and then remain at a nearly constant level. Note that the steady state is independent of the initial configuration of the particles. The initial configuration only affects the shear strain within which the steady state is achieved. Loisel et al. (2015) studied particle migration towards channel walls due to the Segré-Silberberg effect, giving rise to a nonuniform particle concentration profile in the z-direction. This phenomenon is also observed in our simulations. Therefore, the actual volume fraction is measured in the bulk region, one particle length away from the z-walls, called the measurement domain henceforth. Just like the volume fraction, the shear rate for the measurement domain is different than the apparent shear rate applied by the z-walls. The average particle velocity profile in the measurement domain was used to measure the actual shear rate. We explicitly checked and confirmed that the apparent shear rate _ c was constant (i.e. a linear shear flow profile) across the measurement domain for all runs.
The streaming stress tensor (Eq. (13)) can be written more explicitly as (Campbell and Gong, 1986): For a sheared system as used in this study, the opposite velocities at z and Àz walls makes the global average velocity approximately zero. To make sure the velocity fluctuations are correctly calculated, the shear box is split into vertical bins (bins of equal height Dz). Fluctuation of the particle velocities are calculated relative to the average velocity in each bin and averaged over the measurement domain. Test simulations have been performed to find out the minimum number of bins necessary to estimate fluctuations independent of the number of bins in the domain. Based on this study, 16 bins were used to determine the granular temperature.

Model verification
In this section, the implementation of our code and stress measurements is verified by comparing shear flow simulations of smooth granular spheres against the predictions of kinetic theory. As discussed by Guo et al. (2012a) and Campbell (2002), the kinetic theory can be safely applied to dilute and moderately dense sys-tems in which binary instantaneous collisions are dominant. However, it might fail for very dense systems where multi-particle long-lasting collisions become significant and the stresses show a dependence on the interparticle contact stiffness. For validation, our simulations are compared with kinetic theory at low volume fraction. The spherical particles have a diameter such that they have the equivalent volume as the spherocylindrical particles studied in the rest of this paper (d p ¼ 0:0053 m).
The results are compared with the kinetic theory predictions by Lun et al. (1984). According to their theory, the steady-state shear and normal stresses for smooth, hard, and slightly dissipative spheres in the plane of shear flow is given by: a is a function of the coefficient of restitution e, while g 0 is radial distribution function for spheres which depends on the solid volume fraction () (Carnahan and Starling, 1969). The predicted steady-state granular temperature (T) is given as follows: where F ; e ð Þ can be split into a collisional F c ; e ð Þ and a streaming The pressure in the system is a result of generated normal stresses and is measured as 1=3 rd of the trace of the total stress tensor. As can be seen in Fig. 4, the measured granular temperature and pressure both show good agreement with the granular kinetic theory developed by Lun et al. (1984). The slight mismatch, notably the slightly lower granular temperature, can be attributed to the fact that the simulations were carried out for particles which have somewhat dissipative collisions. As pointed out earlier, the assumption of kinetic theory is that the particle collisions are nearly elastic (1 À e ( 1). The comparison for collisional and streaming stresses also show a good match with theoretical predictions, thus verifying the numerical implementation of the model and measurements (see Fig. 5).

Results and discussion
In the following sections, the behaviour of granular spherocylinder particle suspensions is investigated, looking at shear and normal stresses, pressure, granular temperature, particle orientation, apparent friction and viscosity. The results are shown for three cases: (a) dry granular shear flow (in the absence of fluid) where only particle-particle interactions shape the rheology, (b) an air suspension where elongated particles are suspended in an airlike gas and (c) water suspension where particles are suspended in a water-like fluid. For reference purpose, the predictions from granular kinetic theory (Eqs. (17), (19)-(21)) for volumeequivalent spheres are shown in red. Our aim by presenting these comparisons is not just to highlight the difference in the measured values against granular kinetic theory but also to highlight the similarity of the scaling with shear rate observed for rods and that predicted by kinetic theory for spheres.

Granular temperature and pressure
In this section, the results of granular temperature and pressure as a function of shear rate are presented for 5 volume fractions. The granular temperature is proportional to the particle fluctuating velocity squared. As the solid volume fraction increases, the collisional frequency increases, leading to more dissipation and hence lower granular temperature. This prediction is in accordance with our simulation predictions for dry granular flow (see Fig. 6a). The granular temperatures are highest for low volume fractions since they have more space to move in between collisions. Also, the granular temperature is lower for spherocylindrical particles than for volume-equivalent spheres because their projected area in the direction of flow is smaller. Similar behaviour is seen under the influence of air, as shown in Fig. 7a. However, the magnitude of the granular temperature is consistently lower. Even though the particle inertia is high relative to that of air, the air still acts as a momentum sink and takes away some fluctuating energy from the particles. Fig. 8a shows the granular temperature under the influence of water. As expected, water strongly dampens the velocity fluctuations and therefore the resultant granular temperature. In the presence of water, the granular temperature shows the reverse behaviour as a function of solid volume fraction, when compared with the dry granular and air cases. The granular temperature is higher for higher volume fraction of the particles. Due to the higher exchange of momentum with water (as compared to the dry granular and air cases), the particles loose their fluctuating energy fast and only gain more fluctuations when they collide with other particles. In fact, when velocity fluctuations are dampened so strongly, the particles mainly collide because they are forcefully brought into contact by the shear flow. Under such conditions, a higher volume fraction leads to more frequent collisions, leading to a higher granular temperature.
In this study, just like the granular temperature, pressure evolves and is an emergent property of the system (at the given shear rate and solid volume fraction). The pressure shows the expected behaviour: it is increasing with the shear rate and solid volume fraction. A higher pressure at higher shear rate is a direct result of the increase in granular temperature, while a higher pressure at higher volume fraction is a direct result of more collisions per unit time. The measured pressure values are compared with the pressure predicted by kinetic theory for volume-equivalent spheres. It is found that the measured pressure for rod-like particles does not agree with kinetic theory predictions for spheres:  (Lun et al., 1984). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) for the dry and air cases the measured pressure is lower, while for the water case the measured pressure is lower at low volume fractions, but higher at high volume fractions. When comparing spherical and rod-like particles at volume fraction 0.1 (Figs. 4a and 6a), it is clear that rod-like particles attain lower granular temperature and pressure. This effect is a result of shape of these particles and preferred alignment.

Collisional and streaming shear stress
Next the shear stresses for spherocylinder particles are investigated. The collisional and the streaming components of the shear stress are plotted as a function of shear rate for 5 volume fractions in Figs. 9-11. Fig. 9a shows that the kinetic theory for spherical particles agrees well with the measured collisional stress of dry spherocylinder particles at intermediate volume fractions. However, this must be a fortunate cancellation of effects, because Fig. 6a showed that the granular temperature was underpredicted. Moreover, the deviation is large for lower and higher volume fraction. Fig. 9b shows that the streaming stresses are lower for spherocylinder particles than predicted by kinetic theory of spherical particles. This is due to the preferential alignment of the particles, as will be shown in Section 4.5. The preferential alignment allows particles to avoid intense collisions. Similar behaviour was also reported by Guo . When oriented in the direction of shear, the projected area of elongated particles in the direction of flow is smaller compared to volume-equivalent spherical particles. As a result, particle collisions are less probable, resulting in less fluctuations (as evidenced by the lower granular temperature in Fig. 6a) and therefore lower stresses.
Compared to the dry granular simulations, the suspensions in a fluid medium (Figs. 10 and 11) show different stress levels. A higher fluid viscosity translates to less intense collisions. Thus, the fluid medium takes away some of the fluctuating energy from the particles. The collisional stresses are slightly lower in the presence of air (Fig. 10a) when compared to the dry granular case. However, the effect is small because the particles are still highly inertial. Also the streaming stresses are considerably lower in the presence of a fluid because the fluid medium reduces the velocity fluctuations and thereby acts as an additional dissipation channel.

Apparent friction and viscosity
The shear stress measurements can be interpreted in two different ways, either as an apparent coefficient of friction or as a quasi-Newtonian (shear rate dependent) total suspension viscosity. The apparent friction coefficient is defined as the ratio of the total particle shear stress to particle pressure (l s ¼j r t j =P) where r t ¼ r c þ r s is the total particle shear stress. Fig. 12a shows the apparent friction coefficient for dry granular shear flows at different volume fractions. Kinetic theory predicts that for spherical particles the apparent friction coefficient does not change with shear rate and has a non-monotonous dependence on volume fraction (red symbols). For spherocylinder particles, our results show that the apparent friction coefficient has a dependence on shear rate, but that this dependence is only weak, especially in view of the fact that the measured stresses and pressure vary over many orders of magnitude. Contrary to spheres, the apparent friction coefficient for spherocylinder particles displays a monotonous (and stronger) decrease with increasing volume fraction. The generally lower apparent friction coefficient for spherocylinder particles is due to the alignment of particles at higher volume fractions, leading to less resistance to the flow. Due to the absence of tangential particle-particle friction, the smooth surface allows particles to slide past each other, therefore reducing the effective (macroscopic) friction coefficient. The fluid medium effect on the apparent friction coefficient can only be seen in the case of water (Fig. 14a). This effect is largely caused by the wide variation in pressure with volume fraction as shown in Fig. 8b. Comparatively, the shear stresses (Fig. 11) do not change as much with volume fraction, resulting in considerable change in apparent friction.
The total suspension viscosity is calculated as the sum of particle viscosity and fluid viscosity: It is generally observed that non-Brownian particles suspended in a Newtonian fluid raise the viscosity of the suspension. It is also observed that such a suspension has a shear-rate dependent rheology (see e.g. Mari et al. (2014a)). The granular and air cases in Figs. 12b and 13b clearly show shear thickening behaviour where the total suspension viscosity increases with increasing shear rate. This is seen for a wide range of shear rates and volume fractions. The shear thickening behaviour can be associated with the inertial nature of the particles. A number of studies have previously associated particle inertia to the existence of shear thickening behaviour in suspensions (Fall et al., 2010;Kawasaki et al., 2014). Particularly, for suspensions with particle Stokes number larger than 1, shear thickening due to inertia has been widely reported (Fernandez et al., 2013). Perhaps surprisingly, for the case of water almost Newtonian behaviour is seen for high volume fractions while low volume fractions still exhibit weak shear thickening behaviour. At higher volume fraction, particles slide past each other due to preferential alignment of particles, consequently avoiding head-on types of collisions. At the same time, due to the similar densities of particles and fluid, the particles hardly affect the flow behaviour of the overall suspension. It must be noted that this behaviour is a special case and a direct result of the small density difference between fluid and particles. Another key observation here is that the viscosity is increasing with shear rate much more strongly for air than for water case. This is a direct consequence of the dominance of particle stresses shown in the previous section. The collisional particle stresses for water, shown in Fig. 11a, vary over 3-4 orders of magnitude, while for air shown in Fig. 10a, they vary over 6 orders of magnitude.

Normal stress
Under shear flow conditions, there are two independent normal stress differences, the so-called first and second normal stress difference. Normal stress differences are important because together with the shear stress they completely define the shear rheological behaviour of the suspension. In other words, these constitute all the rheological information that can be obtained from measuring the stress components under shear flow. Non-zero normal stress differences are generally an indicator of non-Newtonian behaviour of the suspension. For shear flow in the xz-plane (Fig. 3), they are defined as follows: Their sign and value are dependent on the particle micro-structure of the suspensions. For the case of non-inertial non-Brownian sphere suspensions, N 2 has a negative value, and this value is increasing in the negative direction as the volume fraction of spheres increases (Zarraga et al., 2000;Sierou and Brady, 2002).
As for the case of N 1 , there is still an intense discussion in the scientific community, especially about its sign. For example, Sierou and Brady (2002) showed that the sign of N 1 is always negative and always increases with increase of volume fraction and that hydrodynamic interactions are the main contributor to its value. Instead, Mari et al. (2014b) found that the sign of N 1 changes from negative to positive, especially at volume fractions near the maximum packing fraction. This is due to the formation of frictional contacts, and the transition from hydrodynamically dominated N 1 to frictional dominated N 1 . For the case of non-spherical particles, the literature is more sparse. What is known is that for fiber and spherocylinders of sufficiently high aspect ratio, the magnitude of N 1 is greater than that of N 2 . Specifically, it is predicted that N 2 % À0:5N 1 . This relation is holding both for rigid fibers in the Stokes flow regime (Snook et al., 2014) and for inertial granular spherocylinders (Nagy et al., 2017), and is found to be independent of the aspect ratio up to 30 (Snook et al., 2014). The ratio À N 2 N 1 for the three different cases of granular, air, and water is shown in Fig. 15 for different volume fractions. For the case of pure granular flow (Fig. 15a) our results are similar to that of Nagy et al. (2017). The main difference between Nagy et al. (2017) and our work is that for our simulations the normal stress differences have two contributions: one streaming and the other one collisional, while Nagy et al. (2017) only considered the collisional component. This shows that the requirement to get N 2 % À0:5N 1 is that the collisional stress has to be dominant. Fig. 15b and c show that also in the presence of a fluid (air or water), the relation N 2 % À0:5N 1 still holds, which is because also under these conditions the collisional stress is the dominant one. Our results confirm that the specific relation between N 1 and N 2 holds in a regime which lies between the purely noninertial viscous (St ( 1; Re p ( 1) and purely inertial (St ) 1; Re p ) 1) regimes. It will be interesting to see if the same role holds when another interparticle, force such as friction, is added.

Flow alignment angle
One of the well known features of granular shear flows of elongated particles is the preference in alignment with the flow direction. This alignment of the major particle axis along the flow direction allows the system to achieve a low energy state by reducing the obstruction to the flow of particles. In a non-isotropic system, like in the present study, the eigenvector belonging to the largest eigenvalue of the order tensor preferential orientation of the particles. The angle between this eigenvector and the flow direction (x) is the flow alignment angle, which can also be expressed directly as Þ . Fig. 16 shows the steady state flow alignment (also called the extinction angle) for spherocylinder particles under shear flow. Fig. 16a and b show that for the dry and air cases the flow alignment angle is relatively unaffected by the shear rate (within the range of shear rates studied here). The flow alignment angle decreases strongly with increasing solid volume fraction, as shown in Fig. 16d. These observed alignment angles are in good agreement with the predictions by Nadler et al. (2018). The decrease of the alignment angle with volume fraction can be explained in terms of the reduction in free volume available for random reorientations, which can be appreciated from the snapshots in Fig. 17.
Finally note that for the water case, the alignment angle is smaller than for the other 2 cases. This is a consequence of the strong damping by the fluid phase, leading to a suppression of random reorientations by fluid interactions, rather than interactions with neighbouring particles.

Conclusions and outlook
Simulations of inertial spherocylinder particles under shear flow were performed. The rheological behaviour for dry granular flow and for particles in two fluids, air and water, were studied. 3D CFD-DEM simulations were performed to measure the generated granular temperature, pressure, and stress tensor components for particles of aspect ratio 4. In this work, the effect of friction was neglected. The main contribution of this work was to demonstrate the effect of shape of the particles and fluid medium on the generated particle stresses. Compared to volume-equivalent spheres, a drastic drop in magnitudes of shear stresses for elongated particles was observed. It is evident that particle shape has an influence on the developed stresses and therefore the rheology. We demonstrated that elongated particles tend to align along the shear direction. This alignment allows the system to attain a lower energy state by avoiding intensive collisions. It was also demonstrated that the apparent (or macroscopic) friction coefficient does not change with shear rate for dry granular flows and that it decreases with increasing volume fraction. For rod-like particles suspended in air, the results show a strikingly similar scaling as predicted by the granular kinetic theory for volume equivalent spheres, but with somewhat lower stress values. This explains the validity of previous researchers in using granular kinetic theory for simulating fluidized beds while neglecting both the effect of interstitial fluid and particle shape. Such similarity of scaling is not found for the case of water.
Note that hydrodynamic torque and lift forces were not accounted for in this paper. While there are a few torque and lift models available in literature, the accuracy and applicability of these models to rod-like particles is questionable. We expect that for the volume fractions investigated here, the particle-particle collisional interactions are much stronger than hydrodynamic torque and lift. Additionally, once the steady state is achieved almost all particles are oriented in the direction of shear. Thus the resulting torque and lift forces are even smaller in magnitude. It should be acknowledged that these forces are probably non-negligible for the case of water. We expect a change in results for the case of water, if simulations are re-done including accurate expressions for lift and torque in multi-particle environments, when such correlations become available.
In our future work, we will take into account the effect of particle-particle friction on the shear flow properties of granular rod-like particle suspensions. Together, these studies will help us understand the rheology of rod-like inertial particles suspension. Furthermore, they can be used to develop stress closures for use in coarse-grained models of fluidized bed biomass gasifiers, which typically contain rod-like particles.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.