Numerical Study of Particle Margination in a Square Channel Flow with Red Blood Cells

Abstract

Red blood cells flow near the axis in a small vessel, known as axial accumulation. This causes a region called the cell-free layer, which does not contain red blood cells near the wall. Then, small particles such as platelets come out to the cell-free layer. This phenomenon is called platelet margination. In this study, related to this phenomenon, direct numerical simulations were conducted using the immersed boundary method. The effects of the shear rate, channel size, and hematocrit value were investigated on the pressure-driven flow in a straight tube with a square cross-section. The simulation results indicated that the margination rate, which is the ratio of the distance traveled in the flow direction to the margination distance in the wall direction, is independent of the shear rate. The effect of the channel size on platelet margination was found to be well scaled by introducing a dimensionless parameter, which included the shear rate and effective area of the particle movement. It was also found that the margination rate varied nonlinearly with the tube hematocrit. This was due to the volume exclusion effect of red blood cells, which facilitated or hindered the motion of particles depending on the hematocrit. The relationship between the stable position of the particles near the corner and the width of the cell-free layer was also found. Furthermore, velocity fluctuations normalized by wall shear rate in a cross-section collapsed to one curve in the presented simulations. This indicates that the lateral force acting on the particles increases linearly with the shear rate.

1. Introduction

Blood flow consists of approximately 55% plasma, 45% Red Blood Cells (RBCs), and less than 1% platelets or white blood cells. Plasma can be treated as a Newtonian flow, but the presence of the dispersed phase, especially red blood cells due to their high-volume concentration and easily deformable property, give blood flow non-Newtonian characteristics. One of the well-known characteristics of blood flow is the axial accumulation of RBCs in a small vessel. Unlike smooth spherical particles, deformable particles such as RBCs are pushed away from the wall due to hydrodynamic interaction [1], leading to a change of the apparent viscosity depending on the size of the vessel [2,3].

Radmacher et al. [4] measured the elastic modulus of platelets using atomic force microscopy and reported that it ranged from 1–50 kPa. In the case of spheroidal rigid particles [5] or platelets, which are much less deformable than RBCs, the lift force acting on platelets is much smaller in the shear flow, and the interaction between RBCs and platelets causes the lateral migration of platelets, resulting in a nonhomogeneous distribution of platelets. Platelets play an important role in hemostasis, and they are activated by direct contact with damaged vessel walls [6]. Therefore, a high platelet concentration near the vessel wall is important in that it increases the chance of platelets contacting the injured vessel [7,8].

Plenty of experimental and numerical studies have been reported about the near-wall excess of the platelet concentration. Tangelder et al. [9] conducted an in vivo experiment showing that the platelet concentration near the arteriovenous wall was higher than that in the central area using intravital microscopy. Subsequently, in vitro experiments were conducted to determine the dependence on rheological parameters such as the shear rate, hematocrit (Ht), and particle size. Tilles and Eckstein [10] examined the effects of Ht and the shear rate on platelet migration by platelet-sized particles and found that platelet migration only occurs when the volume concentration of red blood cells is greater than 7%. It has also been observed that platelet margination occurs when the shear rate is above 430 s${}^{-1}$. Observations by the laser Doppler method [11] reported a similar trend to the previous results. They proposed the platelet transport theory based on convection–diffusion and calculated the platelet diffusivity from the power law function of hematocrit and the shear rate. Seki et al. [12] found that the lateral movement of platelet margination across the channel with Y-shaped bifurcated channels occurs in the Cell-Free Layer (CFL) region.

A number of numerical models and direct simulations have been conducted to elucidate the mechanism of platelet or platelet-sized particle margination. Eckstein and Belgacem [13] described the mathematical modeling of platelet transport as a combination of the diffusion and convection model, and a fundamental analogy was well reproduced by simple modeling. Almomani et al. [14] explored the interaction between RBCs and platelets by simulating the two particles as spherical particles. The results showed that the relative size of RBCs and platelets is one of the key factors that causes platelet margination. Crowl and Fogelson [15] performed a two-dimensional simulation using the lattice Boltzmann and immersed boundary method and calculated the lateral diffusion of platelets from their trajectories. They pointed out that the rate of diffusion was large in the core of the red blood cells and almost negligible in the area close to the wall. Mehrabadi et al. [16] simulated the flow of RBCs and platelets using the lattice-Boltzmann-spectrin-link (LB-SL) coupling method, drew a raw scaling for platelet margination in parallel plates, and reported that platelet diffusion is mainly enhanced by red blood cells. Several groups [17,18] have examined and reported the effects of erythrocyte size and deformability on platelet margination. Chang et al. [17] reported that when RBCs become less deformable due to diseases such as type 2 diabetes mellitus, there is less collision between red blood cells and platelets, resulting in less platelet margination. Similarly, Czara et al. [18] reported that the lesser deformability of RBCs leads to a decrease of the CFL, which in turn leads to a decrease of platelet margination.

Krüger [19] investigated the effect of tube diameter and capillary number on platelet margination and showed the reversibility of platelet trapping in the CFL, changing based on the capillary number $0.2$, resulting from the change of the motion of RBCs. There is also the aspect of the different deformability of the two types of capsules. Kumar and Graham [20] evaluated the effect of hydrodynamic forces and the collisions induced by the shear rate between particles of the same size, but with different deformability. Seki and Takinouchi [21] recently reported that in experiments using a 50 × 50 μm rectangular channel under Ht = 20% and $Ht$ = 40%, a non-uniform distribution of platelet-sized particles was observed. More particles were located near the corner of the channel, and this tendency was more pronounced at $Ht$ = 40%. Zhao et al. [22] reported that in a flow containing RBCs in a 10 × 15 μm rectangular channel, 3 μm beads were marginated to the corner of the wall. Coupling continuum mechanical simulation with the bio-chemical reaction of platelet has been reported. Ii et al. [23] presented the full Eulerian coupling method between the blood flow and hyperelastic capsule solving both the governing equation for the FSI problem and the stochastic differential equation for cell adhesion. In this present study, such a chemical reaction was not considered for the numerical simulation. For several years, platelet margination in pure shear flow [16,24] has been examined comprehensively. However, since blood vessels have a closed shape such as circular or rectangular, a detailed study of the motion of particles with RBCs in such a shape is necessary for practical purposes. As discussed by Seki and Takinouchi [21], understanding the migration behavior of particles in a square channel can be applicable, for example to cell sorting using microfluidic channels. They reported that in a square channel, more particles were located near the corner of the channel. In a circular tube [25], on the other hand, the particles are located along the circular walls, which shows the essential difference of their characteristics. Due to this difference, for example, in a flow with bifurcation, the particle distribution and margination behavior are different between circular and square channels. Unlike circular channels, few studies have been reported on platelet margination in a square channel.

In the present study, numerical simulations using the immersed boundary method were performed to reproduce the 3D hydrodynamic interaction between red blood cells and platelet-sized particles. The main objective of this study was to understand the effect of several parameters that affect the margination of particles in a square channel flow. Since most of the previous studies were based on pure shear flow, the present study focused on pressure-driven straight square channel flow with various shear rates, channel sizes, and hematocrit levels. The relationship between the hematocrit and the stable position of particles is also discussed.

2. Numerical Methods

In this section, we describe the numerical method employed in this study. Blood plasma was treated as a Newtonian fluid, and so can the fluid inside blood cells. The immersed boundary method [26] was used to solve the fluid–membrane interaction problem. RBCs and particles were modeled as vesicles having hyperelastic membranes.

2.1. Governing Equations for Fluid–Membrane Interaction Problems

The continuity and momentum equations of incompressible flow with the membrane force are given by:

$$\begin{array}{cc}& \nabla ·\mathbf{u}=\mathbf{0},\end{array}$$

(1)

$$\begin{array}{cc}& \rho \left({\displaystyle \frac{\partial \mathbf{u}}{\partial t}}+\left(\mathbf{u}·\nabla \right)\mathbf{u}\right)=-\nabla p+\nabla ·\left\{\eta \left(\nabla \mathbf{u}+\nabla {\mathbf{u}}^T\right)\right\}+\mathbf{f}\end{array}$$

(2)

where $\rho$ is the constant density, $\eta$ is the viscosity, p is the pressure, u is the velocity vectors, and f is the membrane force vectors applied from the elastic stress of the cell with the Immersed Boundary Method (IBM) [26]. The difference of the viscosity between the inside and the outside of the cells was treated in the same way as Tryggvason’s front-tracking method [27]. Let us denote by $\varphi (\mathsf{x},t)$ a smoothed phase indicator function, and it assigns the fluid domain as 0, the inside domain as 1, and $\in (0,1)$ for the transition layer. Then, the viscosity is calculated from:

$$\eta =\varphi {\eta }_i+\left(1-\varphi \right){\eta }_f$$

(3)

where ${\eta }_i$ and ${\eta }_f$ are the viscosities of the inside capsule domain and fluid. Note that we did not consider the nonlinear convection term $(\mathbf{u}·\nabla )\mathbf{u}$ from (2). Since the Reynolds number ($Re$) was quite low ($Re<1$) in our simulation, the inertia effect was negligible.

RBC membranes are characterized by two major elastic contributions: in-plane shear and bending resistance. In this study, we adopted the Skalak model [28] as the in-plane stress ${\mathbf{T}}_s$ of red blood cells:

$${\mathbf{T}}_s={\displaystyle \frac{G_s}{J_s}}\left(\left(\mathrm{tr}\left({\mathbf{B}}_s\right)-1\right){\mathbf{B}}_s+J_s^2\left(C\left(J_s^2-1\right)-1\right)\mathbf{P}\right)$$

(4)

where $G_s$ is the shear modulus, C is the area dilation modulus, ${\mathbf{B}}_s$ is surface left Cauchy–Green tensor, $J_s$ is the surface Jacobian, the second invariant of ${\mathbf{B}}_s$, and $\mathbf{P}=\mathbf{I}-\mathbf{nn}$ is the surface projection tensor. We employed the Pozrikidis model [29] for the bending modulus ${\mathbf{T}}_b$,

$$\begin{array}{cc}& {\mathbf{T}}_b=\mathbf{qn},\end{array}$$

(5)

$$\begin{array}{cc}& \mathbf{q}=\left({\nabla }_s·\mathbf{m}\right)·\mathbf{P},\end{array}$$

(6)

$$\begin{array}{cc}& \mathbf{m}=E_b(\mathbf{\kappa }-{\kappa }_0\mathbf{P})\end{array}$$

(7)

where $E_b$ is the bending modulus and ${\kappa }_0$ is the spontaneous curvature. In this study, we employed the infinitesimal strain assumption, namely $\mathbf{\kappa }=\kappa \mathbf{P}$, where $\kappa ={\nabla }_s·\mathbf{n}$ is the curvature. The static local equilibrium equation of the membrane is written as:

$${\nabla }_s·\left({\mathbf{T}}_s+{\mathbf{T}}_b\right)+\mathbf{F}=0$$

(8)

The Finite-Element Method (FEM) [30] was used to solve (8) and obtain the elastic force vector F acting on the Lagrangian points of the membrane. The membrane force f from the exerted force on the Lagrangian points X $\in \Gamma$ is given by:

$$\mathbf{f}\left(\mathbf{x}\right)=\sum _{\mathbf{X}\in \Gamma }{D}_h^{\left(3\right)}\left(\mathbf{x}-\mathbf{X}\right)\mathbf{F}\left(\mathbf{X}\right)$$

(9)

where $D_h^{\left(3\right)}\left(\mathbf{x}\right)$ is the three-dimensional approximate Dirac delta function,

$$D_h^{\left(3\right)}\left(\mathbf{x}\right)={\displaystyle \frac{1}{h^3}}{\delta }_h\left({\displaystyle \frac{x}{h}}\right){\delta }_h\left({\displaystyle \frac{y}{h}}\right){\delta }_h\left({\displaystyle \frac{z}{h}}\right)$$

(10)

where h is the computational mesh size and ${\delta }_h\left(s\right)$ is the continuous function,

$${\delta }_h\left(s\right)=\left\{\begin{array}{cc}\frac{1}{8}\left(3+2\left|s\right|+\sqrt{1+4\left|s\right|-{4\left|s\right|}^2}\right),\hfill & \mathrm{for}0\le \left|s\right|<1,\hfill \\ \frac{1}{8}\left(5-2\left|s\right|+\sqrt{-7+12\left|s\right|-{4\left|s\right|}^2}\right),\hfill & \mathrm{for}1\le \left|s\right|<2,\hfill \\ 0,\hfill & \mathrm{for}2\le \left|s\right|\hfill \end{array}\right.$$

(11)

The location of Lagrangian points $\mathbf{X}$ was updated following the IBM,

$$\frac{d\mathbf{X}}{dt}=\sum _{\mathbf{x}\in \Omega }{D}_h^{\left(3\right)}\left(\mathbf{x}-\mathbf{X}\right)\mathbf{u}\left(\mathbf{x}\right)$$

(12)

2.2. Discretization

The Simplified-Marker-And-Cell method (SMAC method) [31] was used to couple the velocity and pressure. In the prediction phase, the tentative velocity ${\mathbf{u}}^p$ was updated by:

$$\frac{{\mathbf{u}}^p-{\mathbf{u}}^n}{\Delta t}+\frac{\nabla p^n}{\rho }-\frac{1}{\rho }\nabla ·\left\{{\eta }^n\left(\nabla {\mathbf{u}}^p+{\left(\nabla {\mathbf{u}}^p\right)}^{\mathbf{T}}\right)\right\}+\frac{{\mathbf{f}}^n}{\rho }=0$$

(13)

where the superscript n represents a discrete time step and $\Delta t$ is a discrete time interval. The pressure term was decided to satisfy the continuity equation of $n+1$ step, described as:

$$\nabla ·{\mathbf{u}}^{n+1}=\nabla ·\left({\mathbf{u}}^p-\Delta t\nabla \varphi \right)=0$$

(14)

where $\varphi$ is the scalar potential. Equation (14) could be transformed in to a Poisson equation of $\varphi$:

$${\nabla }^2\varphi =\frac{1}{\Delta t}\nabla ·{\mathbf{u}}^p$$

(15)

Therefore, the pressure was determined by $p^{n+1}=p^n+\rho \varphi$, and in the projection phase, the velocity was updated as:

$${\mathbf{u}}^{n+1}={\mathbf{u}}^p-\Delta t\nabla \varphi$$

(16)

Note that the second-order central difference method was employed for the spatial discretization in (1) and (2), and a staggered grid was used for the numerical calculation. For a more detailed description and a numerical validation, see [32].

2.3. Simulation Conditions

Cellular flows consisting of plasma, RBCs, and platelets with hematocrit ($Ht$) ranging from $0\%$ to $40\%$ were numerically solved in this study. The computational domain was a square channel that was periodic in the x-direction and wall-bounded in the y- and z-directions. The pressure-driven body force was exerted in the flow direction x. The height H and width W of the channel domain were the same, but the specific values were different depending on the simulation case. The length of the domain L was chosen to 7-times longer than the lateral direction. The shape of the undeformed RBC was biconcave with a diameter of $7.92$ μm. The viscosity of the inner RBC was set to 5-times higher than that of the plasma, and this is similar to the physiological condition. The shear modulus and area dilation modulus of the RBC were $G_s=5.0$ μN/m and $C=50$, respectively. These values are almost identical to those of Takeishi’s numerical experiments [33], which successfully reproduced the deformation of red blood cells in shear flow. The bending modulus $E_b=1.8\times {10}^{-19}$ N·m was also considered. Note that all the mechanical properties of the RBC were obtained experimentally [34,35,36], except for C; the area dilation modulus C was adopted as a reasonable value to prevent the surface area of the RBC from changing. The particle was modeled as sphere-shaped with a diameter of $2.9$ μm. To reproduce the less deformable property of platelets, the membrane shear elasticity and area dilation modulus of the particle were considered 10-times higher than those of the RBC. It is known that the volume fraction of platelets is 0.1–0.2% under physiological condition [37]. In the present simulation, the volume fraction of platelets was set to $0.5\%$ in order to have enough samples for statistical analysis. In the initial state, RBCs and particles were regularly arranged with little random variation. Here, we adopted time step size $\Delta t=2.5$ (μs) and mesh size $\Delta x$ = $D/32$, where D is the diameter of the RBC. See Appendix A for the effect of the grid size on platelet margination.

2.4. Evaluation of Particle Migration toward the Near-Wall Region

It is known that platelet or platelet-sized particles move toward the near-wall region with the hematocrit level over $7\%$ in blood flow [10]. In the current simulation, the initial shape of RBCs was modeled as biconcave. It took time to fully develop and move toward the center to produce the cell-free layer. Figure 1 illustrates how we analyzed the simulation results.

In order to minimize the influence of the initial configuration, we quantified the time history of particle motion and evaluated the change of the ensemble-averaged distance between the particle and the centerline of the computational domain. All particles that passed through a square region of width $H/4$ and height $H/4$ near the center were counted as samples for the ensemble average. Since platelets eventually migrate near the corners, the path of platelets from the central region to the channel wall can be investigated by calculating the trajectory of platelets passing through the center of the rectangular channel in each lateral direction. One example of comparing the time history of the path of platelets in each wall direction is shown in Figure 2.

3. Results

Several studies have been reported on how parameters such as the shear rate and channel size affect the platelet margination in parallel plates [16,22] or circular tubes [25]. As we mentioned in the Introduction, Zhao et al. [22] reported that $3$ μm spherical beads exhibit a high concentration at the corner of a rectangular channel of $10\times 15$ μm${}^2$ in cross-section in RBC flow. In the present study, we investigated how parameters such as the shear rate, channel width, and hematocrit value affected the temporal evolution of particle migration in a square channel.

3.1. Direction Dependency of Particle Migration

Figure 2 shows the time history of the ensemble-averaged distance $\Delta \overline{d}\left(t\right)$ between a particle and the centerline of the straight channel for $Ht=20\%$, $\dot{{\gamma }_w}=972$ s${}^{-1}$ and $L=224$ μm, $H=32$ μm, along with the error bars indicating the standard deviation. $\Delta \overline{d_y}\left(t\right)$ and $\Delta \overline{d_z}\left(t\right)$ represent the ensemble-averaged distances in the height and width direction, respectively. For example, if particles are entering the central square region of width H/4 and height H/4 and their y-position at time t is $d_y\left(t\right)$, the ensemble-averaged distances are described by the following equation,

$$\Delta \overline{d_y}\left(t\right)=\left|\overline{d_y}\left(t\right)-H/2\right|$$

(17)

Although each case showed a large standard deviation, the time history of ensemble-averaged distances showed a similar pattern, suggesting that the overall path for how the particles migrate in the square duct was almost identical for each lateral direction. Hereafter, $\Delta \overline{d}\left(t\right)$ was calculated from $\Delta \overline{d}\left(t\right)=\left(\Delta \overline{d_y}\left(t\right)+\Delta \overline{d_z}\left(t\right)\right)/2$. Since Figure 2 shows that the average path of platelets migrating in both wall directions was almost equal, we present the results such as the particle or RBC concentration profiles in one direction(y) to discuss the physical phenomenon from now on.

Figure 3 shows the snapshots of the present simulation results. In order to focus on the lateral movement of platelets and RBCs, the results of the 3D simulation were visualized by projecting them onto the yz plane (two-dimensional) using the volume fraction information of RBCs and platelets on a fixed grid. Therefore, the shape of the red circle at t = 0 corresponds to RBCs and the shape of the white sphere corresponds to platelets. The simulation conditions corresponded to RBCs and particles in a square channel flow under the conditions of $Ht=20\%$, $\dot{{\gamma }_w}=972$ s${}^{-1}$ at different times of $t=0$, $0.5$ s, $1.0$ s, $1.5$ s, $2.0$ s, and $3.0$ s. The red cloud near the center after t = 0.5 s corresponds to the RBCs. The intensity of the color means the probability of the presence of RBCs or platelets. RBCs gradually move toward the center from the initial position and generate the cell-free layer. All the particles move to the corner with a stochastic fluctuation. Figure 4 shows the change in the particle concentration profile, and each graph corresponds to the result of Figure 3. As we can confirm from Figure 4b–e, many particles did not move directly to the corner, rather migrating toward the wall first, then approached the corner while moving along the wall. After 3 s, most of the particles migrated to the corner while moving 11 mm on average in the flow direction, equivalent to about 340-times the channel width H. Figure 5 shows the change in RBCs’ concentration profile. The axial accumulation of RBCs occurred within 0.1 s, a much faster time-scale than platelet margination.

3.2. Shear Rate Dependency

We investigated the effect of the shear rate on particle margination. The computational domain was set to $L=224$ μm, $H=32$ μm, and hematocrit $Ht=20\%$. For a square duct, the average wall shear rate $\dot{{\gamma }_w}$ was calculated from the following equation,

$$\dot{{\gamma }_w}=\frac{U\lambda }{2H}$$

(18)

where U is the average velocity and $\lambda$ is the shape factor, known as $\lambda =14.3$ for a square channel of $Newtonian$ fluid [38]. The simulation conditions are listed in

Table 1. Conditions for numerical simulation to compare the effect of the shear rate.

U (m/s)	$\mathbf{Ca}$	$\mathbf{Re}$	$\dot{{\mathit{\gamma }}_{\mathit{w}}}$ (s${}^{-1}$)
0.0044	0.66	0.116	972
0.0059	0.90	0.159	1329
0.0089	1.35	0.238	1990

. Note that the average velocity U was calculated from the simulation results. The Reynolds number $Re$ and capillary number $Ca$ were calculated from $Re=\rho UH/\eta$ and $Ca=\eta \dot{{\gamma }_w}{a}_{\mathrm{RBC}}/G_s$ where $a_{\mathrm{RBC}}={(V^3/4\pi )}^{1/3}$ is the equivalent diameter of the RBC. In the present simulation, $a_{RBC}=2.82$ μm and $V=93.58$ μm${}^3$ were used as the volume of the RBC.

Figure 6a shows the time history of the average distance from the centerline $\Delta \overline{d}$, and Figure 6b represents the change of $\Delta \overline{d}$ as a function of the distance traveled in the flow direction $\overline{x}$. As can be seen from Figure 6a, particles moved faster toward the wall direction at a higher shear rate.

Comparing $\Delta \overline{d}$ as a function of distance traveled toward the streamwise direction $\overline{x}$, the path taken by the particles to approach the wall from the central region was almost the same at $\dot{{\gamma }_w}=972$ s${}^{-1}$, 1329 s${}^{-1}$, and 1990 s${}^{-1}$. This result is in good agreement with the reference results in parallel plates or bifurcations [12,16]. Usually, the CFL thickness increases as the shear rate increases [39]. However, with the parameter range of the current simulations, the CFL thickness was almost constant, and the value of $\Delta \overline{d}$ was the same when they reached a stable position near the wall, as shown in Figure 6a.

3.3. Channel Size Dependency

We examined the dependency of particle margination on the channel size. We varied the size of the square channel $H=16$ μm, 24 μm and 32 μm and adjusted the average velocity U for each case, so that the wall shear rate $\dot{{\gamma }_w}$ was equal to 1329 s${}^{-1}$, respectively. As we mentioned, in (21), the diffusion of particles depends on factors such as the hematocrit level, the wall shear rate, and the particle size. Since the simulation conditions were the same except the channel width $H=16$ μm, 24 μm, and 32 μm, the time history of $\overline{d}$ overlapped in RBC-core region from 0 s to 0.5 s. With reference to the derivation of Mehrabadi et al. [16], we combined (18), (20), and (21) to derive the following equation.

$${\left(\frac{\Delta \overline{d}}{d_r}\right)}^2=k\lambda \frac{\overline{x}{a}^2}{H{d_r}^2}$$

(19)

Note that t is replaced by $t=\overline{x}/U$ and $d_r$ stands for the distance from the center of the channel to a stable position of the particle, which had different values in each case. As shown in Figure 7b, $\Delta \overline{d}$ normalized by $d_r$ was well scaled by $\overline{x}$ normalized by $a^2/H{d_r}^2$. The disagreement between the curves can be attributed to the fact that the aspect ratio of the platelet particles to the rectangular tube was different in each case, resulting in different hydrodynamic interactions with the wall and platelets, and also, the local distribution of hematocrit was not totally the same.

3.4. Hematocrit Dependency

Numerical simulations were performed for hematocrit $Ht=0\%$, $5\%$, $10\%$, $20\%$, $30\%$ and $40\%$ to explore the relationship between hematocrit and its effect on particle margination. The wall shear rate $\dot{{\gamma }_w}$ in the plasma was 1910 s${}^{-1}$, and the same pressure gradient was applied in all cases. The simulated results were affected by many complicated factors, and it is not easy to discuss suitable dimensionless parameters. As shown in Figure 8c,d, the stable position $d_r$ for a particle, defined as the distance between the centerline and the center of the particle at equilibrium, varied with the hematocrit level. Another aspect is that rate of $\Delta \overline{d}$ as a function of $\overline{x}$ was not constant because the margination rate of the particles varied with the distance from the centerline. This was because the local shear rate and local hematocrit were different from the location in the plane. If we focus on the case of hematocrit = $20\%$, the margination rate of $\Delta \overline{d}$ was nearly constant up to $\Delta \overline{d}=12$ μm, after which the margination rate slowed down. We confirmed that this region of the constant margination rate of $\Delta \overline{d}$ was consistent with the RBC-core region.

Despite the difficulties in the analysis described earlier, several points can be drawn. Figure 8a shows that in the absence of RBCs ($Ht=0\%$), particle margination did not occur. This result confirmed that the presence of RBCs was the main factor that caused the non-uniform particle concentration in the domain. We confirmed from the simulation result in the case of $Ht=5\%$ that the particles did not reach the corners, but stayed in a circular shape near the boundary between the RBC-core region and the CFL. In other cases, over $Ht=10\%$, the particles gradually migrated toward the wall. Once they reached the wall, the particles moved along the wall to near the corner of the channel. Figure 8b shows the ratio of $\overline{d}$ normalized by the stable position $d_r$ for each case as a function of $\overline{x}$. If we focus on the region where the ratio of $\overline{d}/d_r$ was constant, we can see that the ratio of $\overline{d}/d_r$ reached its peak at $Ht$ = $10\%$ and gradually decreased as the hematocrit increased. Figure 9, Figure 10 and Figure 11 show the results of platelet margination in a flow with $Ht$ = 40% for $\dot{{\gamma }_w}=712$ s${}^{-1}$. Figure 8c,d indicates the relationship between hematocrit and the width of the CFL and relationship between the width of the CFL and the stable position $d_r$. There are several ways to define the thickness of the CFL [16,19]. In the present study, the width of the CFL ${\delta }_{CFL}$ was defined as the distance at which the local hematocrit at equilibrium was less than 0.5% from the wall, following the definition from Mehrabadi et al. [16]. We confirmed that the distance of the stable position of the particles decreased in proportion to the width of the CFL. This can be understood if we consider that the margination of the particle was attributed to the presence of RBCs and the width of the CFL was determined by the hematocrit.

Figure 9f and Figure 10f indicate that many particles could not escape from the RBC-core region and stayed in the region in 5 s. Particles moved 12 mm on average in the flow direction during 5 s, almost the same distance as Figure 3f, but many more particles still remained in the RBC-core region. This can be explained by the fact that particle margination was caused not only by the spatial gradient of the velocity fluctuations resulting from the presence of the RBCs and their axial accumulations, but by the direct interaction between the RBCs and the particles, which pushed the particles in the wall direction. Due to their volume effect, RBCs that were excessively present in the RBC-core layer may act as an obstacle and prevent the particles from moving. Vahidkhah et al. [40] reported that the motion of platelets varies with the density of RBCs present near a platelet. For example, a platelet shows discontinuous and quick motion when very few RBCs are located around the platelet. Such a movement rarely occurs when the hematocrit value is high. Figure 11 shows that RBC concentration profiles in $Ht$ = $40\%$. Similar to $Ht$ = 20%, the axial accumulation of RBCs took place on a much faster time scale than the particles’ margination, within 0.2 s. We can confirm the results of the present simulations through the comparison with the experimental data by Seki and Takinouchi [21]. As mentioned earlier, their experiments were conducted in a rectangular tube with a cross-section of 50 × 50 μm, n average wall shear rate of 950 s${}^{-1}$, a particle size of 2.9 μm, and $Ht$ = 20% and $Ht$ = 40%, which are almost the same conditions as the current simulations except for the channel size. Compared to $Ht$ = 20%, platelet-sized particles moved closer to the corners, and more platelets were observed near the center region at $Ht$ = 40%, which is consistent with the results of the current simulations.

4. Discussion

As for the platelet margination, the platelets moved from the region with a large velocity fluctuation to the region with a small velocity fluctuation, which indicates the corner of the square channel. Here, we define the velocity fluctuations as the squared mean in-plane velocity in the lateral direction. The motion of platelets in the RBC-core region can be treated as a diffusion process. Many studies [8,41,42] compared how fast platelets move by calculating the diffusion coefficient. Among the components of the diffusion tensor, we calculated the diffusion coefficient $D_{yy}$ from the following equation,

$${(\Delta \overline{d})}^2=2D_{yy}t$$

(20)

Since particles approached the wall direction and stayed in the corner of the square domain in the end, $\Delta \overline{d}$ was not constant, so nor was ${(\Delta \overline{d})}^2$, as we can confirm from Figure 12a. Therefore, we determined $D_{yy}$ from the region where the change of ${(\Delta \overline{d})}^2$ was constant. The result is plotted in Figure 12b, and the average wall shear rate $\dot{{\gamma }_w}$ had a linear relationship with the diffusion coefficient $D_{yy}$. This result was supported by finding that the particle diffusivity is expressed as:

$$D_{yy}=k\dot{\gamma }{d}^2$$

(21)

where k is a dimensionless parameter related to hematocrit, $\dot{\gamma }$ is the shear rate, and d is the particle diameter. Seki et al. [12] estimated the diffusivity of $2.9$ μm-diameter particles to be $40$ μm${}^2$/s at a shear rate of 1000 s${}^{-1}$, and in the present simulation, the diffusivity of the particle $D_{yy}=30.12$ μm${}^2$/s at $\dot{{\gamma }_w}=1000$ s${}^{-1}$, which is a comparable value.

To explore the effect of velocity fluctuations on particle migration, velocity fluctuations were calculated from the time-averaged and spatially averaged velocity components of the plasma in $Ht$ = 20% at $\dot{{\gamma }_w}=1990$ s${}^{-1}$. The profile of the velocity fluctuations is shown in Figure 13. We confirmed that the lateral velocity component was largest near the boundary between the red blood cell layer and the CFL and smallest near the wall. It can also be confirmed that each velocity component in the transverse direction was high near the wall due to the fluctuations generated by the particle movement. In the present simulations, the Reynolds number was much smaller than one, which corresponds to laminar flow. However, it showed the so-call pseudo-turbulence effect caused by the motion of the disperse phase, which was RBCs in this study. Seeing Figure 3 and Figure 13, one can easily find that platelets gradually moved from a region of high velocity fluctuation to a region of low velocity fluctuation, that is the corners near the channel. This phenomenon is analogous to thermophoresis [43]: large molecules accumulate in a lower-temperature region, which corresponds to the lower molecular fluctuation. This is expressed by the so-called thermophoretic force, which depends on a temperature gradient. A similar example is turbophoresis [44,45], which is particle migration toward a lower-velocity-fluctuation region in turbulent flow. Therefore, platelet migration in the present simulations was considered as turbophoresis in pseudo-turbulence flow.

Figure 14 shows the velocity fluctuations for different shear rates in the z = 16 μm plane. The velocity distribution normalized by the wall shear rate collapsed into one curve, indicating that the amplitude of velocity fluctuation has a linear relationship with the shear rate. As explained in Section 3.2, the particle trajectories did not depend on the shear rate. Since the shear rate was proportional to the streamwise flow velocity, this indicates that the lateral force was proportional to the streamwise flow velocity. Overall, platelet margination was mainly caused by the spatial gradient of velocity fluctuations, which was generated by the motion of the RBCs and was proportional to the flow velocity.

5. Conclusions

In this study, we investigated numerically the effects of physical parameters such as the shear rate, channel size, and hematocrit level on platelet-sized particles’ migration toward the wall in a square channel. The simulation results showed that the migration of particles becomes faster as the shear rate increases. The results also showed that the margination rate—the ratio of the distance traveled in the flow direction to the margination distance in the wall direction—does not depend on the shear rate. This indicates that the lateral force acting on the particles increases linearly with the shear rate, directly related to the velocity in a square channel.

Particle margination is mainly caused by two effects: first, the force pushing the particles directly toward the wall because of the interaction between the RBCs and the particles; second, the velocity fluctuations generated by the flowing motion of the deformable red blood cells. The velocity fluctuations are the highest near the border between the red blood celllayer (RBC-layer) and the cell-free layer (CFL) and the lowest around the corner. Due to the velocity fluctuation differences, the particles gradually move from the RBC-layer, with large fluctuations, to the vicinity of the corner, with small fluctuations, and stay there. Moreover, the wall-shear-rate-normalized velocity fluctuations collapse into one curve as well, showing that the margination rate does not depend on the shear rate. This result implies that the particle margination is mainly dependent on the velocity fluctuation differences.

In a square channel flow, we investigated the relationship between the width of the square channel and the path of platelet migration. We found that the margination distance normalized by the particle size, shear rate, and effective area of the particle movement showed good agreement independent of the channel size.

We found that the margination distance normalized by the particle size, shear rate, and effective area of the particle movementshowed a similar trend regardless of the channel size.

Finally, we confirmed that the margination rate of particles varies nonlinearly with the hematocrit value. In this study, we observed that under 0–40% hematocrit, the margination rate was the highest at 10% of hematocrit and gradually decreased beyond 10% due to the presence of red blood cells, which inhibited the movement of the particles toward the wall.