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

The transport of cells and substances in dense suspensions like blood heavily depends on the microstructure and the dynamics arising from their interactions with red blood cells (RBCs). Computer simulations are used to probe into the detailed transport-related characteristics of a blood-like suspension, for a wide range of volume fractions and shear rates. The shear-induced diffusion of RBCs does not follow the established linear scaling with shear rate for higher volume fractions. The properties directly related to RBC deformability —stretching and flow orientation— are not sufficient to explain this departure according to the model of Breedveld, pointing to the dominance of collective effects in the suspension. A cluster size analysis confirms that collective effects dominate high volume fractions, as the mean cluster size is above 2 and the number of “free RBCs” is significantly decreased in denser suspensions. The mean duration of RBC contacts in clusters is increased in the high volume fraction and shear rate cases, showing that these clusters live longer.

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.
Methods. -We employ the combined immersedboundary 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 μm ≡ 2r and a thickness d = 1.8 μ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 × 400 μm 2 and the flow profile has the form {u x , u y } = {γ(y − y 0 ), 0} with y 0 = 200 μm. The kinematic viscosity of the fluid is ν 0 = 1.7 × 10 −6 m 2 /s and the density ρ = 1025 kg/m 3 , close to that of blood plasma [19]. The LB relaxation time is set to τ = 1 and the lattice constant to Δx = 1 μm, yielding a timestep Δt = 9.8 × 10 −8 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 rep = −C rep /he ij , where h is the gap between surface points, e ij the connecting vector and C rep a constant set to 0.05 in lattice units. The force becomes zero for distances larger than r cutoff = 0.77 μm. Similar power-law forces have been used by [20,21].
The shear rates studied range from 3 s −1 to 1000 s −1 and the volume fractions from φ = 0.1 to 0.5, yielding an absolute number of RBCs of 1600 to 8000, respectively. In all cases the simulation time wasγt > 10 with a typical duration ofγt ∝ 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 platelet = 320 with φ platelet 1%) and small size, their contribution to viscosity and diffusion is insignificant (diffusion data not shown). The particle Reynolds number Re p = 4γr 2 ν0 ranges from 8.3 × 10 −5 to 2.8×10 −2 forγ = 3 s −1 to 1000 s −1 , respectively. Thermal fluctuations are neglected in our study, since the thermal Péclet number P e t =γ r 2 Dt = γ=3 s −1 10 4 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 = ρν0γr Cspr , where C 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. Results and discussion. -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 φ = 0.0, 0.1, 0.2 the fluid is Newtonian 1 and the effective viscosity of the suspension ν eff = νapp−ν0 ν0 remains relatively constant for all the examined shear rates ( fig. 1 and inset). The suspension starts exhibiting non-Newtonian characteristics for φ = 0.3 and strong shear-thinning behavior in high volume fractions, φ = 0.4, 0.5.
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 φ = 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 hardsphere suspension should scale linearly withγ as whereD(φ) 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γ −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][31][32]. We measure diffusion in the shear-gradient direction D y from the mean square displacement (MSD) as It is well known that in non-Brownian suspensions MSD typically exhibits two distinct regimes: over short shear strainsγΔt the velocity of a particle is self-correlated leading to a super-diffusive quadratic behavior ΔR 2 y (Δt) ∝ (Δt) 2 (ballistic regime) and at large shear strains to the diffusive regime ΔR 2 y (Δt) ∝ (Δ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γt, until the ballistic regime is exceeded and MSD is well into the diffusive regime.
As shown in fig. 3, D y /(γr 2 ) is increasing with φ, with diffusion following the linear scaling withγ in low volume fractions. Contrary to the predictions of (1) D y /(γr 2 ) scales withγ −0.5 for higher volume fractions.
Gross et al. [15] made a similar observation for threedimensional suspension of deformable particles in Couette flow. Using kinetic and scaling arguments, they established a connection between viscosity ν and diffusivity D in which D/γ ∝ const in the Newtonian regime, where q indicates the scaling of the viscosity with respect to the capillary number (ν ∝ Ca −q ), and D/γ ∝ 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 /γ ∝γ −0.5 and not ∝γ −q/2≡−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γ. 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 crosssection, rendering r ≡ r(γ). 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 ≡ r(γ) is computed by integrating D max over all angles R = 2π 0 P (θ)(d + D max sin (θ))dθ, where P (θ) is the distribution of angles and d the RBC thickness. The collisional cross-section R is almost independent ofγ, slightly decreasing with logγ for φ = 0.4, 0.5, as the fittings of fig. 4(c) show. This change of r(γ) 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.
The pair distribution function g(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 μm indicates that neighboring RBCs are often aligned in the flow. RBCs appear stretched at higher shear rates of φ = 0.5 ( fig. 4(b)).
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  g(x, y), as calculated from the distances between the centers of the particles. Lighter gray corresponds to higher probability, units of distance are in μm and the color scale is per subplot. RBCs are deformable with a biconcave shape and can have neighboring particles at distances shorter than their larger diameter, which explains the blurred RBC shape, contrary to hard spheres where the spherical shape is clearly visible in PDF plots. An asymmetry and a preferred orientation with increasing shear rate and volume fraction can be distinguished.
surface points (SPs) is less than a cutoff distance r c 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) = 1 T t0+T t=t0 n(s, t) [36], with T = 7γ −1 s · pn(s) with respect to shear rate (main) and typical cluster size s0 (inset) obtained from fitting cluster distribution to pn(s) = Cs 3/2 e −s/s 0 . A hyper-cluster dominates φ = 0.5, with s φ=0.5 and being many times larger than s φ=0. 4 , is not visible within the range plotted. (e) Mean duration tc of a contact between two RBCs and (f) distribution of contact durations between two RBCs, along with their respective fits to an exponential distribution. Symbol notation is shown in (c) and applies to all plots except for the main figure of panel (d).
being several times larger than the correlation time τ 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)/ ∞ si=1 n(s i ). Figure 6(b) shows the cluster size distribution p n (s). Cluster sizes s > 1 are not frequent for volume fraction φ = 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 φ = 0.3 and 0.5 p n (s > 2) is higher indicating multiparticle collisions and collective motion. The extreme case of φ = 0.5 indicates that one big cluster dominates the flow. With the mean gap being comparable to 2 · r 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γ result in lower single-cell percentages, revealing higher-order collisions. The average cluster size s is a way to measure the order of RBC collisions. Figure 6(d) shows that s remains close to 1 for φ = 0.1 and 0.2, a slight increase withγ is observed for φ = 0.3, while φ = 0.4 is increasing withγ, hinting on many-particle collisions. Case φ = 0.5 is dominated by a hyper-cluster resulting in very large values for s .
An attempt to fit the data to the kinetic clustering model (KCM) using non-linear least squares and obtain a typical cluster size s 0 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) ∝ s −3/2 exp (−s/s 0 ), with s 0 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 φ = 0.1, the typical cluster size is below 1 RBC, in which due the sporadic interaction of RBCs s 0 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 cutoff . If two RBCs escape this distance, t s is stored and a new measurement starts. The mean duration t c of a contact is increasing with φ andγ ( 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 φ = 0.1, the extent of p c (t c ) indicates the duration of a collision in the dilute limit, where approximately afterγ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.
Conclusions. -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][31][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 /(γ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 higherhematocrit 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.