How do the finite-size particles modify the drag in Taylor-Couette turbulent flow

The TC facility used in the present work is shown by the sketch in figure 1(a). The inner cylinder is made of aluminum, while the outer one is made of PMMA (polymethyl methacrylate) to perform optical measurements. The outer cylinder is surrounded by a PMMA cubic circulating bath, which keeps the system temperature at $T=22\pm 0.1$ ℃ during the experiments. The radius of the inner and outer cylinder is $r_{i}$ = 25 mm and $r_{o}$ = 35 mm, respectively, giving the gap width $d$ = $r_{o}-r_{i}$ = 10 mm and the radius ratio $\eta=r_{i}/r_{o}=0.714$. The height of the inner cylinder is $L$ = 75 mm, which gives the aspect ratio of the system $\Gamma=L/d=7.5$.

During the experiments, the outer cylinder is fixed, while the inner cylinder is connected to and driven by a rheometer (Discovery Hybrid Rheometer, TA Instruments), and the overall torque of the system is measured at the same time by the torque sensor of the rheometer. The control parameter of TC flow is the Reynolds number defined by:

where $\omega_{i}$ is the angular velocity of the inner cylinder and $\nu$ is the kinematic viscosity of the fluid. In the present work, the $\omega_{i}$ ranges from $50$ to $200$ rad/s, which gives the $Re$ ranges from $6.5\times 10^{3}$ to $2.6\times 10^{4}$.

The overall torque of the system, $\tau_{total}$, that needed to drive the inner cylinder at a constant angular velocity $\omega_{i}$, can be divided into two parts: (i) the torque, $\tau$, due to the sidewall of the inner cylinder (TC flow), (ii) the torque, $\tau_{end}$, contributed from the von Kármán flow at the end plates. Since the aspect ratio of the system is relatively small here, the rotation of the end plates generates secondary vortices that contribute to the torque measurements (Bagnold, 1954; Hunt et al., 2002). The calibration of the end effect due to the von kármán flow between the bottom and top lids of the inner and outer cylinder can be found in Appendix A. Here, only $\tau$ is used, and it can be non-dimensionalized as:

where $L$ is the length of the inner cylinder and $\rho_{f}$ is the fluid density.

In the experiments, the gap between the bottom and top lids of the inner and outer cylinders is 2 mm and larger than the particle diameter, the particles therefore would enter in and escape from the bottom and top gaps between the inner and outer cylinders. For each experiment, the system is initially started from a static state, and then the inner cylinder is imposed a constant angular velocity by the motor of the rheometer. Before measuring the torque, the inner cylinder is kept rotating to ensure the establishment of a statistically stationary state. For each $\omega_{i}$, the torque measurements are repeated 3 times, and the averaged value is used, of which the standard deviations are less than $1\%$.

The particles we used are printed with a 3D printer using photosensitive resin (Elastic Resin of Formlabs. Inc) and are then cured at 60 ℃ for 60 mins. To check their water-absorbing ability, the particles are immersed in the static glycerin-water solution for 2 hours, and afterward the change in weight is found to be negligible. The density of the particle is $\rho_{p}=1.06\times 10^{3}$ kg/m³ (averaged value based on more than 1000 particles). The particle volume fraction, $\Phi$, ranges from $0\%$ to $10\%$ with a spacing of $2\%$. In addition, optical measurements and additional torque measurements are performed at $\Phi=0.5\%$. To study the shape effect of the finite-size particles, we change the particle aspect ratio, $\lambda=h_{p}/d_{p}$, where $d_{p}$ is the length of the symmetric axis of the particle, and $h_{p}$ is the length that perpendicular to it. In this work, $\lambda$ is designed to be equal to 1/3, 1, and 3, corresponding to the oblate, spherical and prolate, respectively (see figure 1(b)). We fix the volume of each particle the same for three cases of aspect ratio. The equivalent diameter of the particle based on the volume, $d_{e}$, is equal to the diameter of the spherical particle, i.e., $d_{e}=d_{s}=1$ mm. Note that the particles are verified to be finite-size particles. For example, the dissipative length scale of the TC flow, $\eta_{K}$, can be estimated as around $0.075$ mm when $\omega_{i}=50$ rad/s, which gives $d_{e}\textgreater 10\eta_{K}$ and satisfies the assumption of the finite-size model (Voth & Soldati, 2017).

To eliminate the effect of the gravity and the buoyancy force, the density of the particle and the fluid are matched using glycerin-water solution (the weight fraction of the glycerin is $w_{t}=25\%$). The solution density is $\rho_{f}=1.0598\times 10^{3}$ kg/m³ and the kinematic viscosity is $\nu=1.913\times 10^{-6}$ m²/s at $T=22$ ℃  (Cristancho et al., 2011), giving a 0.02% density mismatch from the particle density. Since the temperature of the system is well controlled by the circulating water bath, the effects of the temperature variation on the density and viscosity of the solution can be neglected in the present work.

The particle motion is recorded by a CCD camera (MD028MU-SY, Ximea.Inc) at a frequency $f_{c}=15$ Hz, and then the particle positions are obtained by performing particle detection. The flow domain is illuminated by a light source, which is shown in figure 1(a). The details of the detection methods and examples can be found in Appendix C.

Firstly, we study the effects of changing the Reynolds number on the drag of the system. For the global analysis of the drag modification by bubbles or particles, the friction coefficient, $c_{f,\Phi}$, which evaluates the skin frictional drag of the flow system, has been extensively used in various flow geometries (Sanders et al., 2006; Verschoof et al., 2016; Olivucci et al., 2021). In the TC system, the friction coefficient is defined as (Lathrop et al., 1992; Van Gils et al., 2011):

Figure 2 shows the normalized friction coefficients of the system, $c_{f,\Phi}/c_{f,\Phi=0}$, where $c_{f,\Phi=0}$ is the friction coefficient of the single-phase case. It is found that the drag modification increases with increasing $\Phi$ for a given $Re$. In the investigated ranges, the greatest drag enhancement is found to be around $20\%$, which is observed in the spherical case at $Re=6.5\times 10^{3}$ and $\Phi=10\%$ (figure 2(b)). Additionally, $c_{f,\Phi}/c_{f,\Phi=0}$ decreases as the $Re$ increasing, which has also been found and is referred to as the shear-thinning effect in the suspension of deformable particles (Adams et al., 2004) and the emulsions (Yi et al., 2021; Rosti & Takagi, 2021). Note that, the overall drag of the system ($G$ or $\tau$, which are not shown here), increases with increasing $Re$ and is consistent with the literature (Stickel & Powell, 2005; Fall et al., 2010; Picano et al., 2013).

The modification of the friction coefficient is related to the changes of the wall stress, which can be decomposed into three parts when no external force and torque is applied on particles:

where $\tau_{v}$ the viscous stress of the fluid phase, $\tau_{p}$ the particle-induced stress, and $\tau_{T}=\tau_{T_{f}}+\tau_{T_{p}}$ the turbulent Reynolds stress of the combined phase, $\tau_{T_{f}}$ and $\tau_{T_{p}}$ the turbulent Reynolds stress of fluid and particulate phase, respectively. Following the studies in Zhang & Prosperetti (2010) and Picano et al. (2015), the explicit expression of each term above in a Couette flow (Batchelor, 1970; Wang et al., 2017a) can be written as:

where $\mu$ the dynamic viscosity of the fluid, $U_{f}$ the mean fluid velocity in the azimuthal direction, $u^{\prime}$ and $v^{\prime}$ the velocity fluctuation in the azimuthal and radial direction with the subscripts $f$ and $p$ denoting fluid and particulate phase, respectively, $\sigma^{p}_{xy}$ the general stress in the particulate phase, normal to the cylindrical surface and pointing in the radial direction. Here we note that the contribution to fluid viscous stress due to velocity in the axial direction is neglected since it is a minor effect.

At a fixed $Re$, the wall stress of single-phase case, $\tau_{w,\Phi=0}=\mu\frac{\mathrm{d}U}{\mathrm{d}r}|_{r=r_{i}}$, is a constant, and the particle-induced stress ($\tau_{p}$) increases with increasing $\Phi$. In addition, it has been reported that in channel flow (Picano et al., 2015), for a fixed $Re$ in the current low particle volume fraction ranges, the viscous stress ($\tau_{v}$) weakly depends on $\Phi$, and the turbulent Reynolds stress of the combined phase ($\tau_{T}$) increases with increasing $\Phi$. Though there are differences between the channel flow and the TC flow, the result in the channel flow is instructive for the understanding of the current study. Considering the dependence above of the stress on $\Phi$ at a fixed $Re$, the normalized friction coefficient increases with increasing $\Phi$, which is consistent with the trends in figure 2. On the other hand, for a fixed low $\Phi$ (as in this work), the turbulent Reynolds stress of fluid phase in single-phase flow ($\tau_{T_{f},\Phi=0}$) increases greatly as $Re$ increases, which makes the contribution of particulate phase ($\tau_{T_{p}}$ and $\tau_{p}$) relatively insignificant and therefore results in the decreasing trends of normalized friction coefficients with increasing $Re$. Another feature of the normalized friction coefficients in figure 2 is that, as $Re$ increases, the differences between various $\Phi$ decrease. Specifically, when $Re$ is small, the contributions of particulate phase to the wall stress are non-negligible, hence the normalized wall stress (friction coefficients) highly depends on $\Phi$ through $\tau_{T_{p}}$ and $\tau_{p}$. However, when $Re$ is high enough and $\Phi$ is low (as in this work), the turbulent Reynolds stress due to fluid phase in particle-laden flow and single-phase flow ($\tau_{T_{f}}$ and $\tau_{T_{f},\Phi=0}$) plays the dominant role while $\tau_{T_{p}}$ and $\tau_{p}$ become minor, therefore the contributions from particles become less and less important as $Re$ increases. Hence, the effect of changing $\Phi$ on the normalized friction coefficients becomes smaller at higher $Re$ in the current parameter regimes.

On the other side, taking the particle-laden flow as a continuous effective fluid, the interaction between the particles and the resulting particle distributions are related to the rheology of a particle-laden flow, which can be quantified by the effective viscosity, $\nu_{eff}$. In the Stokes regime ($Re\ll 1$), the rheology of dense granular suspensions has been extensively investigated, and the dependence of the effective viscosity on the particle volume fraction has been discussed (Guazzelli & Pouliquen, 2018). In the current semi-dilute regime ($\Phi\leq 10\%$), the normalized effective viscosity of the granular flow in the Stokes regime, $\nu^{S}_{eff}$, can be approximated using the Krieger & Dougherty (K-D) formula (Krieger & Dougherty, 1959):

where $\nu_{f}$ the viscosity of the fluid, $\Phi_{m}$ the maximum packing particle fraction ($\Phi_{m}=0.65$ is used here as done in Stickel & Powell (2005)) and $[\eta]={5\over 2}$ for rigid spheres. In this work, the $\nu_{eff}$ was calculated using the same method as in previous works (Bakhuis et al., 2021; Yi et al., 2021), and the details can be found in Appendix B. As shown in figure 3, we compare the experimental results with the K-D formula. Though the model captures the increasing trend of the effective viscosity, the $Re$-dependence of the effective viscosity is missing, indicating that the relation in the Stokes regime is no longer valid in the current situation. The disparity found between the experimental data and the model can be understood since equation 9 was obtained based on the assumptions that the flow is inertialess and the particles distribute uniformly in the Stokes regime. However, for the current turbulent flow ($Re\sim 10^{4}$) and the suspended finite-size particles, the particle inertia is non-negligible and can be measured by the Stokes number, which is far beyond unity as shown in the later section. Moreover, as has been reported in previous numerical studies (Picano et al., 2015; Ardekani & Brandt, 2019), the particles show collective effects, which are also found in our experiments and will be discussed later. Note that when $Re=6.5\times 10^{3}$, the flow could be in the classical regime with laminar boundary layers and a turbulent bulk with Taylor vortex (Grossmann et al., 2016), which accounts for the more rapidly increasing trend than the model and other experimental data. Using the effective viscosity, the shear-thinning effect of the particle-laden flow is also examined, which can be well described by the Herschel–Bulkley model (Herschel & Bulkley, 1926), and the details can be found in Appendix B.

Next, we investigate the effects of changing the particle aspect ratio on the drag of the system. Particles with three kinds of aspect ratio are used: $\lambda=1/3$ (oblate), $\lambda=1$ (spherical) and $\lambda=3$ (prolate). The results are already given in figure 2 and figure 3. It is found that the particles could increase the drag of the system regardless of their aspect ratio, which is attributed mostly to the hydrodynamics interactions and rarely to the frictional contact between the particles in the current low $\Phi$ range (Guazzelli & Pouliquen, 2018). The shear-thinning effect and the increasing trend of the drag modification with increasing $\Phi$ are found to be ubiquitous in all cases.

It is remarkable, however, that particles with different aspect ratios could affect the drag of the system to different degrees, even for the same $Re$ or $\Phi$. As shown in figure 2, for a given $\Phi$, the $c_{f,\Phi}/c_{f,\Phi=0}$ of the oblate case decrease much more rapidly than that of the spherical and prolate cases, indicating that the suspended oblate particles result in a more pronounced shear-thinning effect, which can also be verified by the effective viscosity (see Appendix B). Moreover, for a given $Re$ and $\Phi$, the suspensions of the spherical, prolate and oblate particles result in the greatest, moderate and minimal drag, respectively. In the investigated ranges, the largest discrepancy of $c_{f,\Phi}/c_{f,\Phi=0}$ between suspensions of particles with different aspect ratios is $4\%$, which occurs at $\Phi=10\%$ and $Re=1.3\times 10^{4}$ in the spherical and oblate cases. Given the impressive accuracy of the drag measurements (less than 1%) and the relatively low $\Phi$, this disparity in the drag modification is quite noticeable. In addition, the difference of $c_{f,\Phi}/c_{f,\Phi=0}$ increases with increasing $\Phi$, which is due to the increasing importance of particle-induced stress at higher $\Phi$ (Ardekani & Brandt, 2019). Therefore, it is reasonable to speculate that, as $\Phi$ further increases, the differences of the drag modification between different aspect ratios would become larger. Note that, this result is different from the particle-laden channel flow in the previous study (Ardekani & Brandt, 2019), where they found that only the spherical particles increase the drag and the other two types of particles reduce the drag.

The dependence of the drag on the particle aspect ratio provides an experimental clue to the understanding of the mechanisms of bubbly drag reduction. Flow with dispersed bubbles can, under certain conditions, result in significant drag reduction (Van den Berg et al., 2005; van den Berg et al., 2007; Ezeta et al., 2019). However, bubbles are deformable, making it impossible to fix the bubble shape and therefore to isolate the effects of the bubble shape. van Gils et al. (2013) and Verschoof et al. (2016) have reported that the bubble deformation is crucial for achieving significant drag reduction, but the contributions of (i) the deformation process and (ii) the ultimate bubble shape after the deformation remain unknown. Obviously, the rigid particle suspensions are significantly different from bubbly flow, including the slip conditions at the surface, the resistance of the dispersed phase to straining motion, the polydispersity, the presence/absence of contaminants (surfactants), etc, and therefore a one-to-one comparison between these two types of flow is also unrealistic. However, the one that should be emphasized is that, as done in this work, the shape effect of rigid spheroidal finite-size particles on the overall drag provides a relatively effortless implement to study the shape effects of bubbles during the deformation process. Notwithstanding how simplistic the rigid particles are compared to bubbles, these particles (rigid particles and bubbles) behave somehow similar in such as the kinematic motion (Mathai et al., 2020), making the study of the rigid particles a reasonable choice to provide a reference case for bubbly flow studies. In this work, the rigid spheroidal particles are used to study their shape effect on drag. Due to the frictional contact between solid surfaces, rigid particles dispersed in the flows result in the drag enhancement (Guazzelli & Pouliquen, 2018), which is contrary to the drag reduction in the bubbly flow. However, the fact that the particle aspect ratio affects the drag modification in the current particle-laden flow, also hints that the bubble shape could affect the drag of the bubbly flow. Given that the shape effect on the overall drag is relatively small (compared to the drag modulation amplitude in bubbly flow, which is around $40$% at $\Phi=4\%$ (Verschoof et al., 2016)), one would expect that the bubbles can also modulate the turbulence through other mechanisms next to the change of their shapes. Therefore, the unique features of bubbly flow relative to rigid particle-laden flow are worthy of attention for future studies.

To reveal the physical mechanism of the drag modification, we look into the particle distribution by performing optical measurements. It should be noted that, even at the minimal particle volume fraction ($\Phi=2\%$), the particles cause violent light scattering and cannot be accurately detected by algorithms. Therefore, we start with a smaller volume fraction ($\Phi=0.5\%$, the number of particles $\simeq 1700$) to reduce the intensity of the light scattering, so that the particles in the images can be detected using the ellipse detection algorithms (Lu et al., 2020). One may question that whether the particle distribution obtained at this low volume fraction could qualitatively represent that at higher volume fractions since the particle dynamics change greatly with increasing $\Phi$. To assuage this misgiving, we perform torque measurements at $\Phi=0.5\%$ and the results have been shown in figures 2,3,8,9,10. It can be seen that the quantities related to the drag of the system at $\Phi=0.5\%$ (including $c_{f}$, $\nu_{eff}$ and $G$) show similar trends to that at higher volume fractions, hinting that the particles behave in similar ways at the low and high volume fractions in the current parameter regimes. Here, as depicted in figure 4(a,c), we define the region where the plumes are ejected from the inner (outer) boundary layer to the bulk as the ejected (injected) region, and the region between them is defined as the vortex-center region. The spherical particle positions which are obtained from different frames when $\Phi=0.5\%$ at $Re=6.5\times 10^{3}$ and $Re=1.3\times 10^{4}$ are given in figure 4(a) and figure 4(b), respectively. When $Re=6.5\times 10^{3}$, it is found that the particles distribute non-uniformly, and the pattern of the particle distributions reminisces about the Taylor vortices in TC flow (Grossmann et al., 2016). The particle distributions show a lower clustering near the inner wall in the ejected region, while in the injected region, the lower clustering emerges near the outer wall. However, when the $Re$ is increased to $1.3\times 10^{4}$, the particle distribute nearly homogeneously in the entire system (figure 4(b)).

To quantitatively evaluate the collective effects of the particles, we choose three typical regions of the flow structures (i.e., the injected region, ejected region, and the vortex-center region) and calculated the PDF of the particle radial positions. As shown in figure 5, the PDF is consistent with the pattern of the particle distributions in figure 4. When $Re=6.5\times 10^{3}$, the PDF has the minimal value near the inner wall in the ejected region, indicating the emergence of the lower clustering region. While in the injected region, the PDF peaks near the inner wall (i.e., highly clustering region) and approaches the lowest value near the outer wall. Not surprisingly, as shown in figure 5(b), the PDF curves tend to overlap with each other when the $Re$ increases, indicating that the distribution of the particles becomes less non-uniform in the entire system at higher $Re$. Note that, the inner wall peak appearing in the ejected region at higher $Re$ (figure 5b), which is opposite to the behavior found at low $Re$ (figure 5a), and this might result from the slight shift of flow structure positions since the positions of the Taylor vortex could be not exactly the same at low and high $Re$.

The difference of the particle distributions at different $Re$ can be interpreted from the perspective of the evolution of flow structures and the finite-size effect of the particles. On the one side, though the TC flow displays flow structure at high $Re$ (Huisman et al., 2014), the Taylor vortices become turbulent and the velocity fluctuations increase as $Re$ increasing (Grossmann et al., 2016), resulting in the more vigorous particle motion. On the other side, the finite-size effect of the particle can be characterized by the particle Reynolds number. Since the slip velocity of the particles is not accessible in the current experiments, the particle Reynolds number is estimated using the bulk velocity of the flow $u_{b}$, i.e., $Re_{p}=u_{b}d_{e}/\nu$, where $u_{b}$ can be estimated as $u_{b}\simeq 0.4\cdot\omega_{i}r_{i}$ in the present Reynolds number regime (Grossmann et al., 2016). Considering that the flow Reynolds number is defined by equation 1, one can obtain the particle Reynolds number as:

Hence, the $Re_{p}$ ranges from $2.6\times 10^{2}$ to $10^{3}$ in the experiments, and the finite-size effects of the particle become more pronounced at higher $Re$.

Further, the particle inertia can also be measured by the particle Stokes number, $St=\tau_{p}/\tau_{\eta}$, where $\tau_{p}$ and $\tau_{\eta}$ are the inertia response time of the particles and the Kolmogorov timescale of the flow, respectively. For spherical particles suspended in TC flow, the Stokes number is reduced to:

i.e., the $St$ (or, the particle inertia) increases with increasing $G$ and $Re$. For non-spherical particles, $St_{non}=c\cdot St$, where $c=f(\lambda)$ is a constant only depends on $\lambda$ (Voth & Soldati, 2017). Since the suspension is dilute here ($\Phi=0.5\%$), the non-dimensional torque $G$ can be approximated to that of the single-phase case, $G_{sp}$, giving the $St$ ranging from $10$ to $60$. Therefore, in the current $Re_{p}$ and $St$ ranges, the suspended particles could induce unsteady wakes and exhibit inertia clustering (Toschi & Bodenschatz, 2009; Mathai et al., 2020). Specifically, for the case of $Re=6.5\times 10^{3}$ in figure 4(a) (correspondingly, $Re_{p}=2.6\times 10^{2}$ and $St=10$), the particles are principally driven by the strong plumes released from the boundary layer, which are depicted by the arrows in figure 4(c). While for the case of $Re=1.3\times 10^{4}$ in figure 4(b) (correspondingly, $Re_{p}=5.2\times 10^{2}$ and $St=25$), the particles with the increased inertia can escape more easily from the flow structures and distribute more uniformly.

Since the particle distributions are inhomogeneous in the axial direction of the system, we use the averaged PDF, which is obtained by calculating the arithmetic average value of the three typical regions, to quantitatively represent the collective effects of particles in the entire system. As shown in figure 6(a-c), for a given $\lambda$, the trends of the averaged PDF curves remain when the $Re$ increases from $6.4\times 10^{3}$ to $Re=1.3\times 10^{4}$, suggesting that the increasing turbulent intensity has negligible effects on the averaged particle migration in the radial direction of the system in the current situation.

Nevertheless, for different $\lambda$, the averaged PDF curves show distinct differences regardless of the $Re$. In other words, particles with different $\lambda$ show different collective behaviors near the walls and thereby are expected to affect the boundary layer to different degrees. Indeed, the formation of particle layers has been reported in the particle-laden channel flow in the previous simulation works (Costa et al., 2016; Ardekani et al., 2017; Ardekani & Brandt, 2019). Ardekani & Brandt (2019) have shown that the spherical particles form a particle layer, which could enhance the near-wall Reynolds stress. While for the cases of the oblate and prolate particles, the near-wall Reynolds stress is attenuated due to the absence of the particle layer. In the present work in TC flow, as indicated in figure 5 and figure 6(b), the spherical particles preferentially cluster and form a particle layer near the inner wall, where the boundary layer exists. Therefore, the maximum drag modification is observed in the flow laden with spherical particles. However, the near-wall clustering phenomenon disappears in the oblate case (figure 6(a)), of which the PDF peaks in the bulk and has the lowest value near the walls. The averaged PDF curves indicate that the oblate particles prefer to cluster in the bulk and thereby have a smaller effect on the boundary layer, accounting for the observations of the minimal drag modification in figure 2 and figure 3.

However, one more question remains since so far the analyses of particle distribution are based on the optical measurements at $\Phi=0.5\%$: will the particle preferential clustering persist at even higher $\Phi$? To bridge the gap of volume fractions between the torque measurements ($\Phi\geq 2\%$) and optical measurements ($\Phi=0.5\%$), we further conduct optical measurements at higher volume fractions. Since the domain is illuminated from the back, the light intensity decreases greatly at higher volume fractions, making the particles hard to be distinguished from the background fluid. At higher volume fractions, the feasible method to detect the particles is manual detection, which is of relatively low accuracy but provides reliable information. The averaged PDF curves of particle distribution at $\Phi=2\%$ are shown in figure 6(d-f). For each case, 100 frames are used and the total number of detected particles is more than 3000. For even higher volume fractions ($\Phi\geq 4\%$), the enormous amount of unfocused particles make the particles inside the focal plane invisible, therefore no useful data can be obtained.

One can see in figure 6(d-f) that the spherical particles still preferentially cluster near the inner wall and therefore form a particle layer, which is consistent with the result obtained at $\Phi=0.5\%$. Additionally, at $\Phi=2\%$, the distribution of spherical particles becomes more non-uniform in the radial direction than that at $\Phi=0.5\%$, which has also been reported in Fornari et al. (2016), hinting that the spherical particles tend to move toward the walls as $\Phi$ increases. This growing non-uniform distribution possibly results from the stronger shear-induced particle-particle interactions at higher $\Phi$, and could partially account for the increasing drag differences compared with other particle shapes as $\Phi$ increases. Moreover, as $\Phi$ increases, the particles might relaminarize the bulk flow, and thus suppress the Reynolds stress (Fornari et al., 2016). This Reynolds stress suppression could therefore yield a stronger $Re$-dependence of the drag at higher $\Phi$, which is verified by the stronger shear-thinning effects of normalized friction coefficient (figure 2) and the normalized effective viscosity (figure 10). On the other hand, for oblate and prolate cases, it is clear that most of the particles distribute in the bulk region. The near-wall clustering phenomenon disappears and therefore the particle layer is absent, which is similar to that of $\Phi=0.5\%$ and accounts for their similar overall drag.

So far, the optical measurements performed at $\Phi=0.5\%$ and $\Phi=2\%$ lead to the same conclusion that the larger drag modification is connected to the near-wall particle clustering in the current system. Though the result and conclusion of particle preferential clustering might not be directly extrapolated further to higher volume fractions, one reasonable speculation would be that it will play a significant role in the drag modulation at higher $\Phi$. Of course, further studies are needed to confirm it in future work.