The nature of blood as a suspension of red blood cells makes computational haemodynamics a demanding task. Our coarse-grained blood model, which builds on a lattice Boltzmann method for soft particle suspensions, enables the study of the collective behaviour of the order of 106 cells in suspension. After demonstrating the viscosity measurement in Kolmogorov flow, we focus on the statistical analysis of the cell orientation and rotation in Couette flow. We quantify the average inclination with respect to the flow and the nematic order as a function of shear rate and haematocrit. We further record the distribution of rotation periods around the vorticity direction and find a pronounced peak in the vicinity of the theoretical value for free model cells, even though cell–cell interactions manifest themselves in a substantial width of the distribution.
Human blood can be approximated as a suspension of deformable red blood cells (RBCs; also called erythrocytes) in a Newtonian liquid, the blood plasma. The other constituents like leukocytes and thrombocytes can be neglected owing to their low volume concentrations. Typical physiological RBC concentrations are around 50 per cent. In the absence of external stresses, erythrocytes assume the shape of biconcave discs of approximately 8 μm diameter . An understanding of their effect on the rheology and the clotting behaviour of blood is necessary for the study of pathological deviations in the body and the design of microfluidic devices for improved blood analysis.
Well-established methods for the computer simulation of blood flow either consist of an elaborate model of deformable cells [2,3] or restrict themselves to a continuous description at larger scales . Our motivation is to bridge the gap between both classes of models by an intermediate approach: we keep the particulate nature of blood but simplify the description of each cell as far as possible. The main idea of the model  is to distinguish between the long-range hydrodynamic coupling of cells and their complex short-range interactions. Our method is well suited for the implementation of complex boundary conditions and an efficient parallelization on parallel supercomputers. This is important for the accumulation of statistically relevant data to link bulk properties, for example, the effective viscosity, to phenomena at the level of single erythrocytes. In this paper, we present an alternative method of viscosity measurement and a study of the orientational dynamics of cells in sheared suspensions. The improved understanding of the dynamic behaviour of blood might be used for the optimization of macroscopic simulation methods.
2. Coarse-grained model for blood-flow simulations
We apply a D3Q19/Bhatnagar–Gross–Krook (BGK) lattice Boltzmann method for modelling the blood plasma . The single-particle distribution function nr(x,t) resembles the fluid travelling with one of r=1,…,19 discrete velocities cr at the three-dimensional lattice position x and discrete time t. Its evolution in time is determined by the lattice Boltzmann equation 2.1 being the BGK-collision term with a single relaxation time τ. The equilibrium distribution function is an expansion of the Maxwell–Boltzmann distribution. and can be identified as density and momentum. In the limit of small velocities and lattice spacings, the Navier–Stokes equations are recovered with a kinematic viscosity of ν=(2τ−1)/6, where τ=1 in this study.
For a coarse-grained description of the hydrodynamic interaction of cells and blood plasma, a method similar to the one by Aidun et al. modelling rigid particles of finite size is applied . It can be seen as a bounce-back boundary condition, which is enhanced by a correction term C=2αcrρ(x+cr,t)crv/c2s that accounts for a possible local boundary velocity v. The lattice weights αcr and the speed of sound cs are constants for a given set of discrete velocities. The resulting bounce-back rule 2.2 and specifying the direction opposite to r is applied to all fluid distributions that, according to equation (2.1), would travel along a link that crosses the theoretical cell surface. The momentum, , which is transferred during each time step by each single bounce-back process, is used to calculate the resulting force on the boundary.
Instead of the biconcave equilibrium shape of physiological RBCs, we choose a simplified ellipsoidal geometry that is defined by two distinct half-axes R∥ and R⊥ parallel and perpendicular to the unit vector , which points along the direction of the axis of rotational symmetry of each cell i. Since the resulting volumes are rigid, we allow them to overlap in order to account for the deformability of real erythrocytes. We assume a pair of mutual forces at each cell–cell link. This is exactly the momentum transfer during one time step owing to the rigid-particle algorithm for a resting particle and an adjacent site with resting fluid at equilibrium and initial density . While these surface forces solely serve to compensate the lack of fluid pressure inside the particles, we additionally implement phenomenological pair potentials for neighbouring cells that model the volume exclusion of soft and anisotropic particles in an effective way. For the sake of simplicity, we use the repulsive branch of a Hookean spring potential 2.3 for the scalar displacement rij of two cells i and j. With respect to the disc shape of RBCs, we follow the approach of Berne & Pechukas  and choose the energy and range parameters 2.4 and 2.5 as functions of the orientations and of the cells and their normalized centre displacement . We achieve an anisotropic potential with a zero-energy surface that is approximately that of ellipsoidal discs. Their half-axes parallel σ∥ and perpendicular σ⊥ to the symmetry axis enter equations (2.4) and (2.5) via and , whereas determines the potential strength.
Figure 1 shows an outline of the model. Two cells i and j surrounded by blood plasma are displayed. For clarity, the three-dimensional model is visualized in a two-dimensional cut. Depicted are the cell shapes defined by the zero-energy surface of the cell–cell potential equation (2.3) with equation (2.4) and equation (2.5) that can be approximated by ellipsoids with the size parameters σ∥ and σ⊥ as half-axes. The cells are free to assume continuous positions and orientations oi and oj. As a consequence, the centre displacement vectors rij between cells are also continuous. Still, for the cell–plasma interaction, an ellipsoidal volume with half-axes R∥ and R⊥ is discretized on the underlying lattice. For additional clarity, the discretization is drawn only for cell j. The forces and torques emerging from the interaction of the cells with other RBCs and the fluid are integrated by a classical molecular dynamics code in order to evolve the system in time. The conversion from lattice units to physical units is done by multiplication with δx=2/3 μm, δt=6.80×10−8 s and δm=3.05×10−16 kg for length, time and mass, respectively. As a convention, primed variables are used whenever we refer to quantities specified in physical units. The model parameters are chosen as R′⊥=11/3 μm, R′∥=11/9 μm, σ′⊥=4 μm, σ′∥=4/3 μm and . For more detailed information see Janoschek et al. .
(a) Apparent bulk viscosity measurement in Kolmogorov flow
In our previous work, the apparent viscosity of the bulk of the suspension was measured in simulations of unbounded Couette flow . Though this setting is easy to comprise analytically, its efficient and precise implementation proves a non-trivial challenge since the suspended particles would have to be enabled to stretch across the sheared boundary. We therefore demonstrate an alternative method to determine the apparent viscosity that is compatible with completely periodic boundaries and is based on the so-called Kolmogorov flow: a sinusoidally modulated shear flow. We proceed as in Benzi et al.  and apply to the full suspension a sinusoidal body force 3.1 pointing into the z-direction. f0 is the amplitude and k the wavenumber in the x-direction. At steady state, the spatial variation of the shear stress 3.2 matches the external forcing. Integration of equation (3.2), together with equation (3.1), yields an analytical expression for the shear stress 3.3 The y- and z-averaged local shear rate can be evaluated numerically as 3.4 After equilibration, a simulation of Nx×Ny×Nz lattice sites results in Nx numbers for the apparent viscosity 3.5 covering a range of shear rates that depends on f0 and k.
Compared with a viscosity measurement in Couette flow , the above procedure  has the benefit of taking advantage of the whole simulation volume since there are no possibly unphysical boundary regions. Additionally, each single measurement yields data for many different shear rates. On the other side, the non-constant shear stress causes inhomogeneities in the local cell volume concentration Φ(x), which are much more pronounced than in the case of Couette flow. Therefore, we disregard all viscosity data (x,μapp(x)) for which Φ(x)∉[Φ*−ΔΦ,Φ*+ΔΦ] with Φ* being the volume concentration of interest and 2ΔΦ a small interval of tolerance. Figure 2 compares resulting from such measurements at varying f0 and Nx with k=2π/Nx after 68.0 ms with Couette flow measurements as described in Janoschek et al.  (Nx=128). In the plot, Φ*=43% and ΔΦ=0.5%. We further average the data within bins with a width of in order to reduce statistical noise. The good agreement proves the feasibility of the method.
(b) Statistical analysis of cell orientations in Couette flow
Owing to its applicability to huge numbers of particles and the fact that cell orientations are accessible directly in the model, our simulation method is particularly suited for the statistical analysis of RBC orientations at high volume concentrations. In the following, we simulate Couette flow as described in Janoschek et al. . However, with a size of Nx=Nz=512 and Ny=64 lattice sites, the volume is considerably larger and contains about 3×104 cells at a haematocrit of Φ=45%. The velocity gradient is aligned with the x-axis and the flow with the z-axis. To exclude boundary effects in the x-direction, we restrict ourselves to the analysis of those cells that stay at a lateral position 512/3<x<2×512/3 during the whole simulation. The visualization of the cell volumes defined by the model potential is shown for a smaller system in the inset of figure 3. It is clear that also at physiological volume concentrations, there is a preferential orientation of cells. In order to quantify this observation, the orientation vector of every cell i is multiplied with −1 where necessary to achieve . The angle θ is measured between the positive x-axis and the resulting normalized orientation vector. Figure 3 displays the most probable value θ* that is determined by fitting a Gaussian to the distribution of angles. For Φ=45%, compared with θ=0, which would mean alignment parallel to the flow, there is always a positive inclination θ* that is directed against the vorticity of the shear (see inset of figure 3). When increasing the shear rate from around 300 to 6000 s−1, this inclination becomes considerably smaller: θ* decreases from 14.5° to 3.4°.
Studying θ allows us to quantify the orientation with respect to the direction of the flow. To quantify the actual amount of orientational order, the nematic-order parameter known from liquid crystal physics (e.g. ) is a better choice. This parameter is usually defined as the largest eigenvalue λ+ of the nematic-order tensor , which is obtained as the average over different time steps t and cells i within the volume of interest. Possible values for λ+ are between 0 and 1, indicating complete disorder and order. Figure 4 depicts λ+ for different shear rates. A decrease of nematic order with increasing shear rate is visible. The inset shows the dependency on the haematocrit Φ at a fixed shear rate of . For comparison, λ+ resulting for a single ellipsoidal particle that tumbles with perpendicular to the vorticity direction is also shown. Apparently, the higher nematic order at Φ=45% is caused by hindered tumbling motion, while the reason for the lower λ+ at Φ=11% is the lack of alignment in the vorticity plane.
(c) Statistical analysis of cell rotations in Couette flow
The above conclusions make clear that preferential orientations result from the averaged tumbling of many cells. The time evolution of the continuous tumbling angle θ is plotted for an arbitrary selection of cells in the inset of figure 5. The respective volume concentration is Φ=45% and the shear rate . The cells keep a largely constant alignment for varying periods of time and, occasionally, flip over by an angle of π along the vorticity direction. Thus, the angular velocity is strongly time dependent. The probability distribution function of the time required for a rotation by π, which is measured as the time T between two flipping events, is shown in the main part of figure 5. For decreasing volume concentrations, the distribution becomes narrower and its peak is shifted towards shorter times. However, with a width comparable to its average value, even at Φ=11%, the distribution deviates substantially from a delta peak. Therefore, even though typical tumbling periods in suspension do not differ much from the value for a freely tumbling ellipsoid , the effect of cell–cell interactions on our observables is significant for all Φ studied.
Our mesoscale model for blood  was applied to bulk flow of up to 3×104 cells. We demonstrated the feasibility of viscosity measurement in Kolmogorov flow and studied the statistics of the rotational cell positions and velocities at different values of the haematocrit. These observables are extendable also to more elaborate (e.g. deformable) cell models. Their quantitative study and measurement are key to building more reliable cell-based or continuum descriptions for blood.
The authors acknowledge fruitful discussions with T. Krüger, financial support from the TU/e High Potential Research Programme as well as computing resources from JSC Jülich, SSC Karlsruhe, CSC Espoo, EPCC Edinburgh and SARA Amsterdam, the latter three being granted by DEISA as part of the DECI-5 project. J.H. acknowledges The Netherlands Organization for Scientific Research (VIDI grant no 10787).
One contribution of 26 to a Theme Issue ‘Discrete simulation of fluid dynamics: methods’.
- This journal is © 2011 The Royal Society