Influence of each cylinder’s contribution on the total effective viscosity of a two-dimensional suspension by a two-way coupling scheme

Einstein’s viscosity formula is sometimes strongly limited for viscosity estimation of suspensions; that is, it is only applicable for low-concentration suspensions in which hydrodynamic interactions are sufficiently negligible. In particular, hydrodynamic interactions between particles (cylinders in two dimensions) should be taken into consideration when finite-size particles are suspended. Therefore, change in the microstructure, i.e., spatial arrangement of particles in the flow field, is important for understanding mechanism of suspension rheology. In order to provide better practical applications for viscosity estimation instead of Einstein’s formula, we investigated the influence of each cylinder’s contribution on the total effective viscosity of a suspension with finite-size cylinders considering the microstructure, especially in terms of cylinder-wall and cylinder-cylinder distances. Two-dimensional pressure-driven flow simulations were performed using the regularized lattice Boltzmann method and a two-way coupling scheme. The rigid circular cylinders suspended in a Newtonian fluid were assumed to be neutrally buoyant and non-Brownian. As a result, we found that both distances between cylinders and cylinder-wall are significant for viscosity estimation. In addition, the effective viscosity can be estimated accurately when the confinement is sufficiently low ( C ≈ 0.04). It can be stated that the microstructure of the suspension is one of the promising factors to estimate and control suspension rheology.


Introduction
It is important to understand the rheological properties of a suspension since its flows are ubiquitous in our everyday life, such as chemical production, food processing, pharmaceutics, mining engineering and biological systems (blood) (Stickel and Powell, 2005). When it comes to discussing the rheological properties of a suspension, effective viscosity has been recognized as an essential factor. According to Einstein's viscosity formula (Einstein, 1906), the effective viscosity ηeff of a suspension can be assessed by the volume fraction φ of the suspended particles as ηeff = η0[1 + [η]φ], where η0 is the solvent viscosity (i.e., the viscosity of a particle-free fluid) and [η] is the intrinsic viscosity. It is well known that the intrinsic viscosity is to be [η] = 2 for two-dimensional (2D) systems (Brady, 1984) and [η] = 2.5 for threedimensional (3D) systems (Einstein, 1906). When applying Einstein's viscosity formula, there are some assumptions to satisfy the state where hydrodynamic interactions are sufficiently negligible, i.e., low concentration (dilute suspension), small suspended particles (low confinement), and steady and uniform flow (homogeneous flow). For example, finite-size of the suspended particle has a significant effect on the particle dispersion. In other words, hydrodynamic interactions between particles should be taken into consideration when finite-size particles are suspended. In fact, Mueller et al. Okamura, Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) (2010) mentioned that the most important factors for rheology of a particle suspension are various types of bulk flow fields (e.g., shear thinning and shear thickening) and particles (in terms of shape, deformability, interactions and spatial arrangement), as well as the particle concentration.
It has been suggested that change in the microstructure-which is generally defined by the positions of the suspended particles relative to each other-is important for understanding mechanism of suspension rheology (Stickel and Powell, 2005;Morris, 2009). Inertial migration of particles is a typical factor for changing the microstructure. Silberberg (1961, 1962) first observed that a single rigid sphere in pipe flow migrated to an equilibrium position with its center located around 0.6R, with R being the pipe radius. The phenomenon of radial migration driven by inertia was termed as the "Segré-Silberberg effect" or the "tubular pinch effect". According to this inertial migration, many studies have been performed from experimental (Matas et al., 2004a(Matas et al., , 2004b, numerical (Feng et al., 1994;Yang et al., 2005) and theoretical (Asmolov, 1999;Schonberg and Hinch, 1989) approaches. Matas et al. (2004aMatas et al. ( , 2004b experimentally examined the migration of dilute suspensions of neutrally buoyant particles in Poiseuille flow. The equilibrium positions obtained in the experiments were found to migrate toward the channel wall as the Reynolds number Re increases. By contrast, they shift toward the center of the channel when the confinement, which is the ratio of the particle size to the channel width, is increased. Di Calro (2009), Martel and Toner (2014) reviewed that the inertial migration of particles in microfluidics and summarized the primary two key factors which control the phenomenon. The first one is the shear gradient lift force, which arises due to the asymmetry in the relative velocity between the fluid and the particle on either side of the particle. The other prominent factor is the wall repulsion force, which is due to a pressure increase between the particle and the wall.
Regarding the relationship between the effective viscosity and the microstructure of the suspension, Doyeux et al. (2016) numerically examined the role of particle-wall and particle-particle interactions in a 2D linear shear flow in the presence of walls. They reported that the effective viscosity increases when the particle approaches the wall. Fukui et al. (2020) evaluated the inertial effects of the suspended particles on temporal changes in the microstructure of a dilute suspension, in which particles were randomly suspended, for the different Reynolds numbers: Re = 4; Re = 16; Re = 256. For the case Re = 256, particles flowed on the centerline of the channel due to inertial effects, and the thixotropic behavior of the suspension was observed. This phenomenon is usually observed in microcirculation known as "Fåhraeus-Lindqvist effect" (Fåhraeus and Lindqvist, 1931), in which red blood cells flow on the centerline and the apparent effective viscosity decreases. These results suggest that the radial position of the particles in the channel is one of the most important factors in terms of assessing the rheological properties of a suspension. Consequently, Einstein's equation cannot be applied to such non-homogeneous (heterogeneous) flows in which particles are locally clustering, for example, in the center of the channel or near the wall surface. Doyeux et al. (2016) suggested that the total effective viscosity of suspensions can be simply expressed as a summation of each particle's contribution; however, there seems to exist little understanding of this suggestion. It should be added that few studies have been quantitatively evaluated by focusing on the effects of particle-wall and particleparticle distances on the rheological properties of a suspension. In spite of these circumstances, application of this idea is promising, for example, it may enable us to directly estimate the local effective viscosity without using rheometer from snapshots of suspensions obtained in experimental studies such as by Kawaguchi et al. (2019) and Nordstrom et al. (2010). Therefore, it is important to investigate the influence of each particle's contribution on the total effective viscosity in detail. This research is also expected to expand the range of application from microfluidics to nanofluidics in the future. However, the properties of nanofluids are much more complicated because they depend on many parameters, and theoretical formula for predicting nanofluids viscosity seems to be still limited and not practical at this moment (Mishra et al., 2014).
When the microstructure of a particle suspension is numerically considered, 2D simulations are useful in terms of greatly reduced computational time. They have been used to study multiparticle systems such as suspensions of red blood cells (Thiébaud et al., 2014), vesicles (Ghigliotti et al., 2010) or capsules (Li and Ma, 2010). Since our research performs 2D simulations, hereinafter, it will be referred to as "cylinder" instead of "particle". Regarding the governing equation, the lattice Boltzmann method (LBM) has been widely applied in the simulation of solid-liquid suspensions because it is ideally suitable for length scales related to microfluidics, i.e., flow fields with low Reynolds numbers (O'Connor et al., 2016). Because of this, we employed a two-way coupling scheme and the regularized lattice Boltzmann method (explained later in detail) for solid-liquid suspensions with finite-size cylinders at Reynolds number Re = 4. The purpose of this research is, therefore, to conduct 2D pressure-driven flow simulations and investigate the influence of each Okamura, Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) cylinder's contribution on the total effective viscosity of a suspension considering the microstructure, which is regarded as cylinder-wall and cylinder-cylinder distances in this paper.

Computational model
We performed pressure-driven flow simulations in a two-dimensional (2D) channel using a two-way coupling scheme-in which fluid and cylinders influence each other-so as to investigate the influence of each cylinder's contribution on the total effective viscosity of a suspension, as shown in Fig. 1 (results presented in Figs. 5-9 at Sec. 3.2 and Sec. 3.4). The computational domain is Lx × Ly = D × 2l (0 < x < D, −l < y < l) and 800 × 800 grid points (at least 32 grid points per cylinder diameter) were employed. Rigid circular cylinders with radius r suspended in a Newtonian fluid are neutrally buoyant (i.e., the density of the cylinder is the same as that of the fluid) and non-Brownian. The cylinders were initially placed at regular intervals in the region of 0 ≤ y < l, considering the flow symmetry. It should be also noted that the cylinders flowed steadily, and no collisions or contacts between cylinders were included. The center of the cylinder's coordinate is expressed as (xc, k, yc, k), where k represents each cylinder number. The confinement C is defined as the ratio of the cylinder diameter 2r to the width of the channel 2l (Eq. (1)), and gx denotes the regular intervals between cylinders in x-direction (Eq. (2)).
where Nm represents the total number of suspended cylinders. The important points of the computational model to note are that the cylinders were initially arranged regularly and only allowed to move in x-direction with rotational motion, i.e., fixed in y-direction regardless of lift force (Udono and Sakai, 2017). This setting was designed to assess the effects of cylinder-wall distance yc / l and cylinder-cylinder distance gx / 2r on the viscosity estimation without consideration of inertial migration when cylinder inertia was dominant.
The simulation parameters for this study were chosen as follows: the confinement C = 0.04, 0.08, 0.16, the area fraction φ ≤ 12.57%, and the Reynolds number Re = 4.

Relative viscosity estimation
In this study, we evaluated the relative viscosity ηr of the suspension in order to investigate the effect of each Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) cylinder's contribution on its macroscopic rheology. The relative viscosity was determined from the change in the pressure drop (Fukui et al., 2018) as follows: where Δp and Δp0 represent the pressure difference over the channel length D resulting from the suspension and the equivalent to the cylinder-free fluid, respectively. Here, we propose ηestimated as an alternative index when estimating the viscosity in a different way in case if Einstein's formula cannot be applied. ηestimated is defined as a summation of each cylinder's contribution ηk as follows, if ηk can be obtained in advance: When the cylinders are fixed in y-direction, η k (C, yc) (k = 1, Nm) are all equal under the same confinement C and the cylinder position yc because ηk is a function of only C and yc as mentioned by Doyeux et al. (2016). If Eq. (4) could be properly employed, the effective viscosity of a suspension might be directly estimated without using rheometer from the information of confinement C and cylinder position yc obtained in a snapshot of the flow field.

Governing equations 2.3.1 Regularized lattice Boltzmann method
We used the regularized lattice Boltzmann method (Izham et al., 2011;Morinishi and Fukui, 2016) for the 2D 9velocity (D2Q9) model for the fluid analysis in suspension flow. This method is based on the lattice Boltzmann method (Chen and Doolen, 1998;Ladd, 1994), and beneficial in reducing memory usage and improving computational stability.
The distribution function fα with a discrete velocity eα for the lattice Boltzmann equation can be expressed by considering only the advection and collision of the particles. Approximating this distribution function up to second-order moment provides the recovery of the incompressible Navier-Stokes equation, which can be expressed by a0, bi and cij are the parameters which satisfy the following relationships ∑ e αi f α ∑ e αi e αj f α where ρ is the fluid density, ρui is the momentum, and Π ij neq is the non-equilibrium part of the stress tensor. The equilibrium distribution function f α eq is expressed by approximating the Maxwell equilibrium distribution function to the quadratic term, as shown in Eq. (9). The non-equilibrium distribution function therefore can be written as Eq. (10).
Okamura, Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) Accordingly, the time evolution equation in the regularized lattice Boltzmann equation is then where τ is the relaxation time, and 0.8 was applied for all our computations.

Cylinder movements
In this study, interactions between fluid and cylinders should be considered because effect of finite-size of the cylinder on the viscosity estimation is investigated. Therefore, we used a two-way coupling scheme (Fukui et al., 2018), considering not only the translational motion of the cylinders but also the rotational motion. Accordingly, cylinders' motions (in two dimension) are simply expressed by Newton's second law for translation and by the angular equation of motion for the rotation as follows, respectively: where Fc is the external force vector, m is the cylinder mass, xc is the position vector, Tc is the torque, I is the moment of inertia, and θc is the cylinder angle. The external force vector Fc = (Fcx, Fcy) acting on the cylinders are calculated from the hydrodynamic forces (drag force and lift force). These hydrodynamic forces are determined by the pressure and the shear stress on the cylinder surface obtained by performing fluid analysis (LBM). Both Eqs. (12) and (13) were solved by third-order Adams-Bashforth method and Adams-Moulton method, therefore the velocity vector uc = (ucx, ucy) of the suspended cylinders are given by where ρ is the cylinder density. The fixation of the y-direction motion of the cylinders (Udono and Sakai, 2017) at Sections 3.2 and 3.4 is realized by setting ucy = 0, where ucy is the y-direction component of the velocity vector.

Virtual flux method
We used the virtual flux method (Morinishi and Fukui, 2012;Tanno et al., 2006) to describe shape of the cylinders in a Cartesian grid. This method enables us to express the pressure field accurately around the cylinder with its simple algorithm. Figure 2 shows how to arrange the virtual boundaries on the Cartesian grid. In particular, for the lattice Boltzmann method using the D2Q9 model, it is necessary to place virtual boundary points at the intersections of the discrete velocities in eight directions and the cylinder surface.
First, we describe the boundary conditions for obtaining the physical quantity qvb ((uvb, vvb), velocity; pvb, pressure) on the virtual boundary. In order to satisfy the no-slip boundary condition on body surfaces, velocities on virtual boundary points are provided as (uvb, vvb) = (uc, wall, vc, wall), where (uc, wall, vc, wall) are the velocities on the wall surface of the cylinder. Pressures on virtual boundary points pvb are calculated under the Neumann boundary condition (∂p/∂n|vb = 0), which is an approximate pressure condition on the body surfaces.
Briefly, as shown in Fig. 2, we consider the case where the distribution function at point C (i+1, j) moves to point B (i, j) across the virtual boundary so as to find the distribution function at point B in the next time step. The distribution function at point C, however, cannot move to point B because the fluid is separated from the cylinder by the virtual boundary. Thus, we obtain the virtual distribution function at point C from the distribution function at the virtual boundary Okamura, Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021)  point "vb". The equilibrium distribution function f α (t, vb) at the point "vb" is calculated, as shown in Eq. (15), from the distribution function f α (t, B) and the equilibrium distribution function f α eq (t, B).
Lastly, the virtual distribution function f α * (t, C) and the virtual non-equilibrium distribution function f α eq* (t, C) at point C can be extrapolated using internal division ratios a and b.
Consequently, "Fluid" and "Cylinder", as illustrated in Fig. 2, can be computed independently with the cylinder surface as the virtual boundary.

Verification & Validation
To test the influence of the grid resolution, we carried out preliminary computations of a suspension, in which cylinders are homogeneously suspended, and assessed the results in terms of relative viscosity for confinement C = 0.04, area fraction φ = 12.57% and Reynolds number Re = 4, as shown in Fig. 3. It should be noted that the smallest suspended cylinders were chosen to verify our computations, and they were allowed to move freely in a pressure-driven flow, that is, move both in xand y-directions with rotational motion. We employed 200 × 200, 400 × 400, 800 × 800 and 1200 × 1200 lattice grids for the computational domain D × 2l, which corresponds to 8, 16, 32 and 48 grid points per cylinder diameter 2r, respectively. Figure 3(b) shows the relative viscosity of the suspension versus grid resolution. The solid line indicates that from Einstein's viscosity formula. The difference between the case with 32 grid numbers per cylinder diameter and that with 48 grid numbers is less than 0.1%. Note that the numerical relative viscosity slightly exceeds that of Einstein's formula partly due to high concentration condition. Thus, it can be stated that the grid resolution of 32 grid numbers per cylinder diameter provides good accuracy to consider viscosity estimation in this study, and the number of lattice grids for the computational domain was chosen to be corresponding 800 × 800.
To validate our numerical approach for viscosity calculation using the pressure drop represented by Eq.
(3), we compared our numerical results to Einstein's viscosity formula for confinement C = 0.04, cylinder concentration φ ≤ 12.57% and Reynolds number Re = 4, as shown in Fig. 4. As a result, we obtained the intrinsic viscosity [η] = 2.029, and the difference to that from Einstein's equation ([η] = 2.0) was 1.45%. Thus, it was confirmed that the numerical approach for viscosity calculation in this study was physically validated. The data are shown in Table 1 at Sec. 3.3 in detail. Okamura, Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) Normalized axial velocity (a) Normalized axial velocity distributions.

Cylinder
(b) Relative viscosity versus grid resolution.

Viscosity calculation: Single cylinder contribution
In this section, we confirmed the effects of confinement C and cylinder y-axis position yc on the relative viscosity. The viscosity obtained here corresponds to the contribution of the individual cylinder ηk in Eq. (4). Figure 5 shows the relative viscosity ηr depending on the normalized y-axis position y / l for several confinements C = 0.04, 0.08, 0.16. The computational grids per cylinder diameter 2r are 32 (C = 0.04), 64 (C = 0.08) and 128 (C = 0.16). In this case, the flow includes a single cylinder fixed in y-direction in the computational area (see Fig. 1). Note that the center of the channel and the channel wall correspond to y / l = 0.0 and y / l = 1.0, respectively. It can be seen that the viscosity of the system increases when the cylinder approaches the channel wall; moreover, the curve becomes steeper when the cylinder is closer to the wall, especially for high confinement C ≥ 0.08. It should be noted that the relative viscosity of a suspension strongly depends on the suspended cylinder's y-axis position, even if the concentration is the same, i.e., the relative viscosity depends on not only the concentration but also the suspended cylinder's y-axis position. Consequently, Einstein's viscosity formula is only applicable to homogeneously dispersed suspension with low confinement condition. The result here is consistent with the previously reported result (Doyeux et al., 2016). Judging from this finding, the cylinder-wall interaction, i.e., the cylinder-wall distance, could be one of the most important factors in assessing the effect of the microstructure of a suspension on the relative viscosity.

Viscosity estimation: Homogeneous flow
In the present study, we estimate the total effective viscosity of a suspension. Here, the relative viscosity was estimated for suspensions in which cylinders are homogeneously (uniformly) dispersed (i.e., cylinder-wall and cylindercylinder interactions are uniform) for cylinder concentration 0.50% ≤ φ ≤ 12.57%, in which Einstein's viscosity formula is applicable for validation. Table 1 shows the viscosity ηnumerical obtained from numerical analysis (results presented in Fig. 4) and the estimated viscosity ηestimated obtained from Eq. (4) using the data in Sec. 3.2. We represent the cylinder position closest to the channel wall as yc * / l. The reason we describe yc * / l is that the cylinder-wall interaction plays an important role in increasing the effective viscosity, as shown in Fig. 5. Error 1 and Error 2 were calculated as the error of ηnumerical and ηestimated with respect to Einstein's viscosity ηEinstein, respectively: where Error 1 is for validation of our numerical simulation, and Error 2 is for assessment of our proposed estimation using Eq. (4).  Okamura, Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) It is worth noting that Error 2 was less than 1.0% for the cylinder concentration φ ≤ 2.01%; on the other hand, it increased significantly for the cylinder concentration φ ≥ 4.52%. It could be argued that the cylinders approaching the channel wall (see yc * / l in Table 1) and interactions between cylinders have an important effect as the concentration increases. However, few studies have evaluated the viscosity by focusing on such microstructures of the suspension. In contrast, the viscosity estimation proposed in this paper is expected to provide the relative viscosity that takes account of the microstructure of the suspension.
Regarding the relative viscosity of a suspension in homogeneous flow, the viscosity obtained from our numerical analysis is more precise than that estimated from the proposed method in this paper for the concentration φ ≤ 12.57%, as shown in Table 1. Because of this, in the following Section 3.4, we evaluate the applicability of the viscosity estimation proposed in this paper for a suspension, in which cylinders are non-homogeneously dispersed, in terms of the error E calculated by Accordingly, it can be stated that the effective viscosity can be estimated using Eq. (4) when the error E is sufficiently small.

Viscosity estimation: Non-homogeneous flow
As a first step for viscosity estimation, we investigated the effects of cylinder-wall and cylinder-cylinder distances without considering the inertial migration of cylinders. Figure 6 shows the snapshots of the cylinder positions for a representative case with confinement C = 0.08 and cylinder y-axis position yc = 2l / 4. Note that no collisions or contacts between cylinders were considered. Similarly, cylinders flowed under the conditions of the cylinders' position yc = 0-3l / 4 and the cylinder distance gx / 2r = 0.25-5.25 (as an initial state) for different confinements C = 0.04, 0.08, 0.16. As mentioned above in Sec 2.1, the cylinders' movements were fixed in y-direction.  Figure 7 shows the relationship between the error E in viscosity estimation and cylinder distance gx / 2r in x-direction for different confinements C and cylinder y-axis positions yc (center, yc → 0; wall, yc → l). The dotted lines at 1.0% indicate the threshold in error evaluation defined in this study. The error was less than 1.0% at the cylinder distance gx / 2r > 2.0 regardless of any confinements (C = 0.04-0.16) and any cylinder positions in y-direction (yc = 0-3l / 4). On the other hand, the error increases significantly above 1.0% at the cylinder distance gx / 2r ≤ 1.08 when the cylinders approach the wall, especially for yc ≥ 2l / 4. From these results, it was confirmed that the cylinder-wall and cylinder-cylinder distance as well as the cylinder concentration play an important role in viscosity evaluation.
To understand viscosity estimation more clearly for each confinement, three-dimensional error distribution for a representative case with confinement C = 0.08 is depicted in Fig. 8. The distribution is divided into red and blue by the threshold of 1.0%. Deep blue indicates suitable conditions for proposed estimation using Eq. (4). It can be seen that the error E [%] in viscosity estimation increases when the cylinder distance gx / 2r decreases and/or cylinders' position yc / l   , Fukui, Kawaguchi and Morinishi, Journal of Fluid Science and Technology, Vol.16, No.3 (2021) homogeneously dispersed, for various confinements and Reynolds numbers. They reported that the relative viscosity of the suspension decreases with an increase in the particle size even when the concentration of the suspension remains the same. Moreover, they mentioned that the confinement should be less than 0.04 when applying Einstein's viscosity formula. Based on these results , it is considered that the hydrodynamic interactions affected by the confinement are important in viscosity estimation. Figure 9 shows two-dimensionally visualized error distribution (2D visualization of the results of Fig. 8 for simplicity) for different confinements C = 0.04, 0.08, 0.16. The error in viscosity estimation increases as the confinement increases; that is, the applicable range for Eq. (4) becomes more limited for high confinement C ≥ 0.08. It can be stated that the viscosity estimation method proposed here also requires the assumption that the suspended particles are sufficiently small, as in the case of Einstein's viscosity formula.
(a) C = 0.04 (b) C = 0.08 (c) C = 0.16 Fig. 9 Two-dimensionally visualized error distribution in viscosity estimation for the effects of cylinder distance gx / 2r and cylinders' position yc / l for different confinements C = 0.04, 008, 0.16. The distribution is divided into red and blue by the threshold of 1.0%.
In two dimensions, simple visualizations of the entire flow provide clear insight into three-dimensional (3D) problems that are sometimes difficult to understand. In addition, two-dimensional (2D) simulations are important tools for understanding hydrodynamic interactions. In a previous study, it is reported that 2D simulations are capable of capturing several essential physical properties similar to the 3D ones (Ghigliotti et al., 2010), except for the vacillating breathing motion (Danker et al., 2007;Misbah, 2006). Therefore, it can be stated that results in this paper are based on 2D simulations, but are expected to capture some aspects of 3D problems.
Regarding quantitative evaluation of the effective viscosity of a suspension, the method proposed in this paper may be capable of directly estimating its viscosity from the geometrical position data of cylinders and is expected to explain experimental image data without using rheometer. Therefore, we will examine the case in which cylinders (in two dimension) are randomly suspended in more detail for better practical applications for future works.

Conclusions
We conducted two-dimensional pressure-driven suspension flow simulations using a two-way coupling scheme so as to investigate the influence of each cylinder's contribution on the total effective viscosity. As a result, we found that both distances between cylinders and cylinder-wall are significant for viscosity estimation. In addition, the effective viscosity can be estimated accurately when the confinement is sufficiently low (C ≈ 0.04). The viscosity estimation method proposed in this paper is still fundamental, however, is promising to estimate the viscosity because the microstructure of a suspension is also taken into consideration. In our future works, we will investigate the applicability of confinement in viscosity estimation; furthermore, we will estimate the total effective viscosity of a complex suspension in which cylinders are randomly suspended for various Reynolds number conditions.