1 Introduction
The rapid advancement of nanotechnology and nanomaterial manufacturing has led to emerging exposure of human biological systems, such as cardiovascular systems, to nanosized particulate matter (Albanese, Tang & Chan Reference Albanese, Tang and Chan2012; Malysheva, Lombi & Voelcker Reference Malysheva, Lombi and Voelcker2015), ranging from engineered nanocarriers functioning as medicine/imaging agents (Yoo et al. Reference Yoo, Irvine, Discher and Mitragotri2011; Blanco, Shen & Ferrari Reference Blanco, Shen and Ferrari2015; Griffin et al. Reference Griffin, Zhu, Liu, Aidun and Ku2018) to aerosol pollutant particulates potentially causing fatal cardiovascular disease (Brook et al. Reference Brook, Rajagopalan, Pope, Brook, Bhatnagar, Diez-Roux, Holguin, Hong, Luepker and Mittleman2010; Newby et al. Reference Newby, Mannucci, Tell, Baccarelli, Brook, Donaldson, Forastiere, Franchini, Franco and Graham2015; Miller et al. Reference Miller, Raftis, Langrish, McLean, Samutrtai, Connell, Wilson, Vesey, Fokkens and Boere2017). Both the design of optimal nanocarrier systems and the prevention and control of nanoparticle (NP) toxicity rely heavily on the knowledge of NP distribution in vascular systems (Albanese et al. Reference Albanese, Tang and Chan2012). Nevertheless, the biodistribution of NPs remains elusive to date and limits the bioavailability of NP systems to the human bio-environment. One of the primary bottlenecks is the lack of understanding on the NP dispersive mechanisms in blood flow with respect to specific haemorheological conditions.
Multiscale computational models (Lee et al. Reference Lee, Choi, Kopacz, Yun, Liu and Decuzzi2013; Muller, Fedosov & Gompper Reference Muller, Fedosov and Gompper2014; Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018) have been developed and applied to predict the NP distribution in microvessels by directly simulating NPs and red blood cells (RBCs) suspended in blood plasma. Although these models, as particle-based approaches, provide a straightforward means to predict the NP distribution in realistic, microscale vessels, they are still computationally intractable when it comes to organ- and circulation-level applications. Alternatively, continuum models (Eckstein & Belgacem Reference Eckstein and Belgacem1991; Decuzzi et al. Reference Decuzzi, Godin, Tanaka, Lee, Chiappini, Liu and Ferrari2010; Hossain et al. Reference Hossain, Zhang, Liang, Hussain, Ferrari, Hughes and Decuzzi2013; Mehrabadi, Ku & Aidun Reference Mehrabadi, Ku and Aidun2015) have the ability to predict NP distribution in large-scale vascular systems by solving the three-dimensional (3D) convection–diffusion equation
where
$t$
is time,
$\mathit{c}$
is the NP concentration,
$\boldsymbol{u}$
is the local fluid velocity and
$\boldsymbol{J}$
is the flux of NP concentration. Here,
$\boldsymbol{J}$
is often estimated by Fick’s law as
$\boldsymbol{J}=-\unicode[STIX]{x1D63F}^{\infty }\boldsymbol{\cdot }\unicode[STIX]{x1D735}c$
, where the NP diffusion tensor,
$\unicode[STIX]{x1D63F}^{\infty }$
, has so far been treated as isotropic, Brownian diffusivity (Hossain et al.
Reference Hossain, Zhang, Liang, Hussain, Ferrari, Hughes and Decuzzi2013) or solute diffusivity measured in a single principal direction (Zydney & Colton Reference Zydney and Colton1988). Since the particle diffusion tensor is anisotropic in nature even for monodisperse rigid-sphere suspensions (Foss & Brady Reference Foss and Brady1999, Reference Foss and Brady2000b
), an improved constitutive relation capturing the anisotropy of
$\unicode[STIX]{x1D63F}^{\infty }$
subject to local haemorheological properties is therefore necessary to form a better closure of this convection–diffusion problem.
Using Couette-type flow devices, experiments have been conducted to characterize the particle self-diffusivity (Eckstein, Bailey & Shapiro Reference Eckstein, Bailey and Shapiro1977; Breedveld et al. Reference Breedveld, Van den Ende, Tripathi and Acrivos1998) and effective solute diffusivity (Wang & Keller Reference Wang and Keller1985; Zydney & Colton Reference Zydney and Colton1988; Breedveld et al. Reference Breedveld, Van den Ende, Tripathi and Acrivos1998) in non-colloidal particle suspensions under linear shear flow for various shear rates and particle volume fractions. However, due to the difficulties of particle tracking in the presence of the affine flow effect, only the particle diffusivity in the cross-stream directions are reported. Apart from the experimental efforts, particle-scale simulations have become an important tool for characterizing the anisotropic particle diffusivity tensor in both colloidal (Foss & Brady Reference Foss and Brady1999, Reference Foss and Brady2000a ,Reference Foss and Brady b ) and non-colloidal (Sierou & Brady Reference Sierou and Brady2004; Yeo & Maxey Reference Yeo and Maxey2010; Clausen, Reasor & Aidun Reference Clausen, Reasor and Aidun2011) suspensions under shear flow. Owing to the success of those particle-scale simulation techniques, substantial progress has been made in understanding the rheological and hydrodynamic responses of the particle diffusion tensor in sheared monodisperse suspensions. Nonetheless, the bidisperse RBC–NP suspension system remains largely unexplored and entails unique transport phenomenology that is unavailable to the conventional monodisperse particle suspensions.
First, there is a large discrepancy in length scale between NPs (
${\sim}O(10~\text{nm})$
) and RBCs (
${\sim}O(10~\unicode[STIX]{x03BC}\text{m})$
). Consequently, NPs are subject to both molecular-level thermal fluctuations (Brownian motion) and cellular-level interactions with RBCs. On macroscopic scales, the two effects synergistically give rise to an apparent diffusivity contributed by both Brownian diffusivity (BD) and the so-called RBC-enhanced shear-induced diffusivity (RESID) (Mehrabadi et al.
Reference Mehrabadi, Ku and Aidun2015; Mehrabadi, Ku & Aidun Reference Mehrabadi, Ku and Aidun2016; Liu et al.
Reference Liu, Zhu, Rao, Clausen and Aidun2018). Second, the NP phase shows infinite dilution while the RBC phase exhibits a range of physiological concentrations from
${\sim}10\,\%$
to
${\sim}40\,\%$
. Consequently, the BD shows insignificant dependence on shear rate and haematocrit, while the RESID is highly dependent on the haemorheological conditions (Mehrabadi et al.
Reference Mehrabadi, Ku and Aidun2016). Third, RBCs deform considerably under shear, from biconcave shape in equilibrium to large elongation and tank treading of the membrane under high shear (Gross, Kruger & Varnik Reference Gross, Kruger and Varnik2014). Such geometric asymmetry and morphological changes of RBCs could alter RESID substantially.
Therefore, the objective of the present work is to characterize the bulk diffusivity tensor of NP in sheared blood flow and interrogate the NP dispersive mechanism specific to a broad range of haematocrit and shear rate. Given the large length-scale discrepancy (three to four orders of magnitude) between NP and RBC, resolving both particle phases using direct numerical simulation (DNS) is computationally prohibitive. Therefore, a multiscale complex blood flow solver (Aidun & Clausen Reference Aidun and Clausen2010; Reasor, Clausen & Aidun Reference Reasor, Clausen and Aidun2012, Reference Reasor, Clausen and Aidun2013a ; Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019) is employed to treat NPs as effective Brownian particles while directly resolving the RBC phase. Such a multiscale treatment can substantially reduce the computational expense but still preserve the critical suspension physics at distinct scales. As will be shown in the following section (§ 3.3.1), good comparison between the simulation and experimental results can be obtained using this multiscale approach. Since confinement (wall) effects in general lead to spatial heterogeneity of the blood flow (Kumar & Graham Reference Kumar and Graham2012) that forbids the calculation of NP diffusivity related to specific haematocrit and shear rate, the Lees–Edwards boundary condition (LEbc) (Lees & Edwards Reference Lees and Edwards1972) is implemented (MacMeccan et al. Reference MacMeccan, Clausen, Neitzel and Aidun2009) to impose unconfined simple shear flow to obtain the NP bulk diffusive properties.
Note that blood flow typically occurs under confinement (e.g. blood flow in arteries) involving heterogeneous flow structures that could induce cell segregation (Kumar & Graham Reference Kumar and Graham2012; Ahmed et al. Reference Ahmed, Mehrabadi, Liu, Barabino and Aidun2018) and margination (Zhao & Shaqfeh Reference Zhao and Shaqfeh2011; Zhao, Shaqfeh & Narsimhan Reference Zhao, Shaqfeh and Narsimhan2012; Reasor et al. Reference Reasor, Mehrabadi, Ku and Aidun2013b ; Mehrabadi et al. Reference Mehrabadi, Ku and Aidun2015, Reference Mehrabadi, Ku and Aidun2016). Such phenomena are found to be a synergistic outcome of the RBC-induced diffusion in the RBC-laden region and the formation of a cell-free layer (CFL) near the wall (Mehrabadi et al. Reference Mehrabadi, Ku and Aidun2015, Reference Mehrabadi, Ku and Aidun2016). For NPs with negligible inertia and much smaller length scale compared to the CFL thickness (Zhao et al. Reference Zhao, Shaqfeh and Narsimhan2012; Mehrabadi et al. Reference Mehrabadi, Ku and Aidun2016), it is expected that the presence of the wall has an insignificant direct influence on the NP diffusion in the RBC-laden region. Besides, large-scale problems suitable for continuum modelling (such as blood flow through coronary arteries) typically feature much larger length and time scales compared to those considered in the current cellular-scale studies. It is therefore plausible to hypothesize that the long-time NP diffusion tensor evaluated in an unbounded simple shear flow should closely capture the NP diffusive behaviour subject to the same local haemorheological condition in a macroscale heterogeneous blood flow environment.
One unique contribution of this work is the development of a multiscale-simulation-informed empirical expression that links the anisotropic NP diffusion tensor to the local haemorheological conditions. Such 3D NP diffusive information in sheared blood is intractable to measure through either experiment or DNS simulation. On the application side, the developed NP diffusion tensor provides a constitutive relation that can lead to more effective continuum-level models to tackle large-scale NP biotransport problems. On the suspension rheology side, the diffusive phenomenology observed in such a biophysical, bidisperse RBC–NP suspension system could entail novel suspension physics that is unavailable to conventional suspension flows.
The remainder of this article is organized as follows. In § 2, the details of the computational methodology are presented. In § 3, we present results and perform numerical experiments. In § 3.1, the numerical problem is formulated with careful consideration of numerical resolution as well as physiological significance. In § 3.2, the transient dispersive behaviours of NPs are presented. As follows in § 3.3, the long-time diffusive behaviours of NPs subject to a wide range of shear rate and haematocrit are validated with available experimental data and interrogated with various haemorheological scaling behaviours. In § 3.4, the NP–RBC suspension microstructure is analysed to give mechanistic insights into the haemorheological scaling observations. In §§ 3.5 and 3.6, numerical experiments are conducted to shed light on the physical mechanisms governing the nonlinear shear-rate dependence of the NP diffusion tensor. In § 3.7, empirical correlations are proposed based on the haemorheological scalings observed in previous sections. In § 4, we conclude this work with some remarks.
2 Methodology
The numerical method for this study is through a 3D lattice-Boltzmann (LB) based multiscale complex blood flow solver that efficiently resolves both the dynamics and interactions of nanoscale particles and microscale capsules (Aidun, Lu & Ding Reference Aidun, Lu and Ding1998; Aidun & Qi Reference Aidun and Qi1998; Aidun & Clausen Reference Aidun and Clausen2010; Reasor et al. Reference Reasor, Clausen and Aidun2012; Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019), as demonstrated in figure 1. The LB method is a well-established numerical model for hydrodynamics and proves to be a highly scalable method for DNS of dense particulate suspensions (Aidun & Clausen Reference Aidun and Clausen2010; Clausen, Reasor & Aidun Reference Clausen, Reasor and Aidun2010). Modelling of the RBC dynamics and deformation is via a coarse-grained spectrin-link membrane method (Pivkin & Karniadakis Reference Pivkin and Karniadakis2008; Fedosov, Caswell & Karniadakis Reference Fedosov, Caswell and Karniadakis2010) coupled to the LB method (Reasor et al. Reference Reasor, Clausen and Aidun2012), which has been validated against experimental results (Reasor et al. Reference Reasor, Clausen and Aidun2012, Reference Reasor, Clausen and Aidun2013a ). The NP suspension dynamics are resolved via a two-way coupled lattice-Boltzmann Langevin-dynamics (LB-LD) method with both particle Brownian motion and long-range hydrodynamic interactions (HI) directly resolved and validated (Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019). The solver has been successfully applied to several studies of particle and biopolymer transport in cellular blood flow (Reasor et al. Reference Reasor, Clausen and Aidun2013a ,Reference Reasor, Mehrabadi, Ku and Aidun b ; Mehrabadi et al. Reference Mehrabadi, Ku and Aidun2016; Ahmed et al. Reference Ahmed, Mehrabadi, Liu, Barabino and Aidun2018; Griffin et al. Reference Griffin, Zhu, Liu, Aidun and Ku2018; Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018).
(a) The RBC–NP bidisperse suspension in a triply periodic domain at shear rate
$\dot{\unicode[STIX]{x1D6FE}}=1000~\text{s}^{-1}$
and haematocrit
$\unicode[STIX]{x1D719}=40\,\%$
. (b) The RBC membrane immersed in the LB lattice domain. The RBC membrane dynamics and deformation are captured by the course-grained spectrin-link method (Pivkin & Karniadakis Reference Pivkin and Karniadakis2008; Fedosov et al.
Reference Fedosov, Caswell and Karniadakis2010; Reasor et al.
Reference Reasor, Clausen and Aidun2012). (c) The NP suspension dynamics is captured by the LB-LD method (Liu et al.
Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019), which captures the NP Brownian motion and long-range HI in an off-lattice fashion.

2.1 Lattice-Boltzmann method
Simulation of the suspending fluid is based on the LB method developed by Aidun & Lu (Reference Aidun and Lu1995), Aidun et al. (Reference Aidun, Lu and Ding1998) and Aidun & Clausen (Reference Aidun and Clausen2010). The LB method solves the discretized Boltzmann transport equation in velocity space through the streaming–collision process. In streaming, the fictitious fluid particles propagate along discrete velocity vectors forming a lattice space. In collision, the fluid particles at each lattice site collide with each other, causing the relaxation of the particle distribution function (PDF) towards a local ‘Maxwellian’ equilibrium PDF. The collision term is linearized based on the single-relaxation-time Bhatnagar–Gross–Krook (BGK) operator (Bhatnagar, Gross & Krook Reference Bhatnagar, Gross and Krook1954). The temporal evolution of the particle distribution function is given as
where
$f_{i}$
is the fluid PDF,
$f_{i}^{(0)}$
is the equilibrium PDF,
$r$
is the lattice site,
$e_{i}$
is the discrete lattice velocity,
$t$
is time,
$\unicode[STIX]{x1D70F}$
is the single relaxation time and
$f_{i}^{S}$
is a forcing source term introduced to account for the discrete external force effect (He et al.
Reference He, Zou, Luo and Dembo1997). This method has a pseudo speed of sound,
$c_{s}=\unicode[STIX]{x0394}r/(\sqrt{3}\unicode[STIX]{x0394}t)$
, and a fluid kinematic viscosity,
$\unicode[STIX]{x1D708}=(\unicode[STIX]{x1D70F}-\unicode[STIX]{x0394}t/2)c_{s}^{2}$
, where
$\unicode[STIX]{x0394}t$
is the time step and
$\unicode[STIX]{x0394}r$
is the unit lattice distance. The positivity of
$\unicode[STIX]{x1D708}$
requires
$\unicode[STIX]{x1D70F}>\unicode[STIX]{x0394}t/2$
. In the LB method, time and space are typically normalized by
$\unicode[STIX]{x0394}t$
and
$\unicode[STIX]{x0394}r$
, respectively, such that
$\unicode[STIX]{x0394}t_{LB}=\unicode[STIX]{x0394}r_{LB}=1$
are employed to advance equation (2.1). In the near-incompressible limit (i.e. Mach number,
$Ma=u/c_{s}\ll 1$
), the LB equation recovers the Navier–Stokes equation (Junk & Yong Reference Junk and Yong2003) with the equilibrium PDF given in terms of local macroscopic variables as
where
$\unicode[STIX]{x1D714}_{i}$
denotes the set of lattice weights defined by the LB stencil in use. The macroscopic properties such as the fluid density
$\unicode[STIX]{x1D70C}$
, velocity
$\boldsymbol{u}$
and pressure
$p$
are obtained via moments of the equilibrium distribution functions as
respectively. Here,
$\unicode[STIX]{x1D644}$
is the identity tensor, and pressure can be related to density and the speed of sound through
$p=\unicode[STIX]{x1D70C}c_{s}^{2}$
. For the D3Q19 stencil adopted in the current study,
$Q$
is equal to 19. Along the rest, non-diagonal and diagonal lattice directions,
$\unicode[STIX]{x1D714}_{i}$
is equal to
$1/3$
,
$1/18$
and
$1/36$
, and
$|\boldsymbol{e}_{i}|$
is equal to 0,
$\unicode[STIX]{x0394}r/\unicode[STIX]{x0394}t$
and
$\sqrt{2}(\unicode[STIX]{x0394}r/\unicode[STIX]{x0394}t)$
, correspondingly.
2.2 Spectrin-link method
Modelling of the RBC membrane is through a coarse-grained spectrin-link (SL) method (Pivkin & Karniadakis Reference Pivkin and Karniadakis2008; Fedosov et al.
Reference Fedosov, Caswell and Karniadakis2010) coupled with the LB method, which has been extensively validated with experimental results and proved to be a successful tool to capture both single RBC deformation and dynamics (Reasor et al.
Reference Reasor, Clausen and Aidun2012) and rheology of RBC suspensions at physiological haematocrit (Reasor et al.
Reference Reasor, Clausen and Aidun2013a
). In the SL model, the RBC membrane is modelled as a triangulated network with a collection of vertices mimicking actin vertex coordinates, denoted by
$\{\boldsymbol{x}_{n},n\in 1,\ldots ,N\}$
. The Helmholtz free energy of the network system,
$E(\boldsymbol{x}_{n})$
, including in-plane, bending, volume and surface area energy components (Dao, Li & Suresh Reference Dao, Li and Suresh2006), is given by
Here, the in-plane energy,
$E_{IP}$
, characterizes the membrane shear modulus through a worm-like chain (WLC) potential (Bustamante, Bryant & Smith Reference Bustamante, Bryant and Smith2003) coupled with a hydrostatic component (Fedosov et al.
Reference Fedosov, Caswell and Karniadakis2010). The bending energy,
$E_{B}$
, specifies the membrane bending stiffness, which is essential in characterizing the equilibrium RBC biconcave morphology (Dao et al.
Reference Dao, Li and Suresh2006; Fedosov et al.
Reference Fedosov, Caswell and Karniadakis2010). The volumetric constraint energy,
$E_{\unicode[STIX]{x1D6FA}}$
, and the area constraint energy,
$E_{A}$
, preserve the RBC volume and area conservation, respectively, when subject to external forces.
The dynamics of each vertex advance according to Newton’s equations of motion,
where
$\boldsymbol{v}_{n}$
is the velocity of the vertex
$n$
, and
$M$
is the fictitious mass equal to the total mass of the cell divided by the number of vertices
$N$
. The number of vertices used to discretize the RBC membrane is
$N=613$
, which has been shown to yield adequate resolution to resolve hydrodynamic forces (MacMeccan et al.
Reference MacMeccan, Clausen, Neitzel and Aidun2009) and capture single RBC dynamics (Reasor et al.
Reference Reasor, Clausen and Aidun2012) and concentrated RBC suspension rheology (Reasor et al.
Reference Reasor, Clausen and Aidun2013a
) when coupled with the LB method. Force
$\boldsymbol{f}_{n}^{LB}$
specifies the forces on the vertex due to the fluid–solid coupling;
$\boldsymbol{f}_{n}^{CC}$
are the forces due to cell–cell interactions; and the forces due to the Helmholtz free energy based on the SL model are determined by
The SL method is solved by integrating equations (2.5) at each LB time step using a first-order-accurate forward Euler scheme that is consistent with the LB evolution equation to avoid excessive computational expense (Reasor et al. Reference Reasor, Clausen and Aidun2012; Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018).
2.3 Langevin dynamics approach
The nanoscale particle suspensions are resolved through a two-way coupled LB-LD method (Liu et al. Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019), which has been shown to correctly capture the theoretical Brownian diffusivity and long-range many-body HI. This approach treats suspended particles in Stokesian regimes as point particles, while the particles with Brownian effect are coupled to the non-fluctuating LB fluid in a two-way fashion through spatial extrapolation/interpolation schemes (Ahlrichs & Dünweg Reference Ahlrichs and Dünweg1999; Peskin Reference Peskin2002; Mynam, Sunthar & Ansumali Reference Mynam, Sunthar and Ansumali2011). The dynamics of the LD particles is described by the Langevin equation (LE),
where
$m$
is the mass of a single particle. The conservative force
$\boldsymbol{C}_{p}$
, specifying the interparticle and particle–surface interaction forces, is determined by calculating the directional derivatives of the total potential energy
$U_{total}$
as
where the details of
$U_{total}$
are discussed in § 2.5. The frictional force
$\boldsymbol{F}_{p}$
is assumed to be proportional to the relative velocity of the particle with respect to the local viscous fluid velocity (Ahlrichs & Dünweg Reference Ahlrichs and Dünweg1998, Reference Ahlrichs and Dünweg1999),
where
$\boldsymbol{u}_{p}$
denotes the particle velocity, and
$\boldsymbol{u}(\boldsymbol{r}_{p},t)$
is the interpolated LB fluid velocity at the centre of the particle. The friction coefficient
$\unicode[STIX]{x1D701}$
is determined by the Stokes’ drag law,
$\unicode[STIX]{x1D701}=3\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}d_{p}$
, where
$\unicode[STIX]{x1D707}$
is the dynamic viscosity of the suspending fluid. The stochastic force
$\boldsymbol{S}_{p}$
explicitly gives rise to the Brownian motion of the particle and satisfies the fluctuation–dissipation theorem (FDT) (Kubo Reference Kubo1966) by
where
$i,j\in \{x,y,z\}$
,
$\unicode[STIX]{x1D6FC}$
and
$\unicode[STIX]{x1D6FD}$
run through all the particle indices,
$\unicode[STIX]{x1D6FF}_{ij}$
and
$\unicode[STIX]{x1D6FF}_{\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6FD}}$
are Kronecker deltas,
$\unicode[STIX]{x1D6FF}(t-t^{\prime })$
is the Dirac delta function,
$k_{B}$
is the Boltzmann constant and
$T$
is the absolute temperature of the suspending fluid. The angle brackets denote the ensemble average over all the realizations of the random variables. Since only time scales equal to and greater than the Brownian diffusion time scale are of interest, this study solves the overdamped discretized LE as suggested in Liu et al. (Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019).
2.4 Fluid–solid coupling
The coupling between fluid and RBC is accomplished through the Aidun–Lu–Ding (ALD) fluid–solid interaction approach, the details of which are well documented in previous publications (Aidun et al.
Reference Aidun, Lu and Ding1998; Aidun & Clausen Reference Aidun and Clausen2010; Reasor et al.
Reference Reasor, Clausen and Aidun2012). To resolve the NP dynamics subject to the hydrodynamics and the long-range HI among NPs, the LD particle and the non-fluctuating LB fluid phase are coupled in a two-way fashion, as discussed and verified in previous studies (Mynam et al.
Reference Mynam, Sunthar and Ansumali2011; Liu et al.
Reference Liu, Zhu, Rao, Clausen and Aidun2018, Reference Liu, Zhu, Lechman, Rao, Clausen and Aidun2019). Specifically, the hydrodynamic force exerted on an NP,
$\boldsymbol{F}_{p}^{H}$
, is systematically decomposed into frictional and stochastic components as
where the fluid velocity at the particle site,
$\boldsymbol{u}(\boldsymbol{r}_{p},t)$
, is interpolated based on the surrounding LB velocities and applied to update the LD particle dynamics through equation (2.7). The weighting functions,
$w(\boldsymbol{r},\boldsymbol{r}_{p})$
, for interpolation are constructed using a trilinear scheme (Ahlrichs & Dünweg Reference Ahlrichs and Dünweg1998; Liu et al.
Reference Liu, Zhu, Rao, Clausen and Aidun2018). Since
$\boldsymbol{F}_{p}$
and
$\boldsymbol{S}_{p}$
both originate from the ‘collision’ between NPs and liquid molecules,
$\boldsymbol{F}_{p}^{H}$
(instead of
$\boldsymbol{F}_{p}$
) is assigned back to the fluid phase to satisfy momentum conservation. The same weighting function is then applied to constructing the local forcing source term as
which is adopted by equation (2.1) to update the local hydrodynamics. The coupled LB-LD method, similar to the external boundary force (EBF) method (Wu & Aidun Reference Wu and Aidun2010), modifies the conventional LB evolution equation into equation (2.1) by adding the forcing distribution function
$f_{i}^{S}(\boldsymbol{r},t)$
, which is shown to approximate the Navier–Stokes equation in the macroscopic scale (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002).
2.5 Contact modelling
The contact model for RBCs that specifies
$\boldsymbol{f}_{n}^{CC}$
is based on the subgrid contact functions originally formulated in Ding & Aidun (Reference Ding and Aidun2003) and later improved by MacMeccan et al. (Reference MacMeccan, Clausen, Neitzel and Aidun2009) and Clausen et al. (Reference Clausen, Reasor and Aidun2011). It prevents the RBC from overlapping when cell–cell membrane separation is below one LB lattice spacing. In this model, the lubrication term is replaced with an exponential contact function to avoid numerical instability driven by the singular nature of the lubrication hydrodynamics and the discrete nature of the interparticle separation calculation, as explained in detail in MacMeccan et al. (Reference MacMeccan, Clausen, Neitzel and Aidun2009) and Clausen & Aidun (Reference Clausen and Aidun2010). The rheological insensitivity of the contact model to model parameters has been discussed in detail by Clausen et al. (Reference Clausen, Reasor and Aidun2011). This contact model has also been previously applied in the characterization of the rheological properties of concentrated deformable capsule (Clausen et al.
Reference Clausen, Reasor and Aidun2011) and RBC (Reasor et al.
Reference Reasor, Clausen and Aidun2013a
) suspensions with satisfactory agreement with experimental measurements. Although the lubrication force is shown to play some role in particle self-diffusion in sheared monodisperse rigid-particle suspensions (Foss & Brady Reference Foss and Brady2000a
,Reference Foss and Brady
b
), we anticipate such an effect to be insignificant in the current bidisperse suspension system given that NP diffusion is dominated by the RBC–NP interactions.
The NP–RBC contact model that provides
$U_{total}$
is based on the Morse potential documented in Liu et al. (Reference Liu, Zhang, Wang and Liu2004), which is an empirical model for particle–particle interactions that can be calibrated to match experimental measurements (Neu & Meiselman Reference Neu and Meiselman2002; Liu et al.
Reference Liu, Zhang, Wang and Liu2004). Owing to the variety and somewhat lack of statistics for actual NP–RBC short-distance interaction, this study employs the measured cell–cell interaction potential (Neu & Meiselman Reference Neu and Meiselman2002) for the NP–RBC interactions. The potential parameters are specified according to Liu et al. (Reference Liu, Zhang, Wang and Liu2004) but with a cutoff distance selected to preserve only the repulsive effect. The detailed formulation of the Morse potential has been discussed in appendix A with sensitivity analysis performed showing that the computation of NP diffusivity is insensitive to the change of the model parameters up to 60 %. This suggests that the NP diffusion is largely driven by the hydrodynamic interaction rather than by the direct contact between NP and RBC membrane. The short-distance NP–NP interaction is neglected due to the extreme dilution of the NP concentration (
$\ll 1\,\%$
) considered.
This work does not attempt to model any adhesive forces between RBCs since above the shear rate of
$5~\text{s}^{-1}$
the aggregation of RBCs is not significant (Fedosov et al.
Reference Fedosov, Pan, Caswell, Gompper and Karniadakis2011). Although the adhesion or uptake of NPs to cells may be influential to the NP dispersive behaviour (Shang, Nienhaus & Nienhaus Reference Shang, Nienhaus and Nienhaus2014), it is however not within the scope of this study.
2.6 Lees–Edwards boundary condition
Since the primary focus of this study is on the particle bulk diffusive behaviour subject to no wall effect, simulations are performed in an unbounded, triply periodic cubic domain where a constant shear rate is imposed through the Lees–Edwards boundary condition (LEbc) (Lees & Edwards Reference Lees and Edwards1972). This method, originally developed for molecular dynamics simulations, was extended to the LB method by Wagner & Pagonabarraga (Reference Wagner and Pagonabarraga2002) and later applied to deformable suspensions on parallel computing architectures by Clausen et al. (Reference Clausen, Reasor and Aidun2011) and Reasor et al. (Reference Reasor, Clausen and Aidun2013a
). In addition to the operations associated with regular periodic boundary conditions, both the particle (NP and RBC) phase and the fluid phase undergo a shift in position and velocity according to the LEbc scheme as they cross the top (
$+y$
) or bottom (
$-y$
) boundary.
2.7 Characterization of the particle diffusion tensor
The presence of shearing flow imposes a convective effect on the particle suspension and complicates the characterization of particle diffusion tensor. The major difficulty lies in determination of the longitudinal diffusivity (
$\unicode[STIX]{x1D60B}_{xx}^{\infty }$
) and off-diagonal diffusivity (
$\unicode[STIX]{x1D60B}_{xy}^{\infty }$
), which require careful subtraction of the affine particle displacement. The diffusion tensor in sheared monodisperse colloidal suspensions has been successfully quantified by sampling the non-affine particle mean-square displacements (MSDs) (Morris & Brady Reference Morris and Brady1996; Foss & Brady Reference Foss and Brady1999, Reference Foss and Brady2000b
; Zia & Brady Reference Zia and Brady2010). Therefore, this study calculates the long-time NP diffusion tensor,
$\unicode[STIX]{x1D63F}^{\infty }$
, in the form of
$$\begin{eqnarray}\unicode[STIX]{x1D63F}^{\infty }=\left(\begin{array}{@{}ccc@{}}\unicode[STIX]{x1D60B}_{xx}^{\infty } & \unicode[STIX]{x1D60B}_{xy}^{\infty } & 0\\ \unicode[STIX]{x1D60B}_{yx}^{\infty } & \unicode[STIX]{x1D60B}_{yy}^{\infty } & 0\\ 0 & 0 & \unicode[STIX]{x1D60B}_{zz}^{\infty }\end{array}\right),\end{eqnarray}$$
where each non-zero diffusion component is calculated by
3 Results
In this section, we first formulate the simulations based on both physiological and numerical rationales. As follows, the transient growth of the NP mobility is discussed to understand the short-time response of the NP dispersive behaviour at various haematocrits and shear rates. Then, focus will be shifted to understanding the NP long-time diffusive behaviour under different haemorheological conditions with appropriate scaling, where the simulation results are also compared with available experimental data. To gain insight into the mechanisms governing the nonlinear haematocrit and shear-rate dependence of the NP diffusion tensor, we visualize the NP–RBC microstructure and carry out numerical experiments. Eventually, we construct empirical correlations for the complete NP diffusion tensor.
3.1 Problem formulation
Simulation data. Each case is simulated independently in an LB domain of size
$128\times 128\times 80$
at specified
$Pe$
,
$Ca_{G}$
and
$\unicode[STIX]{x1D719}$
. NP long-time diffusivities,
$\unicode[STIX]{x1D60B}_{xx}^{\infty }$
,
$\unicode[STIX]{x1D60B}_{yy}^{\infty }$
,
$\unicode[STIX]{x1D60B}_{zz}^{\infty }$
and
$\unicode[STIX]{x1D60B}_{xy}^{\infty }$
, normalized by the theoretical Brownian diffusivity,
$\unicode[STIX]{x1D60B}^{B}$
, are listed for all cases. The unlisted off-diagonal diffusivities,
$\unicode[STIX]{x1D60B}_{xz}^{\infty }$
and
$\unicode[STIX]{x1D60B}_{yz}^{\infty }$
, are found to be negligible. The measured diffusivities have a standard deviation less than 5
$\,\%$
. RBCs have an effective radius of
$a_{2}=2.9~\unicode[STIX]{x03BC}\text{m}$
. Brownian diffusivity is calculated by
$\unicode[STIX]{x1D60B}^{B}=k_{B}T/6\unicode[STIX]{x1D707}\unicode[STIX]{x03C0}a_{1}$
at temperature
$T=310$
K.

The apparent diffusivity of NPs in unbounded blood flow under simple shear is determined by NP radius
$a_{1}$
, shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
and haematocrit
$\unicode[STIX]{x1D719}$
, with NP concentration in the dilute regime. The relevant dimensionless parameters (besides
$\unicode[STIX]{x1D719}$
) primarily include the NP Péclet number,
expressing the ratio of shear-induced diffusion to Brownian diffusion, and the RBC capillary number,
quantifying the competition between the fluid viscous stress and the membrane elastic stress. Here,
$\unicode[STIX]{x1D707}$
is the dynamic viscosity of the suspending plasma,
$a_{2}$
is the effective radius of the RBC and
$G$
is the elastic shear modulus of the RBC membrane;
$\unicode[STIX]{x1D60B}^{B}$
is the Brownian diffusivity, which is determined by the Stokes–Einstein relation,
$\unicode[STIX]{x1D60B}^{B}=k_{B}T/6\unicode[STIX]{x03C0}\unicode[STIX]{x1D707}a_{1}$
, where
$k_{B}$
is the Boltzmann constant and
$T$
is the absolute temperature. The NP Péclet number quantifies the severity of the NP Brownian effect, while the RBC capillary number determines the deformability of the RBC capsule.
(a) The RBC–NP complex configuration at
$t\dot{\unicode[STIX]{x1D6FE}}=0$
. NPs and RBCs are randomly and uniformly mixed in an unbounded shear flow. RBCs are shown partially in the domain to visualize the NP phase. The athermal system is determined by shear rate
$\dot{\unicode[STIX]{x1D6FE}}$
, haematocrit
$\unicode[STIX]{x1D719}$
and NP radius
$a_{1}$
. (b–f) The RBC–NP complex configuration in the long-time scale regime under a shear rate of (b)
$\dot{\unicode[STIX]{x1D6FE}}=10$
, (c) 100, (d) 1000, (e) 2000 and (f)
$10\,000~\text{s}^{-1}$
with
$\unicode[STIX]{x1D719}=0.4$
.

To obtain appropriate scaling relations, we performed a large number of independent 3D simulations. Table 1 lists all the simulation parameters and the measured NP diffusivities. Wide ranges of shear rates (
$10\leqslant \dot{\unicode[STIX]{x1D6FE}}\leqslant 10\,000~\text{s}^{-1}$
) and haematocrits (
$0\leqslant \unicode[STIX]{x1D719}\leqslant 0.4$
) with physiological relevance (Lipowsky Reference Lipowsky2005; Popel & Johnson Reference Popel and Johnson2005) are covered. Cases with
$\dot{\unicode[STIX]{x1D6FE}}=10\,000~\text{s}^{-1}$
are to match certain vascular pathological conditions, e.g. the high shear induced thrombosis (Casa & Ku Reference Casa and Ku2017). For discussions in this section, NP size is set to
$2a_{1}=100~\text{nm}$
. RBCs are assumed to be in a healthy state with an effective radius
$a_{2}=2.9~\unicode[STIX]{x03BC}\text{m}$
and a membrane shear modulus
$G=0.0063~\text{dyne}~\text{cm}^{-1}$
. The absolute temperature is set to
$T=37\,^{\circ }\text{C}$
, at which the plasma has a viscosity
$\unicode[STIX]{x1D707}=1.2~\text{cP}$
and a density
$\unicode[STIX]{x1D70C}=1.0~\text{g}~\text{cm}^{-3}$
. The viscosity ratio of RBC cytoplasm to plasma is set to the physiological value
$\unicode[STIX]{x1D706}=5.0$
. The density of cytoplasm is set to that of the plasma. The corresponding
$Pe$
and
$Ca_{G}$
lie in the ranges of
$0.0066\leqslant Pe\leqslant 6.60$
and
$0.0055\leqslant Ca_{G}\leqslant 5.52$
, respectively. All simulations are formulated by matching the dimensionless group, i.e.
$\unicode[STIX]{x1D719}$
,
$Pe$
and
$Ca_{G}$
.
Simulations are initiated by imposing steady shear flow on the uniformly, randomly mixed NPs and RBCs at specific shear rate in an LEbc computational domain, as shown in figure 2(a). The domain has a dimension of
$128\times 128\times 80$
(
$42.7~\unicode[STIX]{x03BC}\text{m}\times 42.7~\unicode[STIX]{x03BC}\text{m}\times 26.7~\unicode[STIX]{x03BC}\text{m}$
) in longitudinal (
$x$
), velocity-gradient (
$y$
) and vorticity (
$z$
) directions, respectively. This LB domain size matches the highest resolution applied for the rheological characterization of cellular blood flow under shear by Reasor et al. (Reference Reasor, Clausen and Aidun2013a
). The selected LB grid resolution (
$300~\text{nm}$
per lattice unit) and the equilibrium RBC mesh size (1.5 lattice units per link length) have previously proven to be fine enough to capture both the single RBC dynamics (MacMeccan et al.
Reference MacMeccan, Clausen, Neitzel and Aidun2009; Reasor et al.
Reference Reasor, Clausen and Aidun2012) and the rheological properties of concentrated cellular blood flow (MacMeccan et al.
Reference MacMeccan, Clausen, Neitzel and Aidun2009; Reasor et al.
Reference Reasor, Clausen and Aidun2013a
).
To obtain converged long-time diffusivity, sufficient strains (
$t\dot{\unicode[STIX]{x1D6FE}}\sim 1000$
) and a large number of particles (5000 NPs and up to 208 RBCs) are employed for each simulation. The resolution of these simulations in terms of strain units and number of particles is at the high end compared to other numerical studies on particle diffusion in colloidal/non-colloidal suspensions (Foss & Brady Reference Foss and Brady1999, Reference Foss and Brady2000b
; Sierou & Brady Reference Sierou and Brady2004; Yeo & Maxey Reference Yeo and Maxey2010; Clausen et al.
Reference Clausen, Reasor and Aidun2011; Gross et al.
Reference Gross, Kruger and Varnik2014; Mountrakis, Lorenz & Hoekstra Reference Mountrakis, Lorenz and Hoekstra2016). As will be discussed in § 3.3.1, the selected resolution produces good agreement between the simulation results and the available experimental data.
Simulations are performed on the Intel Xeon Skylake nodes of the TACC (Texas Advanced Computing Center) Stampede-2 system where each node features 48 cores and a 2.1 GHz clock rate. For the case at
$\unicode[STIX]{x1D719}=0.4$
with 208 RBCs and 5000 NPs, each run takes
${\sim}168$
h on 32 cores (
${\sim}5376$
core hours) to accomplish 1000 strains (
$t\dot{\unicode[STIX]{x1D6FE}}$
). The total computational cost for the 35 independent cases listed in table 1 is approximately 140 000 core hours.
3.2 Temporal growth of NP transient mobility
To understand the short-time dispersive characteristics of the NP phase, we examine the transient behaviour of NP mobility by tracking the temporal growth of NP MSDs. Cross-stream MSDs,
$\langle y(t)y(t)\rangle /a_{1}^{2}$
and
$\langle z(t)z(t)\rangle /a_{1}^{2}$
, are sampled from the initial configuration (
$t\dot{\unicode[STIX]{x1D6FE}}=0$
) where NPs and RBCs are randomly and uniformly mixed, as shown in figure 2(a). Longitudinal and off-diagonal MSDs,
$\langle x(t)x(t)\rangle /a_{1}^{2}$
and
$\langle x(t)y(t)\rangle /a_{1}^{2}$
, are sampled starting from a long-time configuration (
$t\dot{\unicode[STIX]{x1D6FE}}=400$
) to avoid any extra convective effect caused by the transient elongation and reorientation of RBCs. Such transient effects are found to introduce an extra affine displacement leading to
${\sim}t^{3}$
growth of MSDs and hence jeopardize the measurement of the long-time diffusivity in the
$x$
-relevant directions. The necessity of eliminating the affine effect for calculating streamwise diffusivities is elaborated by Foss & Brady (Reference Foss and Brady2000b
) in the context of sheared colloidal suspensions.
Normalized MSDs plotted against relative shear strains in logarithmic scale at various haematocrits for different shear rates. Panels (a)–(d) are concerned with evolutions of
$\langle y(t)y(t)\rangle /a_{1}^{2}$
; (e) and (f) depict evolutions of
$\langle x(t)x(t)\rangle /a_{1}^{2}$
and
$|\langle x(t)y(t)\rangle |/a_{1}^{2}$
, respectively; (g) and (h) depict evolutions of
$\langle z(t)z(t)\rangle /a_{1}^{2}$
. The imposed shear rate is denoted on each plot. Inset of (b) shows the RBC–NP distribution at
$t\dot{\unicode[STIX]{x1D6FE}}\sim 1$
under
$\dot{\unicode[STIX]{x1D6FE}}=100~\text{s}^{-1}$
. Insets of (c) show the RBC–NP distribution at
$t\dot{\unicode[STIX]{x1D6FE}}\sim 1$
(upper) and
$t\dot{\unicode[STIX]{x1D6FE}}\sim 100$
(lower) under
$\dot{\unicode[STIX]{x1D6FE}}=2000~\text{s}^{-1}$
and
$\unicode[STIX]{x1D719}=0.4$
. An NP radius of
$a_{1}=50$
nm is used for normalization. Cross-stream MSDs (a–d) and (g,h) are sampled starting from
$t\dot{\unicode[STIX]{x1D6FE}}=0$
to capture the transient NP dispersive behaviour; streamwise and off-diagonal MSDs are sampled starting from
$t\dot{\unicode[STIX]{x1D6FE}}=400$
to avoid the initial extra affine displacement induced by RBC elongation and orientation.

Figure 3 plots the particle normalized MSD growth with respect to the relative strains sampled. At
$\unicode[STIX]{x1D719}=0$
, all diagonal MSDs grow according to the Stokes–Einstein (SE) prediction with no dependence on shear rate, while the off-diagonal MSDs produce zero value (not shown). This is consistent with the isotropic nature of NP Brownian motion in a dilute and unbounded solution, and also serves as a verification of the MSD calculations. At
$\unicode[STIX]{x1D719}>0$
, a deviation from the SE relation occurs in all diagonal MSDs as a result of NP–RBC interactions. At
$\dot{\unicode[STIX]{x1D6FE}}=10~\text{s}^{-1}$
, only a slight deviation from the SE relation is observed given the dominance of Brownian diffusion; see figure 3(a).
The transient deviation of NP mobility from the SE relation can be further interrogated through the evolution of cross-stream MSDs, as depicted in figure 3(a–d, g–h). In the short-time regime (
$t\dot{\unicode[STIX]{x1D6FE}}\ll 1$
), a linear growth of MSD (
${\sim}t$
) is observed particularly under low-shear-rate (e.g.
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 100~\text{s}^{-1}$
) or high-shear-rate and high-hematocrit (e.g.
$\dot{\unicode[STIX]{x1D6FE}}\geqslant 2000~\text{s}^{-1}$
and
$\unicode[STIX]{x1D719}=0.4$
) conditions, as shown in figure 3(b) and (c), respectively. The short-time linear growth of MSD suggests a short-time diffusive mechanism. Since the short-time linear growth of MSD always occurs before the ballistic regime (
${\sim}t^{2}$
) where RBC–NP collisions start to occur, such a short-time diffusive driver should logically be the Brownian effect. Therefore, the initial linear MSD growth at low shear rates (with
$\unicode[STIX]{x1D719}=0\sim 0.4$
) or high shear rates (with
$\unicode[STIX]{x1D719}=0.4$
) can be explained by the Brownian diffusive time scale being much shorter than the RBC–NP collision time scale under such haemorheological conditions. However, under high shear rate (e.g.
$\dot{\unicode[STIX]{x1D6FE}}\geqslant 2000~\text{s}^{-1}$
), the diffusive behaviour (
${\sim}t$
) turns to ballistic (
${\sim}t^{2}$
) as
$\unicode[STIX]{x1D719}$
increases to 0.4; see figure 3(c,d,g,h). This is probably due to the high convective (high-shear-rate) and low inertial (low-haematocrit) effects that reduce the RBC–NP collision time scale to be comparable to the Brownian diffusion time scale. As a result, the Brownian diffusion of NP is overwhelmed by the ballistic behaviour caused by insufficient RBC–NP collisions.
RBC–NP partial pair distribution function,
$g_{21}(\boldsymbol{r})$
, projected on the
$xy$
plane under various haemorheological conditions and averaged in the strain range of
$t\dot{\unicode[STIX]{x1D6FE}}=1{-}3$
.

In the intermediate-time regime,
$t\dot{\unicode[STIX]{x1D6FE}}=O(1){-}O(10)$
, anomalous dispersive behaviour is observed as RBCs start to elongate and rotate to be aligned with the streamwise direction due to shear. For cases at high shear rates (
$\dot{\unicode[STIX]{x1D6FE}}\geqslant 2000~\text{s}^{-1}$
) and high haematocrit (
$\unicode[STIX]{x1D719}=0.4$
), we observe a sublinear growth of MSD, representing a temporary hindrance of the NP mobility. This hindrance effect is caused by a string-ordered microstructure of NP in the shearing plane accompanied by the elongation and alignment of RBCs; see the contrast between the upper inset of figure 3(c) and the inset of figure 3(b) (see also supplementary movie 2, available at https://doi.org/10.1017/jfm.2019.320). Such string-ordered NP distribution can be better visualized by plotting the RBC–NP partial pair distribution function (PPDF),
$g_{21}(\boldsymbol{r})$
, projected on the
$xy$
plane (the computing procedure of PPDF in the short-time regime is discussed in appendix B). As shown in figure 4,
$g_{21}(\boldsymbol{r})$
in the
$xy$
plane exhibits streaks showing an intensified distribution particularly near the RBC surface under high shear rate and high haematocrit, while in contrast no significant string structure is found in low shear rate (
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 100~\text{s}^{-1}$
) (see supplementary movie 1) or/and low haematocrit (
$\unicode[STIX]{x1D719}=0.1$
) cases, as shown in figure 4(a–c). Note that a suspension string-like structure often occurs in sheared monodisperse colloidal suspensions at equilibrium as a consequence of the absence of interparticle lubrication interactions at high concentration (Xue & Grest Reference Xue and Grest1990; Foss & Brady Reference Foss and Brady2000b
) or the presence of long-range repulsive forces at low concentration (Kumar & Higdon Reference Kumar and Higdon2010). However, in the current case, ordering of NPs occurs at a non-equilibrium state that involves the change of the RBC suspension structure from a uniformly distributed and randomly oriented configuration to a streamwise-aligned and elongated configuration under high shear and high haematocrit. The sublinear growth of MSD is followed by a super-ballistic behaviour (
${\sim}t^{3}$
) at high shear rate (
$\dot{\unicode[STIX]{x1D6FE}}=10000~\text{s}^{-1}$
), which might be associated with the large jump of NPs between strings when the NP phase gradually evolves from the string-like structure towards the more uniform structure featuring in the long-time diffusive regime.
In the long-time regime (
$t\dot{\unicode[STIX]{x1D6FE}}>100$
), MSD reaches the second linear growth stage (
${\sim}t$
), where RESID becomes the dominate diffusive mechanism leading to a uniform distribution of NP (see lower inset of figure 3
c). Such a three-stage (diffusive/super-ballistic/diffusive) anomalous dispersion behaviour has previously been reported in the dispersion of Brownian particles subjected to external forces (Siegle, Goychuk & Hanggi Reference Siegle, Goychuk and Hanggi2010). Here, we show that such dispersive anomalies also occur in a sheared NP–RBC bidisperse suspension, where the anomalous NP dispersion is, however, driven by an internal shear-induced mechanism, i.e. the elongation and the alignment of the concentrated RBCs along the streamwise direction under high shear. Such transient shear-induced morphological adaptation of RBCs contributes to extra mobility of the NP phase, seemingly playing a role of external forces exerted on the NP phase.
Figures 3(e,f) depict the longitudinal and off-diagonal MSD evolution at
$\dot{\unicode[STIX]{x1D6FE}}=10\,000~\text{s}^{-1}$
. Results for lower shear rate show similar MSD behaviours and are not presented for discussion. Because the initial transient regime is neglected, the longitudinal and off-diagonal MSDs yield classical ballistic–diffusive transitions similar to the dispersive behaviour of rigid-particle suspensions under shear (Foss & Brady Reference Foss and Brady1999; Clausen et al.
Reference Clausen, Reasor and Aidun2011). The off-diagonal MSD under high shear starts with positive values during the ballistic regime, exhibits a cross-over transition involving a change of sign, and eventually maintains negative values in the diffusive regime. Absolute values are shown for the off-diagonal MSDs with the sign at certain temporal stage denoted in figure 3(f). The change of sign from positive (
$+$
) to negative (
$-$
) is a hallmark of the dominant NP migration direction shifted from along the extensional axes (first and third quadrants) to along the compressive axes (second and fourth quadrants) of the flow, as also observed in sheared colloidal suspensions (Foss & Brady Reference Foss and Brady1999).
Normalized particle long-time RESID in (a)
$xx$
, (b)
$yy$
, (c)
$zz$
and (d)
$xy$
directions against normalized shear rate (
$Pe$
or
$Ca_{G}$
) in a log–log scale. Each curve at specific
$\unicode[STIX]{x1D719}$
forms by connecting
$\hat{\unicode[STIX]{x1D60B}}_{ij}^{R}$
data at shear rate
$\dot{\unicode[STIX]{x1D6FE}}=10$
, 30, 100, 300, 1000, 2000 and
$10\,000~\text{s}^{-1}$
, from left to right;
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
at
$\dot{\unicode[STIX]{x1D6FE}}=10~\text{s}^{-1}$
is not shown for scaling purpose due to its small magnitude. Experimental results for RBC enhanced solute (Diller, Mikic & Drinker Reference Diller, Mikic and Drinker1980; Wang & Keller Reference Wang and Keller1985) and platelet (Grabowski, Friedman & Leonard Reference Grabowski, Friedman and Leonard1972; Antonini et al.
Reference Antonini, Guiffant, Quemada and Dosne1978) diffusivity in the velocity-gradient (
$yy$
) direction are plotted for comparison. NPs of size
$2a_{1}=100~\text{nm}$
yield
$\unicode[STIX]{x1D60B}^{B}=3.78~\unicode[STIX]{x03BC}\text{m}^{2}~\text{s}^{-1}$
.

3.3 Haemorheological dependence of NP long-time diffusivity
In this section, we focus on examining the NP long-time diffusive behaviour. The NP long-time diffusivities are evaluated in the long-time regime (after 100 strain units) and listed in table 1. Experimental statistics from various sources (Grabowski et al.
Reference Grabowski, Friedman and Leonard1972; Antonini et al.
Reference Antonini, Guiffant, Quemada and Dosne1978; Diller et al.
Reference Diller, Mikic and Drinker1980; Wang & Keller Reference Wang and Keller1985) are selected for comparison to gain credibility of the simulation results. Given the distinct time scales associated with the Brownian (
$\unicode[STIX]{x1D70F}_{B}\sim 10^{-4}~\text{s}$
) and the long-time RBC-enhanced diffusion (
$\unicode[STIX]{x1D70F}_{R}\sim 100/\dot{\unicode[STIX]{x1D6FE}}$
) phenomenon, the coupling of BD and RESID follows a simple superposition,
$\unicode[STIX]{x1D60B}_{ij}^{\infty }=\unicode[STIX]{x1D60B}_{ij}^{R}+\unicode[STIX]{x1D60B}^{B}\unicode[STIX]{x1D6FF}_{ij}$
, as confirmed in previous studies (e.g. Liu et al.
Reference Liu, Zhu, Rao, Clausen and Aidun2018), where
$\unicode[STIX]{x1D60B}_{ij}^{R}$
denotes the RESID. Since
$\unicode[STIX]{x1D60B}^{B}$
at infinite dilution follows the SE relation with negligible dependence on flow conditions, the haemorheological response of
$\unicode[STIX]{x1D60B}_{ij}^{\infty }$
is essentially determined by
$\unicode[STIX]{x1D60B}_{ij}^{R}$
.
3.3.1 Shear-rate dependence
Figure 5 plots the normalized RESID,
$\hat{\unicode[STIX]{x1D60B}}_{ij}^{R}=\unicode[STIX]{x1D60B}_{ij}^{\infty }/\unicode[STIX]{x1D60B}^{B}-\unicode[STIX]{x1D6FF}_{ij}$
, against normalized shear rate,
$Pe$
(or
$Ca_{G}$
). To first validate the simulation results, experimental results of solute or cell velocity-gradient (
$yy$
) diffusivity in sheared human or animal blood are selected for comparison; see figure 5(b). Wang & Keller (Reference Wang and Keller1985) measured the augmentation of ferricyanide solute diffusivity in both bovine and human RBC suspensions using a rotating Couette flow device. Diller et al. (Reference Diller, Mikic and Drinker1980) measured the enhanced radial diffusivity of oxygen solute in human blood using a tube oxygenator device. Grabowski et al. (Reference Grabowski, Friedman and Leonard1972) calculated the platelet velocity-gradient diffusivity in cavine (guinea pig) blood flowing through a channel by measuring the platelet rate of adhesion/deposition to a foreign surface attached to the flow chamber. A similar experiment was later conducted by Antonini et al. (Reference Antonini, Guiffant, Quemada and Dosne1978) to measure the platelet radial diffusivity in human blood. These measured diffusivities in 40 % haematocrit normalized based on our notation compares favourably with the calculated
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
at
$\unicode[STIX]{x1D719}=0.4$
based on the simulation, as shown in figure 5(b).
Depending on the level of shear rate imposed at various haematocrits, the diffusion tensor of NP in sheared blood shows different shear-rate dependence. At low shear rates (
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 100$
), all RESID terms exhibit a linear dependence on shear rate (
${\sim}\dot{\unicode[STIX]{x1D6FE}}$
), matching the linear
$\dot{\unicode[STIX]{x1D6FE}}$
scaling of shear-induced diffusivity in rigid-particle suspensions (Foss & Brady Reference Foss and Brady1999; Sierou & Brady Reference Sierou and Brady2004). This is also consistent with the insignificant RBC morphological changes at
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 100~\text{s}^{-1}$
, as shown in figure 2(a–c).
At intermediate shear rates (
$100\leqslant \dot{\unicode[STIX]{x1D6FE}}\leqslant 2000$
), significant streamwise elongation of RBCs occurs with increasing shear rate; see figure 2(c–e). As a result, nonlinear
$\dot{\unicode[STIX]{x1D6FE}}$
scaling is observed in all diagonal RESID terms. Specifically, cross-stream diffusivities,
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
and
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
, show sublinear scales (
${\sim}\dot{\unicode[STIX]{x1D6FE}}^{m}$
,
$m=0.7{-}0.8$
), while streamwise diffusivity,
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
, exhibits superlinear scales (
${\sim}\dot{\unicode[STIX]{x1D6FE}}^{n}$
,
$n=1{-}1.8$
). In contrast to the nonlinear scaling in diagonal diffusivities,
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
maintains largely a linear scale at the intermediate shear-rate regime, as shown in figure 5(d). It is also noted that such nonlinear shear-rate scaling in diagonal RESID terms is most prominent at intermediate to high haematocrit (
$\unicode[STIX]{x1D719}>0.1$
), which implies that RBC deformability plays a less important role in altering the RESID shear dependence at low haematocrits. This observation is consistent with the results of shear-augmented solute diffusivity in model-RBC suspensions at various RBC deformability reported by Wang & Keller (Reference Wang and Keller1985), where they found changing RBC deformability barely affects the augmentation at a particle volume concentration of 0.1 or less. Wang & Keller (Reference Wang and Keller1985) also find a sublinear scaling of
${\sim}\dot{\unicode[STIX]{x1D6FE}}^{\unicode[STIX]{x1D6FD}}$
(
$0.67\leqslant \unicode[STIX]{x1D6FD}\leqslant 0.89$
) for the solute diffusivity in velocity-gradient directions, while our simulations observe an exponent of
$\unicode[STIX]{x1D6FD}=0.7$
.
At high shear rates (
$\dot{\unicode[STIX]{x1D6FE}}\geqslant 2000$
) and intermediate to high haematocrit (
$\unicode[STIX]{x1D719}>0.1$
), all RESID terms except
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
exhibit reduced shear-rate dependence compared to the intermediate shear-rate regime. In the velocity-gradient (
$y$
) direction, the hindrance of diffusion is due to the fact that the concentrated RBCs become more aligned and elongated with flow such that RBCs act as obstacles against the NP cross-stream diffusion in the
$y$
direction; see figure 2(f). Such an effect, however, shows less hindrance on
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
, as it does not forbid the NP migration in the vorticity (
$z$
) direction. In the streamwise (
$x$
) direction, the reduction of the shear dependence of
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
is likely to be associated with the saturation of the RBC elongation. Owing to the compound effects in both
$x$
and
$y$
directions, the off-diagonal diffusivity,
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
, also exhibits a certain reduction of the shear-rate dependence.
3.3.2 Haematocrit dependence
Normalized particle long-time RESID in (a)
$xx$
, (b)
$yy$
, (c)
$zz$
and (d)
$xy$
directions against haematocrit under different shear rate in logarithm scale. Lines are the best power-law fits to the data points at specific
$\dot{\unicode[STIX]{x1D6FE}}$
. The
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
value at
$\dot{\unicode[STIX]{x1D6FE}}=10~\text{s}^{-1}$
(not shown) is positive with negligible magnitude compared to other diffusion terms. NPs of size
$2a_{1}=100~\text{nm}$
yield
$\unicode[STIX]{x1D60B}^{B}=3.78~\unicode[STIX]{x03BC}\text{m}^{2}~\text{s}^{-1}$
.

Figure 6 displays the same data as figure 5 but plotted against
$\unicode[STIX]{x1D719}$
to show the haematocrit dependence of RESID at different shear rate. At low shear rates (
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 100$
), all diffusion terms manifest a linear
$\unicode[STIX]{x1D719}$
scaling. Note that a linear concentration dependence of particle self-diffusivity has been observed in sheared monodisperse suspensions in the presence of surface roughness (Da Cunha & Hinch Reference Da Cunha and Hinch1996) or residual Brownian motion (Brady & Morris Reference Brady and Morris1997), which causes the two-body interaction to be irreversible and hence give rise to a diffusive behaviour. Since the shear-induced diffusion of NP is driven by NP–RBC interaction, it is likely that the
$\unicode[STIX]{x1D719}$
scaling of the NP diffusion tensor results from the irreversible two-body interaction between NP and RBC induced by the RBC membrane roughness/flexibility and the NP Brownian effect. The off-diagonal diffusivity
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
is found to be positive (
$+$
) at very low shear rate (
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 10~\text{s}^{-1}$
) and otherwise negative (
$-$
). The change-of-sign behaviour of the off-diagonal diffusivity designates the dominant displacement direction of NP changing from along the extensional axes to along the compressive axes of the flow as shear rate increases. A similar observation has been reported in sheared monodisperse colloidal suspensions (Foss & Brady Reference Foss and Brady1999).
As shear rate grows above
$100~\text{s}^{-1}$
, various
$\unicode[STIX]{x1D719}$
scaling arises in different RESID components. The
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
component exhibits a quadratic
$\unicode[STIX]{x1D719}$
scaling (
${\sim}\unicode[STIX]{x1D719}^{2}$
), which suggests that in the
$x$
direction additional effects exist to drive the NP diffusive motion besides the irreversible pairwise interactions. The
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
and
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
components scale linearly with
$\unicode[STIX]{x1D719}$
, suggesting that the irreversible pairwise interaction remains the dominant diffusive mechanism. The
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
component exhibits a transition from linear to sublinear
$\unicode[STIX]{x1D719}$
scale (
${\sim}\unicode[STIX]{x1D719}^{0.6}$
). The mostly linear
$\unicode[STIX]{x1D719}$
scale in
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
observed in the current numerical study is consistent with the experimental observation in Wang & Keller (Reference Wang and Keller1985), where they show that the solute diffusivity in the velocity-gradient (
$y$
) direction is augmented by approximately threefold when the particle concentration increases from 0.1 to 0.31. Since suspension viscosity typically increases with particle concentration (Foss & Brady Reference Foss and Brady2000b
), the enhancement of solute (e.g. NP) diffusivity at increased
$\unicode[STIX]{x1D719}$
indicates that the RESID is due to the RBC–NP interaction rather than the secondary flow effect (Wang & Keller Reference Wang and Keller1985).
RBC–NP partial pair distribution function,
$g_{21}(\boldsymbol{r})$
, projected onto the
$xy$
,
$xz$
and
$yz$
planes in the long-time regime under various haemorheological conditions, where the PPDF contour has the horizontal edge aligned with the first axis (e.g.
$xy$
PPDF contour has the horizontal edge aligned with the
$x$
axis). The edge length of the sampling box for computing
$g_{21}(\boldsymbol{r})$
is three times the maximum diameter of the undeformed RBC.

3.4 Microstructure
The rheological properties of particle suspensions are often determined by the mechanistic phenomenon occurring on the particle length scales. To elucidate the physical mechanisms that govern the haemorheological scaling behaviours of the NP diffusion tensor observed in § 3.3, the RBC–NP PPDF,
$g_{21}(\boldsymbol{r})$
, in the long-time regime is plotted to visualize the configurational microstructure of the RBC–NP bidisperse suspension under various haematocrits and shear rates. The techniques used to compute and project
$g_{21}(\boldsymbol{r})$
are discussed in detail in appendix B.
Figure 7 presents a matrix of
$g_{21}(\boldsymbol{r})$
projections onto the velocity–velocity gradient (
$xy$
) plane, the velocity–vorticity (
$xz$
) plane and the velocity gradient–vorticity (
$yz$
) plane under two shear rates (
$\dot{\unicode[STIX]{x1D6FE}}=100$
or
$2000~\text{s}^{-1}$
) and two haematocrits (
$\unicode[STIX]{x1D719}=0.1$
or 0.4). In general, all PPDFs feature a large-scale rhombus shape as opposed to the typical circular shape commonly observed in rigid-sphere-particle suspensions (Foss & Brady Reference Foss and Brady2000a
; Wang & Brady Reference Wang and Brady2016; Pednekar, Chun & Morris Reference Pednekar, Chun and Morris2018). This can be attributed to the disk shape of RBC that causes geometry-specific anisotropy of the microstructure. The average RBC shapes under various haemorheological conditions, as shown in the central low-intensity region of the PPDF contours, are nicely captured through the PPDF sampling procedure. All PPDFs considered in the
$xy$
plane show more intensified distribution (i.e. higher probability of RBC–NP interaction) near the RBC disk surfaces, of which the surface normal directions are more aligned with the compressive axes (in second and fourth quadrants) in accordance with the negative values of
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
. In the
$xz$
and
$yz$
plane, the PPDF distribution tends to be symmetric about the principal axes (
$x$
,
$y$
or
$z$
), which explains the zero values of the
$xz$
and
$yz$
diffusivities.
As the haemorheological condition changes, the detailed configuration of the PPDF within the rhombus structure also varies. At low haematocrit (
$\unicode[STIX]{x1D719}=0.1$
), the
$xy$
PPDF shows fore–aft low intensity similar to rigid-particle suspensions (Kumar & Higdon Reference Kumar and Higdon2011). However, the break of the fore–aft symmetry, being different from the rigid-particle suspensions, seems to be related to the RBC orientation aligned with the extensional flow axes that is further caused by the tank-treading motion of the RBC membrane (Reasor et al.
Reference Reasor, Clausen and Aidun2013a
). As the haematocrit increases to
$\unicode[STIX]{x1D719}=0.4$
, the overall PPDF around the RBC becomes more intensified, meaning the local NP concentration near the RBC surface increases. More interestingly, the fore–aft low PPDF region observed under low haematocrit gets intensified substantially. The change of the NP microstructure with increased haematocrit can be explained by smaller inter-cell separation and hence NP getting squeezed in a smaller inter-cell region. At low shear rate (
$\dot{\unicode[STIX]{x1D6FE}}=100~\text{s}^{-1}$
), RBC shows a close-to-equilibrium biconcave shape. Increasing the shear rate to
$\dot{\unicode[STIX]{x1D6FE}}=2000~\text{s}^{-1}$
results in a significant elongation of the average RBC shape in the flow (
$x$
) direction accompanied by certain contraction in the velocity-gradient (
$y$
) and vorticity (
$z$
) directions.
The above configurational changes of
$g_{21}(\boldsymbol{r})$
under various shear rates and haematocrits provide possible mechanistic explanations for the nonlinear
$\unicode[STIX]{x1D719}$
dependence observed at high shear rates, as discussed in § 3.3.2. The quadratic dependence of
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
(
${\sim}\unicode[STIX]{x1D719}^{2}$
), shown in figure 6(a), is probably due to the occurrence of three-body RBC–NP interactions in the longitudinal (
$x$
) direction caused by the compound effect of the elevated NP concentration at the RBC fore–aft surface and the elongation of the RBC in the
$x$
direction. The reduced
$\unicode[STIX]{x1D719}$
dependence in
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
at high shear rate, as shown in figure 6(c), can be attributed to the relatively large contraction of the RBC in the vorticity direction, which reduces the effective
$\unicode[STIX]{x1D719}$
in the
$z$
direction. The lesser reduction of the
$\unicode[STIX]{x1D719}$
dependence in the
$y$
direction, as shown in figure 6(c), is due to the
$y$
direction contraction of RBC under shear being less significant than that in the
$z$
direction, as clearly indicated in the PPDF contours.
3.5 Role of RBC deformability
The above PPDF analysis shows prominent RBC morphological change with elevated shear rate, which suggests that RBC deformability may play an important role in causing the nonlinear
$\dot{\unicode[STIX]{x1D6FE}}$
dependence of the NP diffusion tensor. In this section, we perform numerical experiments to further quantitatively explain the nonlinear
$\dot{\unicode[STIX]{x1D6FE}}$
scaling of RESID in the intermediate-shear-rate regime (
$100\leqslant \dot{\unicode[STIX]{x1D6FE}}\leqslant 2000~\text{s}^{-1}$
), as observed in § 3.3.1. Changing shear rate alters both the fluid inertia and RBC deformability, quantified by
$Pe$
and
$Ca_{G}$
, respectively. To interrogate the isolated effect of RBC deformability (
$Ca_{G}$
), we fix
$Ca_{G}$
by scaling up
$G$
while increasing
$Pe$
(through increasing
$\dot{\unicode[STIX]{x1D6FE}}$
). Two capillary numbers,
$Ca_{G}=0.055$
and
$0.55$
, are considered corresponding to the
$Ca_{G}$
regime where the nonlinear shear-rate dependence of RESID occurs. For these simulations, we select a fixed haematocrit of
$\unicode[STIX]{x1D719}=0.4$
and an NP size of
$2a_{1}=100~\text{nm}$
. Table 2 lists all parameters and the NP diffusivity values for the cases tested.
The
$Pe$
dependence of the (a) diagonal and (b) off-diagonal RESID terms at fixed capillary number (by rescaling
$G$
) with
$\unicode[STIX]{x1D719}=0.4$
. Two RBC capillary numbers,
$Ca_{G}=0.55$
and
$0.055$
, are considered to see the isolated effect of changing RBC deformability (
$Ca_{G}$
) or adjusting the inertial effect (
$Pe$
).

Simulation data for numerical experiments concerning the isolated RBC deformation effect in NP–RBC suspensions. The RBC membrane shear modulus,
$G$
, is rescaled to numerically match specific
$Ca_{G}$
while changing shear rate. Here, the NP size is fixed to
$2a_{1}=100~\text{nm}$
; RBC has an effective radius of
$a_{2}=2.9~\unicode[STIX]{x03BC}\text{m}$
. Brownian diffusivity is calculated by
$\unicode[STIX]{x1D60B}^{B}=k_{B}T/6\unicode[STIX]{x1D707}\unicode[STIX]{x03C0}a_{1}$
at a temperature of
$T=310$
K.

As presented in figure 8(a,b), instead of exhibiting nonlinear
$Pe$
scaling as observed in § 3.3.1, all RESID components show a linear dependence on
$Pe$
at fixed RBC capillary number of
$Ca_{G}=0.55$
or
$0.055$
across the wide range of shear rates. This recovers the linear
$Pe$
dependence of solute/self-diffusivity in sheared rigid-particle suspensions (Zydney & Colton Reference Zydney and Colton1988; Foss & Brady Reference Foss and Brady1999; Sierou & Brady Reference Sierou and Brady2004) since we intentionally fix the RBC deformation (
$Ca_{G}$
) while adjusting
$Pe$
.
Moreover, an increase of
$Ca_{G}$
(from 0.055 to 0.55) leads to an increase of
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
but a decrease of
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
and
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
, as shown in figure 8(a). These results are direct evidence indicating that the nonlinear shear-rate scaling of the diagonal RESID terms observed in § 3.3.1 is due to the variation of RBC deformability (
$Ca_{G}$
) induced by changing shear rate at fixed
$G$
. More specifically, the increase of
$Ca_{G}$
with increasing shear rate at
$100<\dot{\unicode[STIX]{x1D6FE}}<2000~\text{s}^{-1}$
under fixed
$G$
leads to the superlinear shear-rate dependence of longitudinal diffusivity and meanwhile the sublinear dependence of cross-stream diffusivities. In figure 8(b), the off-diagonal diffusivity shows insignificant change at different
$Ca_{G}$
, consistent with its mostly linear shear-rate dependence shown in figure 5(d).
3.6 Relevant length scale
Although the shear adaptation of the RBC deformability is shown to be responsible for the nonlinear shear-rate dependence of the RESID, it is still unclear in what way it alters the RESID. In this section, we further identify the characteristic length scale relevant to the RESID to gain a more in-depth understanding of this nonlinear phenomenon.
In a series of experiments measuring the particle diffusivity in a concentrated non-colloidal suspension, Breedveld et al. (Reference Breedveld, Van den Ende, Tripathi and Acrivos1998) show that the ratio of the diffusivity of fluid tracers to the self-diffusivity of non-colloidal particles is close to unity. Through DNS of the cellular blood flow in a microchannel, Zhao et al. (Reference Zhao, Shaqfeh and Narsimhan2012) show that the cross-stream diffusivity of platelets (
${\sim}2~\unicode[STIX]{x03BC}\text{m}$
) in blood exhibits similar magnitude to that of passive tracers. These findings imply that the shear-induced diffusivity is insensitive to the size of the scarce (in terms of volume fraction) small particles. Following these results, we hypothesize that RESID is insensitive to the size of NPs. Moreover, the characteristic length scale associated with the RESID should be the size of RBCs. To confirm this hypothesis, we evaluate the RESID at fixed haemorheological condition with various NP size. The NP size is kept in submicrometre scale such that the NP–RBC size ratio satisfies
$a_{1}/a_{2}\ll 1$
. The volume fraction of the NP phase is kept below 0.32 % to satisfy the dilute condition. The set-up and measured diffusivities are listed in table 3.
In figure 9(a), we plot the RESID versus
$Pe$
at
$\dot{\unicode[STIX]{x1D6FE}}=1000~\text{s}^{-1}$
and
$\unicode[STIX]{x1D719}=0.2$
. Here,
$Pe$
is adjusted by NP size in the range of
$50\leqslant 2a_{1}\leqslant 800~\text{nm}$
. A sublinear
$\hat{\unicode[STIX]{x1D60B}}_{ij}^{R}$
–
$Pe$
relationship is observed for all diffusion coefficients. Figure 9(b) further plots the same data as figure 9(a), but with
$Pe$
rescaled by the NP–RBC size ratio as
This simple rescaling leads to a strong linear relationship between all RESID terms and the rescaled Péclet number,
$\widetilde{Pe}$
. Moreover, the observed linear relationship,
$(\unicode[STIX]{x1D60B}_{ij}^{\infty }/\unicode[STIX]{x1D60B}^{B}-\unicode[STIX]{x1D6FF}_{ij})\sim O(\widetilde{Pe})$
, can be deduced to show that
$\unicode[STIX]{x1D60B}_{ij}^{R}\sim O(\dot{\unicode[STIX]{x1D6FE}}a_{2}^{2})$
, which indicates that the dimensional RESID is indeed insensitive to the NP size within the range of NP sizes considered in the current study. This also suggests that NPs at submicrometre size play a secondary role in affecting the NP apparent diffusivity primarily through altering BD.
RESID plotted against (a) Péclet number
$Pe$
and (b) the rescaled Péclet number
$\widetilde{Pe}$
, at
$\dot{\unicode[STIX]{x1D6FE}}=1000~\text{s}^{-1}$
and
$\unicode[STIX]{x1D719}=0.2$
in logarithmic scales. Both
$Pe$
and
$\widetilde{Pe}$
are adjusted with NP size in the range of
$2a_{1}=50{-}800~\text{nm}$
. (c) Linear dependence of RESID on
$\widetilde{Pe}$
at fixed
$Ca_{G}=0.55$
and
$\unicode[STIX]{x1D719}=0.2$
;
$\widetilde{Pe}$
is adjusted by shear rate (with
$2a_{1}=100~\text{nm}$
) or NP size (with
$\dot{\unicode[STIX]{x1D6FE}}=1000~\text{s}^{-1}$
). The diffusion tensor in monodisperse colloidal suspensions (Foss & Brady Reference Foss and Brady1999) is also plotted for comparison, where
$\widetilde{Pe}$
degrades to
$Pe$
.

Simulation data for numerical experiments concerning the relevant length scale in NP–RBC suspensions. NP size is varied in the range of
$2a_{1}=50{-}800~\text{nm}$
. RBC has an effective radius of
$a_{2}=2.9~\unicode[STIX]{x03BC}\text{m}$
. Brownian diffusivity is calculated by
$\unicode[STIX]{x1D60B}^{B}=k_{B}T/6\unicode[STIX]{x1D707}\unicode[STIX]{x03C0}a_{1}$
at temperature
$T=310$
K. For the
$800~\text{nm}$
case, the number of NPs is set to 500 to satisfy the dilute condition for the NP volume fraction. The last three cases are performed with adjusted
$G$
to fix
$Ca_{G}$
.

To further justify the universality of the rescaled Péclet number,
$\widetilde{Pe}$
, i.e. whether the RBC size is the reasonable length scale associated with RBC-enhanced diffusion or not, we plot RESID against
$\widetilde{Pe}$
at
$Ca_{G}=0.55$
and
$\unicode[STIX]{x1D719}=0.2$
, with
$\widetilde{Pe}$
adjusted by either shear rate or NP size. As shown in figure 9(c), the data points for a particular RESID coefficient are found to be well aligned on the same linear
$\widetilde{Pe}$
scaling line. This directly confirms that the RBC size is the characteristic length scale governing the RBC-enhanced diffusion, and
$\widetilde{Pe}$
is a more general non-dimensional term quantifying the RESID.
Identifying RBC size as the RESID length scale helps in further understanding the nonlinear shear-rate dependence of the RESID. In figure 10, we present snapshots of RBC and NP distribution and morphology in the mid-cross-stream cross-sectional plane. As shown in figure 10, the
$x$
direction elongation of RBCs under high shear causes a contraction of RBCs in the
$yz$
plane, which causes a reduction of the average effective RBC size (i.e. the length scale
$\bar{a}_{2,yz}$
) and hence a decrease of the effective
$\widetilde{Pe}$
in cross-stream directions. This explains the sublinear shear-rate scaling of the cross-stream RESID when increasing
$Ca_{G}$
. Likewise, the superlinear shear-rate scaling of the longitudinal RESID can be attributed to the elongation of RBC in the streamwise direction, i.e. the increase of the RESID length scale in the longitudinal direction. The change of the RESID length scales in different principal directions due to shear is the root cause of the nonlinear shear-rate dependence of the NP diffusion tensor.
RBC and NP distributions and morphology in the cross-stream plane (
$y$
–
$z$
) with (a)
$Ca_{G}=0.055$
, (b)
$Ca_{G}=0.55$
and (c)
$Ca_{G}=1.1$
at
$Pe=0.066$
and
$\unicode[STIX]{x1D719}=0.4$
. Large non-circular disks are RBC cross-sectional contours. Small circles are NP cross-sectional contours. The
$\bar{a}_{2,yz}$
value can be estimated from the PPDF analysis.

3.7 Empirical correlations
In this section, we construct empirical correlations for the long-time NP diffusion tensor in sheared blood based on the scaling observations in the previous sections. Since all RESID coefficients scale linearly with
$\widetilde{Pe}$
at fixed
$Ca_{G}$
, and the nonlinear shear-rate dependence of RESID is primarily due to the variation of
$Ca_{G}$
, as demonstrated in §§ 3.5 and 3.6, there hence exists a scaling relation,
$\hat{\unicode[STIX]{x1D60B}}_{ij}^{R}\sim O(\unicode[STIX]{x1D719}^{p_{1}}\widetilde{Pe}Ca_{G}^{p_{2}})$
, such that the exponents
$p_{1}$
and
$p_{2}$
can be estimated through matching the nonlinear
$\unicode[STIX]{x1D719}$
and
$\dot{\unicode[STIX]{x1D6FE}}$
scalings, respectively, as observed in § 3.3. Therefore, empirical correlations of NP long-time diffusivities can be constructed as functions of
$\unicode[STIX]{x1D719}$
,
$\widetilde{Pe}$
and
$Ca_{G}$
in the haemorheological range of
$10\leqslant \dot{\unicode[STIX]{x1D6FE}}\leqslant 2000~\text{s}^{-1}$
and
$0.1\leqslant \unicode[STIX]{x1D719}\leqslant 0.4$
as follows:
where the anisotropic tensor,
$\unicode[STIX]{x1D648}^{q}$
, according to the severity of the RBC deformability (
$Ca_{G}=\unicode[STIX]{x1D707}\dot{\unicode[STIX]{x1D6FE}}a_{2}/G$
) at different level of shear rates, yields piecewise expressions:
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}^{l}=\left(\begin{array}{@{}ccc@{}}{\mathcal{C}}_{xx}^{l} & {\mathcal{C}}_{xy}^{l} & 0\\ {\mathcal{C}}_{yx}^{l} & {\mathcal{C}}_{yy}^{l} & 0\\ 0 & 0 & {\mathcal{C}}_{zz}^{l}\end{array}\right),\quad 10\leqslant \dot{\unicode[STIX]{x1D6FE}}\leqslant 100~\text{s}^{-1}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D648}^{h}=\left(\begin{array}{@{}ccc@{}}{\mathcal{C}}_{xx}^{h}\unicode[STIX]{x1D719}^{Ca_{G}-0.055}Ca_{G}^{0.8} & {\mathcal{C}}_{xy}^{h} & 0\\ {\mathcal{C}}_{yx}^{h} & {\mathcal{C}}_{yy}^{h}/Ca_{G}^{0.3} & 0\\ 0 & 0 & {\mathcal{C}}_{zz}^{h}/(\unicode[STIX]{x1D719}^{2Ca_{G}+0.11}Ca_{G})^{0.2}\end{array}\right), & \displaystyle \nonumber\\ \displaystyle & \displaystyle 100<\dot{\unicode[STIX]{x1D6FE}}\leqslant 2000~\text{s}^{-1}, & \displaystyle\end{eqnarray}$$
Measured NP RESID,
$\unicode[STIX]{x1D60B}_{ij}^{(R,m)}=\unicode[STIX]{x1D60B}_{ij}^{\infty }-\unicode[STIX]{x1D60B}^{B}\unicode[STIX]{x1D6FF}_{ij}$
, as a function of the theoretical RESID,
$\unicode[STIX]{x1D60B}_{ij}^{(R,th)}=\unicode[STIX]{x1D719}\dot{\unicode[STIX]{x1D6FE}}a_{2}^{2}\unicode[STIX]{x1D614}_{ij}^{q}$
, based on the empirical correlations presented in equations (3.6a
) and (3.6b
). For cases with
$\dot{\unicode[STIX]{x1D6FE}}\leqslant 100~\text{s}^{-1}$
(filled symbols),
$\unicode[STIX]{x1D614}_{ij}^{l}$
is applied to calculate
$\unicode[STIX]{x1D60B}_{ij}^{(R,th)}$
. For cases with
$\dot{\unicode[STIX]{x1D6FE}}>100~\text{s}^{-1}$
(open symbols),
$\unicode[STIX]{x1D614}_{ij}^{h}$
is applied to calculate
$\unicode[STIX]{x1D60B}_{ij}^{(R,th)}$
.

Correlation constants for the empirical correlations of NP diffusion tensor in sheared cellular blood flow.
Figure 11 plots the measured NP diffusivities (subtracted by BD) versus the theoretical estimation based on the empirical correlations. The good collapse of the diffusivity measurements on the
$\unicode[STIX]{x1D60B}_{ij}^{(R,m)}=\unicode[STIX]{x1D60B}_{ij}^{(R,th)}$
line demonstrates the empirical correlations can well reproduce the numerically measured diffusivity tensor of NP in sheared blood.
4 Conclusions
The dispersion of NP in cellular blood flow under unbounded, homogeneous shear has been investigated over a wide range of shear rate and haematocrit using an LB–LD–SL multiscale complex blood flow solver. In the short-time regimes, NP dispersive anomalies are observed and attributed to the transient morphology and orientation change of RBCs under high shear and high haematocrit. In the long-time regimes, results for the long-time diffusivity in the velocity-gradient direction agree well with existing experimental data. The long-time NP diffusion tensor has been described as a function of shear rate and haematocrit with various power-law scalings.
By plotting the RBC–NP PPDF, the NP microstructure in sheared blood has been visualized for the first time, and features a rhombus configuration with a detailed inner structure that changes according to specific haemorheological conditions. The RBC–NP PPDF analysis also suggests a novel approach to visualizing the average RBC morphology in concentrated RBC suspensions subject to different haemorheological conditions. Under high shear rate, the
$\unicode[STIX]{x1D719}^{2}$
dependence in
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
is proposed to be related to the
$x$
elongation of RBC and elevated RBC–NP PPDF near the fore–aft region, which together suggests possible more-than-two-body interaction occurring particularly in the streamwise direction. The sublinear
$\unicode[STIX]{x1D719}$
dependence is suggested to be related to the reduced effective
$\unicode[STIX]{x1D719}$
due to the substantial contraction of RBC in the vorticity direction subject to high shear. The RBC–NP bidisperse suspension presents an example of a highly anisotropic microstructure of particle suspensions caused by the compound effect of particle shape/orientation anisotropy and shear-flow anisotropy.
It is also found that there exists a critical shear rate (
${\sim}100~\text{s}^{-1}$
) around which the RESID shear-rate dependence changes from linear to nonlinear scale. Through numerical experiments, the transition to nonlinear shear-rate scaling of RESID has been related to the prominent change of average RBC morphological state between different shear rates. Specifically, the superlinear scalings (
$\dot{\unicode[STIX]{x1D6FE}}^{1-1.8}$
) of
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
are due to the streamwise elongation of RBCs, while the sublinear scalings (
$\dot{\unicode[STIX]{x1D6FE}}^{0.7-0.8}$
) of
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
and
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
are associated with the cross-stream contraction of RBCs (in response to the streamwise elongation). The morphological changes under shear alter the RESID length scale in different principal directions, which has been demonstrated to be the fundamental cause of the nonlinear shear-rate dependence of the RESID. This mechanism is also worth distinguishing from the causes of the nonlinear shear-rate scaling of the self-diffusivity of RBCs (Gross et al.
Reference Gross, Kruger and Varnik2014; Mountrakis et al.
Reference Mountrakis, Lorenz and Hoekstra2016) or deformable capsules (Clausen et al.
Reference Clausen, Reasor and Aidun2011), where the latter has been attributed to the heterogeneous interparticle ‘collision’ due to the cell deformability (Kumar & Graham Reference Kumar and Graham2012). However, the nonlinear shear-rate dependence of NP diffusion in sheared blood, based on our interrogation, is more associated with a ‘one-way’ mechanism, i.e. the RBC morphological adaptation to shear flow changes the RESID length scale, which further alters the NP diffusion rate.
The determination of the rescaled Péclet number being a more general non-dimensional term to describe the severity of RESID enables the comparison between the bidisperse NP–RBC suspension and the monodisperse colloidal suspensions. In the latter scenario where the particle size ratio is one,
$\widetilde{Pe}$
drops to
$Pe$
and ‘RESID’ drops to the particle self-diffusivity. The self-diffusion tensor reported by Foss & Brady (Reference Foss and Brady1999) in a sheared monodisperse colloidal suspension is plotted in figure 9(c) for comparison. In both monodisperse and bidisperse scenarios,
$\hat{\unicode[STIX]{x1D60B}}_{xx}^{R}$
shows the greatest magnitude among all diffusivity terms; and
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
is slightly greater than
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
. In general,
$\hat{\unicode[STIX]{x1D60B}}_{ij}^{R}$
in RBC–NP suspension shows higher anisotropy than the monodisperse case owing to the geometric asymmetry of RBCs (
$x$
being greater than
$y$
and
$z$
dimensions). Such a geometric asymmetry effect can be further increased with
$Ca_{G}$
, leading to higher anisotropy of the diffusivity tensor. In the monodisperse scenario,
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
is the smallest; while in the NP–RBC bidisperse scenario,
$\hat{\unicode[STIX]{x1D60B}}_{xy}^{R}$
is greater than
$\hat{\unicode[STIX]{x1D60B}}_{yy}^{R}$
and
$\hat{\unicode[STIX]{x1D60B}}_{zz}^{R}$
due to the severe diffusive effect in the
$x$
direction.
This work, to the authors’ knowledge, offers the first detailed study of the complete 3D NP diffusion tensor in cellular blood flow over a wide range of shear rates and haematocrits. The proposed empirical correlations for the NP diffusion tensor offers a constitutive relation that can be adopted by effective continuum models to pursue large-scale NP biotransport applications (e.g. in vivo NP drug delivery) with better accuracy.
Acknowledgements
The authors acknowledge the support from Sandia National Laboratories under grant 2506X36 and the computational resource provided by National Science Foundation under grant TG-CT100012. The authors appreciate the suggestion by one of the reviewers to consider the pair distribution function of NPs around RBCs in all three planes. Z.L. acknowledges the constructive discussions with Dr J. B. Lechman, Professor E. C. Eckstein and Professor K. B. Wiesenfeld. Sandia National Laboratories is a multimission laboratory managed and operated by National Technology and Engineering Solutions of Sandia LLC, a wholly owned subsidiary of Honeywell International Inc. for the US Department of Energy’s National Nuclear Security Administration under contract DE-NA0003525. This paper describes objective technical results and analysis. Any subjective views or opinions that might be expressed in the paper do not necessarily represent the views of the US Department of Energy or the United States Government.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2019.320.
Appendix A. Sensitivity to NP–RBC contact model
The NP–RBC short-distance interaction is through a Morse potential that forbids NP from penetrating the RBC membrane. The potential parameters are adjusted to match the measured inter-cell potential energy, as discussed in Liu et al. (Reference Liu, Zhang, Wang and Liu2004). The Morse potential function is given as
where
$r$
is the normal distance between the particle centre to the RBC surface,
$r_{0}$
is a cutoff distance in which no interaction forces are present,
$D_{e}$
is the potential well depth and
$\unicode[STIX]{x1D6FD}$
is a scaling factor. The surface energy is set to
$D_{e}=1\times 10^{7}k_{B}T$
and the equilibrium distance is set to
$r_{0}=a_{1}+10~\text{nm}$
. In figure 12, we explore the sensitivity of the NP diffusivity calculation to the adjustment of the Morse model parameters. Figures 12(a) and 12(b) show that the
$\unicode[STIX]{x1D60B}_{yy}^{\infty }$
and
$\unicode[STIX]{x1D60B}_{zz}^{\infty }$
calculation exhibits less than 5 % variation when changing the magnitude of the energy and equilibrium distance, respectively, by up to
$\pm 60\,\%$
. This indicates that the NP diffusion is largely driven by the hydrodynamic interaction rather than by the direct contact between NP and RBC membrane.
(a) NP long-time diffusivities in the cross-stream directions plotted against the relative variation of the Morse potential well depth (a) and equilibrium distance (b).

Appendix B. Calculation of the partial pair distribution function
In the bidisperse RBC–NP suspension system, the NP–RBC partial pair distribution function (PPDF),
$g_{21}(\boldsymbol{r})$
, quantifies the conditional probability of finding an NP (species 1) at a position of
$\boldsymbol{r}$
with respect to the geometric centre of a single RBC (species 2). This quantity can be calculated as
where
$L_{b}$
defines the cubic box size for PPDF sampling around one RBC, index
$j$
goes through all NPs within the local sampling box, and
$i$
goes through all RBCs in the computational LB domain. Notation
$r_{2}^{i}$
and
$r_{1}^{j}$
denotes the positions of RBC and NP, respectively;
$n_{1}^{i}$
shows the number concentration of NP within the sampling box; while
$n_{2}$
denotes the number concentration of RBC in the entire domain. The angle brackets represent the ensemble average among independent realizations, which in the current case is through time averaging given the ergodic hypothesis. Similar techniques have recently been used in calculating PPDFs in bidisperse and polydisperse rigid-particle suspensions (Wang & Brady Reference Wang and Brady2016; Pednekar et al.
Reference Pednekar, Chun and Morris2018). Projection of
$g_{21}(\boldsymbol{r})$
to principal planes follows the integration procedure discussed in Kumar & Higdon (Reference Kumar and Higdon2011) but with a smaller integration interval, (
$-a_{2}/4$
,
$a_{2}/4$
), to better capture the RBC morphological change.
The PPDF sampling box size for all cases is selected to be three times the RBC maximum diameter at equilibrium in biconcave shape, i.e.
$L_{b}=24~\unicode[STIX]{x03BC}\text{m}$
. To capture the PPDF in the long-time regime (
$t\dot{\unicode[STIX]{x1D6FE}}>100$
) with detailed suspension microstructure, a total number of 300 strain units are employed for time averaging. In the short-time regime (
$t\dot{\unicode[STIX]{x1D6FE}}\sim 1$
), owing to limited time steps associated with specific suspension configuration (e.g. the string ordered configuration), approximately three strain units are adopted for averaging.















