Scaling of shear-induced diffusion and clustering in a blood-like suspension

In dilute suspensions of hard spheres, the collision of two isolated spheres in the Stokes regime is symmetric and reversible. In the absence of long-range interactions diffusion arises from the collision of three or more particles [1,2] or due to surface roughness [3]. However, the collision of two or more elastic particles, such as red blood cells (RBCs), can alone give rise to diffusive behavior at long times [4–6]. Pair collisions of RBCs have been studied under a variety of conditions [7–9], as well as the shear-induced diffusion of RBCs in dilute systems [10], yet the mechanism behind the transport of deformable particles in dense suspensions has received less attention and is not yet fully understood. The hydrodynamic interactions of blood components like platelets, white blood cells and micro-/nano-particles with RBCs determine their distribution inside the blood vessels as well as the efficacy of drug-delivery systems [11,12] and are an active research topic [12–14].

Recently, Gross et al. [15] proposed a connection between viscosity and diffusion in such systems using kinetic arguments, obtaining good agreement with their simulation results. The range of validity for this relation has not been identified. In this letter we present results from simulations of a blood-like suspension in a shear-flow environment, focusing on the shear-induced diffusion from the perspective of the micro-structure and the collective effects, employing a wide range of shear rates and volume fractions.

We employ the combined immersed-boundary lattice-Boltzmann method (IB-LBM) in two dimensions [16]. Fluid is simulated with LBM using the D2Q9 LBGK scheme, while RBCs are represented as closed deformable biconcave membranes with a diameter $D=6.86\ \mu\text{m} \equiv2r$ and a thickness ${d=1.8\ \mu\text{m}}$. The membrane of an RBC consists of 26 neighboring vertices connected by Hookean springs and enhanced with bending resistance. In addition, forces that ensure the conservation of area (2D equivalent of volume) and that consider the cell-cell interactions are specified. This model has been used and validated in previous studies [16,17].

Simulations are performed in a Lees-Edwards periodic boundary conditions (LEbc) domain, providing an unbounded type of shear flow and avoiding wall-induced migration [18]. The size of the domain was set to $400\times400\ {\mu\text{m}^2}$ and the flow profile has the form $\{u_x, u_y\}=\{\dot{\gamma}(y-y_0), 0\}$ with $y_0=200\ \mu\text{m}$. The kinematic viscosity of the fluid is $\nu_0=1.7 \times 10^{-6}\ \text{m}^2/\text{s}$ and the density $\rho=1025\ \text{kg}/\text{m}^3$, close to that of blood plasma [19]. The LB relaxation time is set to $\tau=1$ and the lattice constant to $\Delta{x}=1\ \mu\text{m}$, yielding a timestep $\Delta t=9.8 \times 10^{-8}\ \text{s}$. Inner and outer viscosities and densities are equal, while a soft repulsive force is used to prevent overlap. The form of the repulsive force is $F_{\mathrm{rep}}= -C_\mathrm{rep}/h\textbf{e}_{ij}$, where h is the gap between surface points, $\textbf{e}_{ij}$ the connecting vector and $C_\mathrm{rep}$ a constant set to 0.05 in lattice units. The force becomes zero for distances larger than $r_\mathrm{cutoff} = 0.77 \ \mu\text{m}$. Similar power-law forces have been used by [20,21].

The shear rates studied range from $3\ \text{s}^{-1}$ to $1000\ \text{s}^{-1}$ and the volume fractions from $\phi=0.1$ to 0.5, yielding an absolute number of RBCs of 1600 to 8000, respectively. In all cases the simulation time was $\dot{\gamma}{t} > 10$ with a typical duration of $\dot{\gamma}{t} \propto 10^2$, sufficient to measure the transport properties of RBCs. A constant number of platelets is present in the system as part of a larger study regarding their transport in channel and shear flows, but due to their low relative numbers ($N_\mathrm{platelet}=320$ with $\phi_\mathrm{platelet}\ll1{\%}$) and small size, their contribution to viscosity and diffusion is insignificant (diffusion data not shown). The particle Reynolds number $Re_p = \frac{4 \dot\gamma r^2 }{\nu_0}$ ranges from $8.3 \times 10^{-5}$ to $2.8 \times 10^{-2}$ for $\dot{\gamma}=3\ \text{s}^{-1}$ to $1000\ \text{s}^{-1}$, respectively. Thermal fluctuations are neglected in our study, since the thermal Péclet number $Pe_t=\frac{\dot{\gamma} r^2}{D_t} \underset{\dot{\gamma}= 3\ \text{s}^{-1}}{=} 10^4 \gg 1$.

Performing a reliable analysis of the transport of RBCs requires sufficiently large domain sizes and shear strains [22], which leads to very long simulation times, especially in the lower shear rate regimes. It is common to rescale the timestep of the simulation according to the capillary number $Ca= \frac{\rho \nu_0 \dot{\gamma} r}{C_\mathrm{spr}}$, where $C_\mathrm{spr}$ is the elastic modulus of the cell. IB-LBM [15,23] and dissipative particle dynamics (DPD) [24,25] studies have followed this approach, motivated by the assumption that the shear elasticity provides the most dominant contribution to the suspension stress. Such scaling, however, deprives the timescales defined by the particle Reynolds number, lubrication, and interparticle forces from an appropriate scaling. Lubrication and interparticle forces dominate dense suspensions, where the gap between particles becomes comparable to the extents of these forces [15,26]. In contrast to these earlier studies we avoid such scaling because in our opinion this contributes to a consistent and well-defined model.

As with blood, this suspension exhibits a shear-thinning behavior: its viscosity decreases when subjected to shear strain. In all simulations RBCs were randomly positioned in the LEbc domain and the relative apparent viscosity was computed using Batchelor's method [27]. In low volume fractions ${\phi=0.0, 0.1, 0.2}$ the fluid is Newtonian¹ and the effective viscosity of the suspension $\nu_\mathrm{eff}= \frac{\nu_{\text{app}}-\nu_0}{\nu_0}$ remains relatively constant for all the examined shear rates (fig. 1 and inset). The suspension starts exhibiting non-Newtonian characteristics for $\phi=0.3$ and strong shear-thinning behavior in high volume fractions, $\phi=0.4,0.5$.

Zoom In Zoom Out Reset image size

As the mean particle distance decreases with increasing ϕ, RBCs interact more frequently and the effective volume of the RBCs appears to be practically larger for $\phi=0.4$ and 0.5, due to the extent of the repulsive force [28]. The shear thinning observed is an effect of both repulsion and deformability [26,29].

Shear-induced diffusion is the result of excluded-volume effects [1]. Neighboring particles are displaced from their initial streamlines when passing each other, thus leading to net post-collisional displacements, as opposed to the reversible two-particle collisions in the Stokes flow. Breedveld et al. [1] suggested that the diffusion tensor of a hard-sphere suspension should scale linearly with $\dot{\gamma}$ as

$$\begin{equation} \textbf{D}=\dot{\gamma} r^2 \mathbf{\hat{D}}(\phi), \end{equation} \tag{ 1 }$$

where $\mathbf{\hat{D}}(\phi)$ is a dimensionless symmetric tensor. The arguments are based on a dimensional analysis, in which the only relevant dimensions are the particle radius r as the lengthscale, and $\dot{\gamma}^{-1}$ as the timescale, considering the shear-induced character of the process. Equation (1) is widely used to describe the diffusion of RBCs and platelets with advection-diffusion models [30–32].

We measure diffusion in the shear-gradient direction D_y from the mean square displacement (MSD) as ${D_y=\underset{t\to\infty}{\lim}\frac{1}{2t} \langle\Delta {R_y(t)}^2\rangle}$, where $\Delta{R}_y$ is the displacement in the y-direction. It is well known that in non-Brownian suspensions MSD typically exhibits two distinct regimes: over short shear strains $\dot{\gamma} \Delta{t}$ the velocity of a particle is self-correlated leading to a super-diffusive quadratic behavior $\langle\Delta{R}^2_y(\Delta{t}) \rangle\propto (\Delta{t})^2$ (ballistic regime) and at large shear strains to the diffusive regime $\langle\Delta{R}^2_y(\Delta{t})\rangle\propto (\Delta{t})$, where the interactions of particles have lead to decorrelation [15]. This behavior is also observed in our model as shown in fig. 2. It is consequential that the duration of the simulations has to span a sufficient number of shear strains $\dot{\gamma}{t}$, until the ballistic regime is exceeded and MSD is well into the diffusive regime.

Zoom In Zoom Out Reset image size

As shown in fig. 3, $D_y/(\dot{\gamma} r^2)$ is increasing with ϕ, with diffusion following the linear scaling with $\dot{\gamma}$ in low volume fractions. Contrary to the predictions of (1) $D_y/(\dot{\gamma} r^2)$ scales with $\dot{\gamma}^{-0.5}$ for higher volume fractions.

Zoom In Zoom Out Reset image size

Gross et al. [15] made a similar observation for three-dimensional suspension of deformable particles in Couette flow. Using kinetic and scaling arguments, they established a connection between viscosity ν and diffusivity D in which $D/\dot{\gamma} \propto \mathrm{const}$ in the Newtonian regime, where q indicates the scaling of the viscosity with respect to the capillary number $(\nu\propto Ca^{-q})$, and $D/\dot{\gamma} \propto Ca^{-q/2}$ in the shear-thinning regime. While for the Newtonian cases this connection holds, for the shear-thinning regime our simulations do not confirm their findings, in which $D_y/\dot{\gamma} \propto \dot{\gamma} ^{-0.5}$ and not $\propto \dot{\gamma}^{-q/2\equiv-0.3}$. The difference in the dimensionality of the problem, in the form of missing shear elasticity and reduced mobility, may explain partly the differences in our results. However, the same arguments and assumptions hold in 2D and since Gross et al. do not employ relations or theorems that involve dimensionality, like the equipartition theorem, the explanation should be sought in their assumptions. One of the principal assumptions of the study was that the shear-elasticity provides the most dominant contribution to the suspension stress, rendering capillary number the only relevant quantity under examination. This assumption disregards the contribution of interparticle forces, which is substantial at high volume fractions and can, alone, induce shear thinning even in hard-sphere suspensions [26]. With proper comparison, the shift between our results and Gross et al. [15] could be explained in terms of differences in the resolution and the reach of the repulsive force, highlighting the role of the interparticle force in high volume fractions.

The deviation of diffusion from the linear scaling with the shear rate requires further investigation. Looking into the other parameters of eq. (1), only the radius r is possible to be dependent on $\dot{\gamma}$. Since RBCs are deformable, which stretch and tank-tread under shear [33], it is reasonable to associate stretch and angle of orientation of RBCs with a shear-dependent increase in the collisional cross-section, rendering $r\equiv r(\dot{\gamma})$. A decrease in the collisional cross-section results in smaller post-collisional migrations and eventually to smaller values of diffusion. However, our data do not support this assumption. Even though an increase of the average orientation angle 〈θ〉 is measured (fig. 4(a)), there is no notable difference with respect to the shear rate. The average maximum diameter of an RBC $D_{\max}$, defined as the maximum distance between any two surface points (SPs) of an RBC, shows that cells are stretched in all cases of high shear rates and in high volume fractions (fig. 4(b)). The increase in high shear rates is an effect of the shear forces and the deformability of the RBCs, while the latter can be attributed to the interparticle force and the decreasing mean gap with increasing volume fraction. The collisional cross-section $R\equiv r(\dot{\gamma})$ is computed by integrating $D_{\max}$ over all angles $R=\int_{0}^{2\pi} P(\theta)(d + D_{\max}\sin{(\theta)} ) \mathrm{d}\theta$, where $P(\theta)$ is the distribution of angles and d the RBC thickness. The collisional cross-section R is almost independent of $\dot{\gamma}$, slightly decreasing with $\log\dot{\gamma}$ for $\phi=0.4,0.5$, as the fittings of fig. 4(c) show. This change of $r(\dot{\gamma})$ is not sufficient to explain the non-linear scaling of diffusion in high volume fractions. These results demonstrate that pair collisions, which drive diffusion in dilute systems, are not the principal mechanism for dense systems and higher-order collisions and collective effects should be taken into account.

Zoom In Zoom Out Reset image size

The pair distribution function $g(\mathbf{r})$ (PDF) describes the microstructure from the perspective of the spatial arrangement of particles. It expresses the probability of finding particles at distance r from a reference particle normalized by the number concentration of the suspension [34] (fig. 5). In the dilute case, where particle interactions are less frequent, $g(x,y)$ has a uniform distribution lacking local maxima. Local maxima develop close to the small RBC radius for larger ϕ, while the distinguishable ripples indicate rouleux-like structures. The second peak close the large RBC diameter $r=7.5\ \mu\text{m}$ indicates that neighboring RBCs are often aligned in the flow. RBCs appear stretched at higher shear rates of $\phi=0.5$ (fig. 4(b)).

Zoom In Zoom Out Reset image size

The formation of RBC clusters is a key finding that is extracted from the pair distribution function. These clusters are different from the hydro-clusters of hard-sphere suspensions as studied by Brady and Bossis [35], owing to the repulsive interparticle force keeping RBCs separated and preserving a lubrication layer. We perform a cluster size analysis to obtain a quantitative estimate on collective effects of the suspension and the RBC-RBC interactions.

Two RBCs are considered "in contact" and belong to the same cluster when the distance between two of their surface points (SPs) is less than a cutoff distance $r^c_{\text{cutoff}}$, defined as the peak of the SP-SP distance distribution. The collection of RBCs that are in contact forms a cluster and the number of particles in one cluster defines the size s of the cluster. Figure 6(a) shows a snapshot of the simulation, coloring the RBCs belonging to the same cluster. The total number of clusters with size s at time t is denoted as $n(s,t)$. The average number of clusters n(s) is obtained over a time period T, and is defined as $n(s) = \frac{1}{T} \sum _{t=t_0}^{t_0 +T}n(s,t)$ [36], with $T=7\dot{\gamma}^{-1}$ being several times larger than the correlation time $\tau_c$ extracted from the velocity autocorrelation function and long enough to capture the duration of RBC contacts. The quantity we are interested in is the cluster size distribution $p_n(s)=n(s)/ \sum_{s_i=1}^{\infty}n(s_i)$.

Zoom In Zoom Out Reset image size

Figure 6(b) shows the cluster size distribution $p_n(s)$. Cluster sizes s > 1 are not frequent for volume fraction $\phi=0.1$ (dilute case), signifying that pair collisions are more frequent than three- (or multi-) particle collisions as they are expressed by $p_n(s>2)$, being 10 times less frequent than $p_n(s=2)$. The maximum cluster size $s_{\max}$ remains well under 10 for all dilute cases.

For $\phi=0.3$ and 0.5 $p_n(s>2)$ is higher indicating multi-particle collisions and collective motion. The extreme case of $\phi=0.5$ indicates that one big cluster dominates the flow. With the mean gap being comparable to $2\cdot r_{\text{cutoff}}$ though, it is highly unlikely that it is a finite-size effect which will change by increasing the domain size.

The percentage of "free" RBCs $p_n(s=1)$ is decreasing with ϕ and only slightly depends of the shear rate, as shown in fig. 6(a). For the same ϕ, higher $\dot{\gamma}$ result in lower single-cell percentages, revealing higher-order collisions. The average cluster size $\langle s\rangle$ is a way to measure the order of RBC collisions. Figure 6(d) shows that $\langle s\rangle$ remains close to 1 for $\phi=0.1$ and 0.2, a slight increase with $\dot{\gamma}$ is observed for $\phi=0.3$, while $\phi=0.4$ is increasing with $\dot{\gamma}$, hinting on many-particle collisions. Case $\phi=0.5$ is dominated by a hyper-cluster resulting in very large values for $\langle s\rangle$.

An attempt to fit the data to the kinetic clustering model (KCM) using non-linear least squares and obtain a typical cluster size s₀ is shown in fig. 6(d). KCM was developed by Raiskinmäki et al. [37] using kinetic arguments to study the basic mechanisms of shear thickening in particulate suspensions. The cluster size distribution of KCM follows a Poissonian, $n(s) \propto s^{-3/2} \exp{(-s/s_0)}$, with s₀ denoting a typical cluster size and was also used by Ding and Aidun [36] to uncover a universal scaling relation for the cluster size distribution, independently of particle shape and concentration. The typical cluster size is increasing with ϕ, until the fitting algorithm fails, due to the appearance of big clusters in $p_n(s)$ in the form of tails. For $\phi=0.1$, the typical cluster size is below 1 RBC, in which due the sporadic interaction of RBCs s₀ is underestimated.

The duration of a contact t_c is defined as the maximum time interval in which the distance of any of the surface points between a pair of RBCs, remains under $r_{\text{cutoff}}$. If two RBCs escape this distance, t_s is stored and a new measurement starts. The mean duration $\langle t_c\rangle$ of a contact is increasing with ϕ and $\dot{\gamma}$ (fig. 6(e)), showing that in these cases clusters live longer. The duration distribution $p_c(t_c)$ (fig. 6(f)) appears to follow an exponential distribution describing the time between events in a Poisson process, which seems to be in agreement with the KCM approach. Since the suspension is governed mostly by two-particle collisions for $\phi=0.1$, the extent of $p_c(t_c)$ indicates the duration of a collision in the dilute limit, where approximately after $\dot{\gamma}{t_c}=2$, none of the pair contacts has survived. For larger hematocrits, however, the distribution is more widespread, indicating that more than two-particle collisions are present and sustain a "bond" for larger times.

In this letter we have investigated the transport of RBCs in a shear-flow environment. The shear-induced diffusion of the suspension exhibits a departure from the linear scaling with the shear rate, resulting in a less efficient mixing. This result is important due to the wide use of the linear scaling model in coarse-grained models using the advection diffusion equation [30–32]. The analysis provided by Gross et al. [15] with similar findings for a similar system, is not confirmed by our results. The decrease of the normalized diffusivity $(D_y/(\dot{\gamma} r^2))$ due to a potential decrease in collisional cross-section, is by itself not enough to explain this departure, indicating that pair collisions are not the driving mechanism behind transport in dense suspensions as collective effects are taking place. A cluster analysis shows that multi-particle collisions resulting in collective effects dominate higher-hematocrit simulations. For future work a parametric study covering various interparticle forces and elastic moduli would be enlightening regarding the mixing of RBCs.

This research received funding from the European Union, 7th Framework Programme under grant agreement No. 26996 (THROMBUS). Author AGH acknowledges partial financial support by the Russian Scientific Foundation, grant No. 14-11-00826. We thank Profs. P. M. A. Sloot, G. E. Karniadakis and D. Bonn for fruitful discussions.

Additional remark: Nine movies (hmt10yd1e1.mp4, hmt10yd1e2.mp4, hmt10yd1e3.mp4, hmt30yd1e1.mp4, hmt30yd1e2.mp4, hmt30yd1e3.mp4, hmt50yd1e1.mp4, hmt50yd1e2.mp4, hmt50yd1e3.mp4) from our simulations visualizing the motion of RBCs were added as supplementary material.