Shear-Induced Transport of Platelets in a Red Blood Cell Suspension

(1)

Shear-Induced Transport of Platelets in

a Red Blood Cell Suspension

Ziggy Pleunis (10001615)

August 14, 2014

Bachelor project Physics & Astronomy, 18 EC, research conducted between January and June 2014.

Supervisor : dr. Eric Lorenz Daily supervisor : drs. Lampros Mountrakis Second examiner : dr. Rudolf Sprik Third examiner : prof. dr. Ton van Leeuwen

(2)

Abstract

Platelets play an important role in hemostasis, the healing of wounds, and their transport behavior in blood is still not yet fully understood. In this thesis the transport of platelets is studied with a two-dimensional model of fully resolved blood flow. The flow of the blood plasma is simulated with the lattice Boltzmann method and the mechanical models for red blood cells (RBCs) and platelets are coupled to the plasma with the immersed boundary method (IBM). The speed of the margination of platelets, the viscosity of the suspension and the diffusivity of platelets are measured in channel flow simulations with different setups to get a better understanding of the transport properties of platelets. The distribution profiles of the modeled blood agree with experimental results. The speed of the platelet transport into the RBC-free layer is quantified and compared with the equilibration time of the blood plasma and the RBCs.

This thesis is extended from 12 EC to 18 EC. Roughly the division is 12 EC for chapters 1, 2 and 3 (understanding fluid mechanics and the model), and 6 EC for

(3)

Populaire samenvatting

Natuurkundigen zoeken naar verklaringen voor alles wat zich in de wereld afspeelt en ze proberen wat er gebeurt in modellen te vangen. In dit onderzoek is gekeken naar de beweging van bloedplaatjes binnen bloed. Ziektebeelden die gepaard gaan met ve-randerdende reologie (stroomeigenschappen) van bloed kunnen vervelende consequenties hebben, zoals een zuurstoftekorten, een beroerte of een infarct. Een goed begrip van de reologie van bloed helpt bij het zoeken naar oorzaken en behandelingen van dit soort ziektebeelden. Ook draagt de bestudering van bloed bij aan betere kennis over de eigen-schappen van gelijksoortige mengsels, zoals gel en verf.

Bloed is geen vloeistof maar een suspensie: in het vloeibare bloedplasma (ongeveer 55% van het volume) zitten onopgeloste rode bloedcellen (ongeveer 40-45% van het volume), witte bloedcellen en bloedplaatjes. Wanneer bloed door grote aderen stroomt gedraagt het zich net als water (als een zogenaamde Newtoniaanse vloeistof). Als aderen echter dunner dan 500 micrometer worden gaan de onopgeloste deeltjes een rol spelen en lijkt het gedrag van bloed meer op dat van ketchup: hoe meer druk je er op uitoefent, hoe makkelijker het gaat stromen. De precieze natuurkundige beschrijving van bloed wordt nog ingewikkelder gemaakt doordat rode bloedcellen vervormbaar zijn. Ze kunnen bijvoorbeeld gemakkelijk bewegen door bloedvaten die half zo dun als zijzelf in normale toestand zijn.

Bloedplaatjes spelen een belangrijke rol in de stolling van bloed. Al sinds 1931 is bekend dat bloedplaatjes zich voornamelijk aan de rand van bloedvaten ophouden, terwijl de rode bloedcellen zich in het midden van het vat bevinden. Met het oog op de functie van bloedplaatjes bij bloedstolling is het maar goed dat ze zich aan de rand, waar scheurtjes optreden, bevinden. Natuurkundig is het alleen nog niet duidelijk wat de precieze oorzaak is.

Om het gedrag van vloeistoffen en gassen te beschrijven worden de vergelijkingen van Navier-Stokes gebruikt. Deze vergelijkingen kunnen tot op heden alleen analytisch - met algebra - opgelost worden in simpele gevallen. Om uit te kunnen rekenen hoe vloeistoffen zich gedragen in meer ingewikkelde situaties kunnen de Navier-Stokes vergelijkingen numeriek benaderd worden. In dit onderzoek is de rooster-Boltzmann methode gebruikt als numerieke benadermethode om het gedrag van bloed na te bootsen. Met de rooster-Boltzmann methode wordt vloeistof voorgesteld als een verzameling deeltjes dat zich op een rooster bevindt. Het voordeel van deze methode boven andere methoden is dat deze relatief weinig rekentijd op een computer kost.

Het voordeel van een werkende computersimulatie is dat parameters makkelijker kunnen worden veranderd dan in een laboratorium. Ook kan je in een computersimulatie op elk tijdstip op elke plek in je simulatie bepaalde natuurkundige waarden (zoals snelheid of druk) meten. Door de breedte van de gesimuleerde bloedvaten, de snelheid van het bloed en het volumepercentage rode bloedcellen in de simulaties te veranderen is gezocht naar de oorzaak voor de ruimtelijke verdeling van bloedplaatjes en rode bloedcellen.

(4)

Introduction

Diseases involving changes in the flow properties of blood take a lot of casualties, by resulting in for example strokes or infarcts. Only in the Netherlands, 39048 people died of cardiovascular diseases in 2012 (1). Many of these diseases affect blood flow. When for example a blood vessel is (partly) blocked by sclerosis or a weakened vessel wall led to an aneurysm, the blood flow around the affected parts is different than the flow in a healthy vessel.

Blood consists of red blood cells, white blood cells and platelets suspended in blood plasma. In this work only the physical, and not the chemical properties of blood will be considered, to keep the parameter space limited. There is already a relatively good understanding of the flow of hard sphere suspensions (2), but since red blood cells are deformable soft sphere particles, describing the dynamics of blood is more difficult. For soft sphere suspensions there is still a lack of good understanding of how collective behaviour, deformability, and orientation in the flow impacts rheology and transport properties. By studying blood flow one therefore helps both physicians trying to heal cardiovascular diseases and physicists trying to understand suspension flow.

The ultimate goal of physics is to come up with an analytical solution for the flow profiles studied. As will be shown this can be done for simple situations, such as the flow of a homogeneous liquid (without any suspended particles) in straight tubes (section

1.5). By adding suspended cells the dynamics get more complex, because the cells not only interact with the fluid, but also with each other. When the cells studied are then soft (as is the matter with blood) instead of rigid, things get even more difficult; the deformability gives the cells more degrees of freedom. Attempts have been made to write down equations that describe the flow profile of blood, but they don’t seem to be in good agreement with experimental results (section 2.4).

In the twentieth century numerous numerical techniques have been developed to study the flow of complex fluids with computers. Because the fast increasing power of (su-per)computers many papers have been published in the last twenty years on modeling blood flow (3–6). This work uses Suspensions of Membraneous Objects (SuMO) (7), a computer model developed in the Computational Science research group at the Univer-sity of Amsterdam. SuMO is a two dimensional model of suspensions of deformable red blood cells (RBCs) and platelets, which allows the study of various quantities related to fluid dynamics. The focus in this work is on the dynamics of platelets. Platelets play an important role in hemostasis: the process which causes bleeding to stop.

This thesis starts with the necessary background information. Section 1 considers the basics of fluid mechanics. In section2the components of blood and blood flow properties are described. On what ideas SuMO is based is written down in section3. How all these parts come together in the study of the dynamics of platelets is then described in section

4. The reader already familiar with fluid mechanics, the rheology of blood and SuMO could jump to chapter 4to read about the outcomes of this thesis.

(6)

1 Fluid mechanics

Fluid mechanics concerns itself with the study of fluids (liquids, gases and plasmas) and the forces working on fluids. The dynamics of fluids are described by conservation equa-tions for mass, momentum and energy. This chapter gives the derivation of the Navier-Stokes equations for incompressible fluids in sections 1.1–1.4, as well as an overview of relevant quantities and concepts in fluid mechanics in sections1.5and 1.6.

This chapter is based on text books (8–10).

1.1 The continuum hypothesis

At a microscopic scale, a fluid consists of different molecules with their own set of non-uniformally distributed physical properties (like velocity and density). The phenomena studied in fluid dynamics, however, are macroscopic. Which means that a small volume element in the fluid is still so large that it contains a great number of molecules. By averaging over the fluctuating (Brownian) velocities of the individual molecules we can assign a fluid velocity v(x, t) to an element at point x at time t. Following the same line of arguments we can define averages for the other relevant quantities of the fluid. Because of the continuum hypothesis these quantities all vary smoothly with x on the macroscopic scale, so we can use continuum calculus to calculate properties of the fluid. This makes for easier calculation while still producing very accurate results.

1.2 Conservation of mass: the continuity equation

Consider a volume V bounded by a surface S, fixed in space. The mass inside is given by R_V ρ dV , where ρ(x) is the density of the fluid. The decrease of mass over time in this volume is − Z V ∂ρ ∂t dV. (1.1)

The mass flowing through a surface element dS is ρv · dS. By integrating over the whole surface S, an expression for the total mass flowing out of V is obtained: R

Sρv · dS. Applying Green’s formula to the latter gives

Z S ρv · dS = Z V ∇ · (ρv)dV. (1.2)

To make sure that mass is conserved, the decrease of mass in volume V should be equal to the total mass flux out of V, for any volume. So equations (1.1) and (1.2) should be equal

∂ρ

∂t + ∇ · (ρv) = 0. (1.3)

The integrals vanish because this has to hold for any volume V . Equation (1.3) is called the continuity equation.

(7)

1.3 Conservation of momentum: the Cauchy equation

To derive the conservation of momentum in fluids, Newton’s second law (F = ma) is applied to a volume V with surface S. The rate of change of the momentum (ma) is given by

Z

V ρdv

dt dV. (1.4)

The forces acting on the volume are divided in two categories:

1. Long ranged external body forces that act equally on all the material in a volume element. We only consider gravity, dFg = ρg · dV . (Other examples of body forces are electric and magnetic forces.)

2. Short ranged surface forces. These forces work on a surface layer of the element considered, as is shown in figure 1.1. In three dimensions, a volume is bound by three surface planes with three force components each. This gives a total of 9 components, which form the stress tensor σ:

σ = σij =   σ11 σ12 σ13 σ21 σ22 σ23 σ31 σ32 σ33  ≡   σxx σxy σxz σyx σyy σyz σzx σzy σzz  ≡   σx τxy τxz τyx σy τyz τzx τzy σz  , (1.5)

where σ11, σ22 and σ33 are normal stresses and the other six are shear stresses, noted with τij.

Figure 1.1: A volume V with surface stresses σij, unit vectors ei and total stresses

T(ei)_. _(Source: _{http://upload.wikimedia.org/wikipedia/commons/b/b3/Components_}

stress_tensor_cartesian.svg)

The total force on volume V is given by the body force plus the surface force F = Fg+ FS= Z V ρg · dV + Z S σ · dS (1.6)

(8)

By applying Green’s formula to the latter integral: R_Sσ ·dS =R_V ∇·σdV , the expression for the total force becomes

F = Z

V

(ρg + ∇ · σ)dV. (1.7)

We can now write Newton’s second law by combining (1.4) and (1.7) ρdv

dt = ρg + ∇ · σ. (1.8)

The integrals vanish because this has to hold for any volume V . Equation (1.8) is known as the Cauchy equation.

Newton’s first law states that all moments must balance. Consider a cube with edges x1, x2 and x3 and look at the moment about one of the axes. The force exerted by stress σ12 is F = σ12· x1· x2. The moment arm is the distance to the middle of the cube, so the moment M = σ12· x1· x2·x₂3. The total of forces on the four planes about this axis is σ12· x1· x2· x3 2 + σ12· x1· x2· x3 2 = σ21· x1· x3· x2 2 + σ21· x1· x3· x2 2 , σ12· x1· x2· x3 = σ21· x1· x3· x2, σ12= σ21.

Because the same holds if the moment about the other two axes is considered, it can be seen that the stress sensor is symmetric, and thus the following holds

σij = σji. (1.9)

Thus there are six different surface stresses. The stresses arise from a combination of pressure p and viscous friction. Pressure is always acting normal and inward to the sur-face, and is therefore only apparent in σxx, σyyand σzz. Because it is directed inward the surface, it gets a negative sign. The viscous stresses are generated by velocity gradients according to σij = µ∂v_∂ji, where µ is the dynamic viscosity. Putting this together, and writing v = (u, v, w) gives

σxx = −p + λ∇ · v + 2µ ∂u ∂x, σxy = σyx= µ( ∂u ∂y + ∂v ∂x), (1.10) σyy= −p + λ∇ · v + 2µ ∂v ∂y, σyz= σzy= µ( ∂v ∂z + ∂w ∂y), (1.11) σzz = −p + λ∇ · v + 2µ ∂w ∂z, σzx= σxz = µ( ∂u ∂z + ∂w ∂x). (1.12)

The relationship between stress and velocity gradients is assumed to be linear here (which is true for Newtonian fluids) and it is assumed that the fluid is isotropic and thus has no preferred direction.

Important to note is that fluids with a higher viscosity have more difficulty flowing. Water, for example, has a low viscosity, and honey has a higher viscosity.

(9)

1.4 Incompressible Navier-Stokes equations

Incompressible fluids are fluids with constant density ρ, so the equation of continuity (equation1.3) for incompressible fluids gets

∇ · v = 0. (1.13)

Equations (1.10), (1.11) and (1.12) thus reduce to σij = −pδij + µ( ∂vi ∂xj +∂vj ∂xi ). (1.14)

By combining (1.8) and (1.14) and expanding the total derivative of v as a partial derivative (dv/dt = ∂v/∂t + v · ∇v) the famous Navier-Stokes equation for fluid flow is obtained ∂v ∂t + v · ∇v = − 1 ρ∇p + ν∇ 2_v. _(1.15)

The kinematic viscosity ν is the dynamic viscosity µ divided by the density ρ. The gravity term and the pressure term cancel out. Because we can’t work with infinite viscosities and because experimental results show this, there holds a no-slip boundary condition which prevents the fluid from slipping relative to the wall. Thus velocity v = 0 at the walls.

As we are studying blood flow in 2D, it is useful to write down the Navier-Stokes equa-tions for incompressible fluids in 2D. Writing v = (u, w) equaequa-tions (1.13) and (1.15) reduce to ∂u ∂x+ ∂v ∂y = 0, (1.16) u∂u ∂x+ v ∂u ∂y = − 1 ρ ∂p ∂x+ ν( ∂2u ∂x2 + ∂2u ∂y2), (1.17) u∂v ∂x+ v ∂v ∂y = − 1 ρ ∂p ∂y + ν( ∂2v ∂x2 + ∂2v ∂y2). (1.18)

1.5 Poiseuille and Couette flow

When studying fluid dynamics different types of flow can be considered. This work uses Poiseuille and Couette flow. Poiseuille flow is the flow of a fluid between two parallel rigid1 walls, where a pressure gradient along the direction drives the flow of the viscous fluid. In the case of Couette flow the movement of one of the two walls - with velocity u - is driving the flow of the fluid.

Poiseuille flow considers pressure driven constant flow of a viscous fluid between two infinite rigid parallel plates. In the computer model used in this work the pressure is due

(10)

to gravity, Fg = ρg. Solving the Navier-Stokes equation with these conditions results in flow velocity profile2

v(y) = −ρg

2νy(y − h). (1.19)

Couette flow does also consider the flow of a viscous flow between two rigid parallel plates. Only now one of the plates is moving with velocity u relative to the other plate, resulting in a drag force on the fluid. Solving the Navier-Stokes equation with these conditions results in flow velocity profile3

v(y) = uy

h. (1.20)

Figure 1.2: Schematic representation of Poiseuille flow (left) and Couette flow (right).

In fluids different layers of fluid move over other layers of fluid. Layers can have different velocities (as can be seen in for example the flow profiles in figure 1.2). Shear rate is a measure for the velocity gradient (∂v/∂y) along the diameter of flow, it is used in this work to distinguish different regions in flow profiles.

As can be seen in equation (1.19) Poiseuille flow has fluid velocity vy ∝ y2 and a shear rate ˙γ ∝ y. From equation (1.20) follows that Couette flow has vy ∝ y and ˙γ = c (c is a constant).

Combinations of Poiseuille and Couette flow, in many - often complex - geometries are abundant in nature, but for the sake of simplicity this work only considers simple Poiseuille flow in channels.

1.6 Relevant quantities

The Reynolds number is a dimensionless number, which is the ratio of the inertia force to the viscous force in the flow

Re = vd

ν , (1.21)

where v is the velocity of the flow, d the vessel diameter and ν the kinematic viscos-ity.

2

SeeA.1for the calculation.

(11)

2 Hemorheology

This chapter describes the flow properties of blood and its elements. Section 2.1 gives an overview of the properties of blood plasma, red blood cells (RBCs) and platelets. The flow properties of blood as well as relevant concepts as for example the F˚ ahræus-Lindqvist effect and the margination of platelets are described in section2.2. The phys-ical concept of diffusion is laid out in section 2.3. Finally earlier attempts to describe blood flow mathematically are named in section2.4.

2.1 Elements of blood

Blood consists of red blood cells (RBCs), platelets and white blood cells suspended in plasma. Roughly the volume is divided: 55% for plasma, 40-45% (their is a slight differ-ence between RBC concentrations for men and women) for RBCs and a few percentages for platelets and white blood cells. In this work white blood cells are neglected, because of their small volume percentage (circa 0.7%; i.e. one white blood cell for every thousand RBCs).

Plasma is a Newtonian fluid, which means that its viscosity scales linearly with the strain on the fluid. It is mostly made of water (> 90%).

RBCs - also called erythrocytes - are deformable biconcave discs. Under normal con-ditions they have a diameter of 8µm and a thickness of 2.5µm. RBCs are made of a highly flexible membrane and enclosed cytoplasm (11). The main task of RBCs is to deliver oxygen (via oxygen carrier hemoglobin) to all parts of the body. Because their deformability RBCs can squeeze through small vessels. The smallest vessels in the hu-man body are called capillaries. They are smaller than half the diameter of a RBC. In the capillaries the cells and the blood exchange substances. The hematocrit - or RBC volume ratio - is 40.0 ± 2.4% for women and 45.8 ± 2.7% for men (12). This adds up to around five million RBCs per mm3_.

Figure 2.1: From left to right: a red blood cells, a platelet and a T-lymphocyte (white blood cell), as photographed by a scanning electron microscope. (Source: http://upload.wikimedia. org/wikipedia/commons/2/24/Red_White_Blood_cells.jpg)

(12)

Platelets - also called thrombocytes - are irregularly shaped discs (as can be seen in figure2.1) with a diameter of 2 − 3µm (11). They play an important role in hemostasis; the process which causes bleeding to stop. The platelet concentration in blood is 2.5 − 5.0 × 105 platelets per mm3 (13).

2.2 Blood flow properties

There is a growing interest in the rheology of soft suspensions. More specifically there is broad interest in the rheology of blood - so called hemorheology. A range of diseases is related to modified rheological properties of blood, causing for example undersupply of oxygen, stroke or infarction. A proper understanding of hemorheology helps to treat those diseases. The previous chapter already dealt with fluid dynamics in general, and the previous section with the constituents of blood. How blood flows (or how we think blood flows) is described in this section.

The largest blood vessels in the human body are a few millimeters to even centimeters wide4_{. The smallest vessels (capillaries) and also RBCs have a diameter of a few} mi-crometers. So there is a range of three orders of magnitudes in vessel diameters. This leads to an interesting palette of physics. Reynolds numbers vary from Re ≈ 0.001 in small vessels to Re ≈ 1000 in large vessels.

Blood is a shear thinning suspension for shear rates smaller than ˙γ = 100 s−1. Shear thinning means the viscosity of the fluid decreases with an increase of stress5. For shear rates above 100 s−1 blood behaves as a Newtonian fluid. The shear thinning of blood is related to three mechanisms (14):

1. At low shear rates: aggregation of RBCs takes place; the RBCs clutch together in so called rouleaux, causing higher apparent viscosities. It is not clear whether this is a hydrodynamic or a biomechanical effect.

2. For intermediate and high shear rates: deformation, elongation and tank-treading of RBCs is observed.

3. For high shear rates: rouleaux break down and RBCs align with the flow, which leads to a Newtonian behavior.

The viscosity of blood further seems to depend on the red blood cell concentration (hematocrit), the temperature, the viscosity of the blood plasma and the health of RBCs (15).

The decrease of the diameter of a blood vessel leads to the decrease of the hematocrit. This effect is called the F˚ahræus effect and starts when the channel diameter gets smaller

4_{The aorta of an adult has a typical size of 2-3 cm.} 5

This is why you’ll always end up with too many tooth-paste on your tooth-brush and too much ketchup in your fries.

(13)

Figure 2.2: Distribution profile of red blood cells (—) and platelets (- -) as observed in vitro by (17). A cell-free layer without RBCs near the wall of the channel and an excess of platelets near the walls of the channel is visible.

than 500µm. It was first observed in vitro6 in 1929. A closely related phenomenon now known as the F˚ahræus-Lindqvist effect was first described in 1931 (16). The effect says that the viscosity of blood changes with the diameter of the channel it travels through. This lowering of viscosity is caused by the migration of RBCs towards the center of the channel, leaving a cell-free layer (CFL) near the walls. This small CFL has a lower viscosity than the rest of the blood and serves as a kind of lubrication layer. Because the CFL is only a few µm thin, this effect is only significant in vessels with a small diameter. The motion of RBCs causes platelets to migrate near the walls of the vessels (which leaves them nearby the place where they are needed for hemostasis).

One theorem related to the migration of platelets towards the walls of the vessel is the Segr´e-Silberberg effect. In the early 1960’s Segr´e and Silberberg noted that particles in a flow are subjected to viscous drag forces an inertial lift forces, leaving them with a net force driving them towards the center of flow channel (18). The viscous drag force drives particles along the streamlines of the flow, while the parabolic nature of Poiseuille flow (as can be seen in figure 1.2) gives the particles a shear-induced inertial lift force towards the walls of the channel. The flow around the particles induces a pressure increase between the particles and the walls of the channel. This pressure prevents the particles to move even closer to the wall, causing a cell free layer. How the distribution of platelets in blood without any red blood cells looks like is showed in figure2.3. Originally the experiments with this effect were done with suspensions of solid particles, but it also seems to explain the profile in figure2.3. It might be that the deformability of the RBCs makes it easier for RBCs to move in the center of the vessels pushing the platelets out to the walls, thus pushing the two peaks visible in figure2.3 to the walls of the vessel, resulting in a profile similar to the one in figure 2.2.

6

In biological experiments the important distinction between in vitro and in vivo is made, the first meaning experiments in the lab (outside an organism) and the latter meaning experiments inside an organism.

(14)

Figure 2.3: Distribution of ˆCplt over the channel width d in flow without RBCs. The blue line

is the simulation result and the red dots represent in vitro data from Aarts et al. (17). (Source: Mountrakis et al., 2013 (19).)

2.3 Diffusion

Transport phenomena are all phenomena in nature concerning the transport of mass, energy and momentum. One example of a transport process - fluid flow - is already discussed in chapter1, and another important one is diffusion. Diffusion considers the net movement of randomly moving particles. For this work the diffusion of particles is interesting; another form of diffusion is heat diffusion. A diffusion approach is also used to study for example population dynamics and financial systems.

Particle diffusion is the movement of a concentration c(x, t) from one region to another. Because forming a certain distribution profile in blood involves the concentration move-ment of RBCs and platelets it is interesting to compare the diffusion profile over time with the distribution profile over time. In many cases diffusion is described by Fick’s law and the diffusion coefficient D

∂

∂tc(x, t) = D∇

2_{c(x, t).} _(2.1)

When D is not constant the right hand site becomes ∇(D∇c(x, t)). In most cases diffusion is caused by a difference in concentration gradients. But particle-particle in-teractions as a result of hydrodynamics are also driving diffusion; particles collide with other particles along their path as they are flowing in a suspension. This collision-driven diffusion is called shear-induced diffusion. Shear-induced diffusion is a deterministic pro-cess, but it can more easily be described as a diffusion propro-cess, due to the complexity of the hydrodynamics. Shear-induced self-diffusion (as opposed to shear-induced gradient diffusion) considers the motion of individual particles in suspensions. The shear-induced diffusion caused by suspended particles scales as

(15)

with R the radius of the suspended particles (RBCs), ˙γ the shear rate and φRBC the hematocrit (20). Thus the total diffusion in blood flow is

Deff = Dth+ D( ˙γ, φRBC), (2.3)

where Dth denotes the thermal diffusion.

Transport due to the bulk motion of a fluid is called advection, it is described by the advection-diffusion equation. For incompressible flow this is

∂

∂tc(x, t) = ∇ · (D∇c) − v · ∇c(x, t), (2.4)

where D is the diffusion tensor and v is the flow velocity.

If we consider the RBC concentration as the field RBCs and platelets are diffusing through, the diffusion tensor DRBCfor red blood cells depends on the hematocrit φRBC and the local shear rate ˙γ. An expression for DRBC can be obtained empirically from micro-model simulations by measuring the root mean square displacement of red blood cells. The diffusion tensor Dplt for platelets also depends on the hematocrit φRBC and the local shear rate ˙γ. Adding those dependencies to (2.4) gives

∂

∂tφRBC= ∇ · (DRBC(φRBC, v)∇φRBC) − v · ∇φRBC, (2.5) ∂

∂tφplt= ∇ · (Dplt(φRBC, v)∇φplt) − v · ∇φplt. (2.6) In a simulation or experiment the diffusion can be measured via the root mean square displacement of the platelets

D = hx 2_i

2t . (2.7)

It is important to only start measuring after some time (corresponding to observations), because it takes some time for all the particles to have had at least one collision.

2.4 Mathematical model to describe blood flow

The micro-model used in this work resolves the individual behavior of red blood cells and platelets suspended in plasma. The final goal of this research is to determine the driving forces of the margination of platelets, and possibly help in finding a mathematical model that describes the margination, so it would be possible to calculate the distribution of the platelets at time tx if we know the relevant quantities at time t0.

As is described in2.2the viscosity of blood depends - among other factors - on the shear rate ˙γ and on the local hematocrit φRBC (11), the contribution of the platelets to the

(16)

blood viscosity can be neglected. The shear-thinning behavior of blood seems to be well described by the Carreau-Yasuda model (5)

η( ˙γ) = η∞+ (η0− η∞)(1 + (λ ˙γ)a)

n−1

a , (2.8)

therefore this is one of the most used models to describe blood. The Carreau-Yasuda model describes the flow of non-Newtonian fluids with asymptotic viscosities at zero (η0) en infinite (η∞) shear rates and without yield stress. The latter is a critical stress; after reaching yield stress a fluid behaves as a solid. Fluids obeying the Carreau-Yasuda model behave as a Newtonian fluid at low shear rates and as a power-law fluid at high shear rates. λ is a constant, where 1/λ is the critical shear rate at which viscosity begins to decrease. a represents the width of the transition region between η0 and the power-law region.

The Carreau-Yasuda fails however to account for the dependence of viscosity on the local hematocrit, a better model is thus necessary for a more accurate description of blood.

(17)

3 Model

This work uses Suspensions of Membraneous Objects (SuMO) as a tool to simulate blood flow with fully resolved red blood cells (RBCs) and platelets. The usability of SuMO in studying aneurysmal vessels is already shown (19) and the simulated data is validated with experimental data (7). SuMO uses the Lattice Boltzmann Method (LBM) to solve the Navier-Stokes equations numerically. The immersed boundary method (IBM) is used to couple membranes (e.g. RBCs and platelets) to the fluid. This chapter introduces both the LBM (in section 3.1) and the IBM (in section 3.2) in their historical and mathematical contexts. The possibilities for the setup of simulations with SuMO are described in section 3.3 and finally the validation of the model is discussed in section

3.4.

3.1 Lattice Boltzmann Method

1 6 5

4

3 2

Figure 3.1: The Frisch, Hasslacher, Pomeau hexagonal lattice, particles on the site in the middle can have velocities in one of six directions.

In the late 1950s Stanislaw Ulam and John von Neumann were the first who used dis-crete systems called cellular automata (CA) to model liquids7_{. In 1986 Uriel Frisch,} Brosl Hasslacher and Yves Pomeau made a leap forward by producing a ‘class of deter-ministic lattice gases with discrete Boolean elements [that] simulates the Navier-Stokes equation’ (21). They considered a hexagonal lattice with up to six particles on each site. Particles have discrete velocities in directions i = 1, 2, .., 6 as shown in figure 3.1. Fur-thermore a Fermi exclusion principle holds, which allows only one particle per direction per lattice site. After every time step particles hop to another site on the grid according to their discrete velocity. This model is known as a type of Lattice Gas Cellular Automa-ton (LGCA). This LGCA Frisch, Hasslacher and Pomeau proposed conserves particle number, total momentum, and (the greatest achievement of this model) rotational in-variance. In the abstract of the aforementioned article the authors promise ‘simple, massively parallel computing machines’ by using their method, but unfortunately there were some problems with this approach: because only discrete numbers of particles are allowed on every site there are big (integer-sized) jumps in particle numbers, which does

7_{The early succes of this method lead computer pioneer Konrad Zuse to the publication of Rechnender}

(18)

not give a smooth profile. By looking not only at single particles, but at particle dis-tributions, while still considering discrete lattices the problem was solved. This way of modeling is called the Lattice Boltzmann Method. (A clear and thorough overview of the development from LGCA to LBM is given in a textbook written by Sauro Succi (22).)

SuMO uses a D2Q9 lattice (see figure3.2), where in DxQy the x stands for number of dimensions and the y for number of velocities. Particles have lattice velocities

c0 = (0, 0), c1= (1, 0), c2 = (0, 1), c3 = (−1, 0), c4= (0, −1), c5= (1, 1), c6 = (−1, 1), c7 = (−1, −1), c8 = (1, −1).

A Fermi exclusion principle holds: there can only be one particle with a velocity ci at each lattice site (so there can be up to 9 particles on each lattice site). During each time step ∆t particles travel from site to site. If they meet each other they collide according to certain rules in zero time. Mass, momentum and energy are conserved (the basic physical conservation laws). It is possible to calculate the density and the average momentum at every lattice site for every time step.

c2 c4 c8 c1 c5 c7 c3 c6

Figure 3.2: A cartesian grid with grid spacing ∆x = 1 and nine discrete velocities. This is a D2Q9 lattice.

LBM considers particle distributions instead of the movement of individual particles, therefore a probability density function f (x, v, t) is introduced which gives the probabil-ity of finding a particle at position x with velocprobabil-ity v at time t. Because velocities come from a discrete set, we can write fi(x, t) as the probability density function for velocity ci. To look at the evolution of the system we set up an evolution equation

f_i0(x, t) = fi(x, t) + Ωi(f ), (3.1) where Ωi(f ) is the collision operator; a continuous version of the collision rules for the particles on the lattice. The collision operator depends on the simulated phe-nomenon.

The Lattice Boltzmann Equation describes the dynamics of the probability density and is

(19)

From the probability densities we can calculate the fluid density ρ(x, t) and the fluid momentum ρ(x, t) · u(x, t). The algorithm can be described as follows:

1. Calculate densities and velocities.

2. Stream step: particles travel to the next node, according to their velocities. 3. Collision step: collisions of particles are handled.

4. Boundary conditions: interactions with walls or suspended obstacles.

In the scheme of scales the LBM is on meso-scale, in between the micro-scale of parti-cles and the macro-scale of continuum fields. A more extended overview of LBM and a derivation that shows that LBM resembles the Navier-Stokes equations for incompress-ible fluids is given in (23).

3.2 Immersed Boundary Method

To allow fluid-cell interaction the immersed boundary method is used. With IBM cells (or called membranes or structures) are immersed in the flow. Fluids are described by a discrete lattice system (an Eulerian grid) with LBM where each lattice node stays on a fixed position. The cells however are represented by Lagrangian surface points (LSPs) which can move continuously between the lattice nodes. An example is showed in figure

3.3. In IMB there thus are two coordinate systems, and the trick of IBM is the

com-Figure 3.3: A membrane (in red) immersed in a discrete rectangular lattice (in white). (Source: (24))

munication between the two systems. The communication is based on force interaction between the cells and the fluid. The choice for the number of LSPs is important and depends on the lattice spacing ∆x; a few examples are shows in figure 3.4. A precise mathematical overview of IBM is given in (25).

(20)

Figure 3.4: The choice of number of LSPs is important: f1 has too little LSPs, f2 has too many

LSPs and the number of LSPs of f3 is just right. (Source: (24))

3.3 Suspensions of Membraneous Objects (SuMO)

SuMO is a simulation code that is developed at the Computational Science research group at the University of Amsterdam. It uses the immersed boundary method on top of a lattice Boltzmann method with a Bathnagar-Gross-Krook operator (26)

Ωi = 1

τ(ni− n EQ

i ). (3.3)

In this operator function τ is a relaxation constant related to the viscosity of the fluid.

When setting up a computer simulation there is a trade-off between the computations done and to what extent the model resembles nature. More precision means more com-puter time or a more difficult algorithm. The parameters that SuMO uses in simulations are based on a combination of experimental data and trial-and-error results of simula-tions. The used parameters are listed in table 1. Every computer iteration represents ∆t = 9.8 × 10−8 s in real life. The distances between the nodes on the lattice are ∆x = 1 µm. The number of Lagrangian surface points (LSPs) representing a red blood cell, denoted by NRBC, is 26. Platelets have Nplt= 8. The numbers of LSPs is based on the ∆x and the IBM theory; IBM works best if a membrane has a LSP near all lattice nodes of its surface. The kinetic viscosity of the blood plasma is ν = 1.7 × 10−6 m2 / s. The spring constant Cspr, the bending constant Ctrs, the bending viscosity Dtrs, the area constant Carea and the cell-cell constant Crep govern the behaviour of the membranes representing the RBCs and platelets. How RBCs and platelets look like in SuMO is shown in figure3.5.

The maximum input velocity of the fluid can be tweaked by changing the body force gx, according to

umax=

gx· D2

(21)

Figure 3.5: Representation of a RBC (left) and a platelet (right) in SuMO.

Parameter Value

Time step, ∆t 9.8 × 10−8 s

Spatial resolution, ∆x 1 µm

Number of LSPs per RBC, NRBC 26 Number of LSPs per platelet, Nplt 8

Kinetic viscosity, ν 1.7 × 10−6 m2/s

Spring constant, Cspr 0.2 N/m

Bending constant, Ctrs 9.8 × 10−8 N/rad Bending viscosity, Dtrs 0.0 N·s/rad

Brea constant, Carea 9.8 × 104 N/m2

Cell-cell constant, Crep 1.97 × 10−19 N·m2

Table 1: Model parameters and constants.

The Reynolds number Re (1.21) is dimensionless so quantities can either be given in lattice Boltzmann units or in SI units. Due to the complex nature of suspensions the resulting Reynolds number in a simulation can be different than the input Reynolds number. The value of the real Reynolds number is approached by using umax (the maximum velocity in the channel) after equilibration of the simulation and the imposed kinematic viscosity ν.

3.4 Validation of SuMO

For a simulation model it is important that the results correspond to experimental data. It is validated that SuMO recovers the F˚ahræus-Lindqvist effect and the shear thinning of blood. The size of the RBC-free layer and also the tank-treading frequency of single cells matches experimental data (7).

(22)

4 Results

The results of this thesis are presented in this chapter. The used setup of simulations is described in section 4.1. The equilibration of the simulations is considered first, in section 4.2. Section 4.3 then presents the study of the margination of platelets. In section 4.4the ways of measuring diffusion in the micro-model are discussed.

4.1 Simulation setup and post processing

Simulations with 32 different setups were run on the Lisa Compute Cluster, as described in table 2. For each setup 8 simulations with different initial conditions were run in parallel to increase statistics. This makes for a total of 32 × 8 = 256 simulations. The used parameters were chosen to compare the results with recent studies (3,4,7).

Vessel diameter × length [µm2] Reynolds number Hematocrit [%]

20 × 40 1 15, 30, 40, 50

40 × 80 5, 10 15, 30, 40, 50

100 × 200 5, 25, 50, 100 15, 30, 40, 50

150 × 300 10 15, 30, 40, 50

Table 2: Setup of SuMO simulations. Noted is the imposed Reynolds number, the real Reynolds number is different, due to the presence of particles.

Figure 4.1shows a snapshot of one of the simulations. This snapshot is taken ≈ 0.42 s after initialization: platelets are now mostly living near the walls of the vessel and a cell free layer is apparent.

Figure 4.1: A snapshot after ≈ 0.42 s in a 100×200 µm channel with Re ≈ 10 and H = 40%. The biconcave shaped membranes represent red blood cells and the round shaped membranes represent platelets.

(23)

With Python scripts the velocities, shear rates ˙γ, hematocrit H, platelet concentration ˆ

Cpltand diffusivities D were measured. The processing of the raw SuMO data for every setup was done in two steps:

1. For every case with different initial conditions the 2D data was collapsed (averaged) to a 1D array: i.e. 1 µm × the length of the vessel → 1 × 1 µm. Resulting in arrays of data with a length equal to the vessel width.

2. The data from the 8 cases with different initial conditions were averaged, resulting in just one array for every quantity per setup. Not all 8 cases reached the same tmax so the average is made up of 3-8 different cases. The standard deviations are calculated and added in the plots8.

In SuMO, cells are represented by Lagrangian surface points (LSPs). The raw data contains the velocities and positions of those LSPs for every time step. To calculate the concentration of particles for a certain region the radius and volume ratio are required. RBCs and platelets are irregularly shaped; their volume ratios are defined as their cell volume divided by the volume of a sphere with the same radius. The radii and volume ratios are based on experimental work and are noted in table 3(11).

RBCs Platelets

Radius [µm] 2.66 0.80

Volume ratio 0.45 0.9

Table 3: Quantities for calculating concentrations of RBCs and platelets.

The volume of a particle is

Vparticle = Cπr2, (4.1)

where C is the volume ratio. To calculate the volume per LSP we divide the volume of the cells by the number of LSPs for that particle. The number of LSPs representing a cell, the resulting RBC and platelet volumes and volumes per RBC LSP and platelet LSP (needed for calculations) can be found in table 4.

RBCs Platelets

LSPs per cell 26 8

Volume [µm2] 10.00 1.810

Volume per LSP [µm2] 0.3847 0.2262

Table 4: Volumes of RBCs and platelets and RBC LSPs and platelet LSPs in SuMO.

8_{The standard deviations, and thus the error bars, will be bigger when the average is taken over}

a smaller number of different cases; the transitions between averages from ni to nj files is not noted

(24)

The hematocrit H or φRBC (RBC concentration) for a region is the volume occupied by RBCs divided by the total volume. For practical reasons it can be calculated by counting the number of RBC LSPs in a region and multiplying this number by the volume per LSP. The same holds for the platelet concentration φplt.

When looking at the distribution of platelets it is more useful to plot the normalized platelet concentration ˆCplt, which is calculated with

ˆ Cplt= φplt ¯ φplt , (4.2)

where φplt is the local platelet concentration and ¯φplt is the average platelet concentra-tion. A normalized platelet concentration ˆCplt of 1 is thus equal to the average platelet concentration. To quantify the speed of the margination of the platelets it is useful to look at the particle near wall excess (PNWE), the percentage of platelets inhabiting the two peaks (one on the left and one on the right) near the wall. With PNWE = 100% all the platelets ‘live’ in the peaks near the wall.

Because the Reynolds number (1.21) is dimensionless and thus the same in all systems, we have for SI and lattice Boltzmann units

m ·m_s m2 s = ∆x · ∆x ∆t ∆x2 ∆t . (4.3)

We can therefore use

∆t = νLB νSI

∆x2, (4.4)

to converse between the two unit systems.

For all data the standard deviation is used as uncertainty. There are two cases of propagation of uncertainty used, they are calculated in the standard way

σ(x+y)= q (σx)2+ (σy)2, (4.5) σ_(x/y) = x y s (σx)2 x2 + (σy)2 y2 . (4.6) 4.2 Equilibration

Simulations in SuMO start with a homogeneous distribution of cells. It therefore takes some time for the suspension to reach an equilibrium.

As a start it is useful to check whether the fluid in the simulation data conserves mo-mentum p = mv. Since SuMO results live on a discrete 2D lattice the total momo-mentum is the sum over the momenta of the lattice nodes

ptot = X

ij

(25)

with ρ the density and u the velocity.

In figure 4.2 the development of the total momentum for respectively 3 and 4 runs of two different setups is presented. As seen in 1.3 the momentum should be conserved in fluid dynamics. Though it seems like the total momentum equilibrates to a certain value, it still fluctuates, which is not in resemblance with nature. Red blood cells are deformable, thus they store energy. This may cause the fluctuations.

Figure 4.2: The total momentum of the fluid in simulations with a vessel of 100 × 200µm, and Re ≈ 1.5 and H = 40% and Re ≈ and H = 40%. Different colours represent different runs with different initial conditions.

In (27) Nott and Brady propose an equation for the length Leqa suspension has to travel before reaching equilibrium

Leq d ∼ 1 12 ˆD(φ) d2 R, (4.8)

where d is the width of the channel, ˆD(φ) the scaling of the diffusion and R the particle radius. Trying to apply this to blood Leq can be written as teqv and following (28)

ˆ D ∼ ˙γR2φRBC teq ∼ 1 12v ˙γφRBC d3 R3. (4.9)

The equilibration times for the fluid velocity, RBC-free layer and the margination of the platelets (more extensively discussed in section 4.3) are shown in the plots on the right side in appendixB. The quantities that change are the vessel diameter d, the fluid velocity v (by changing the imposed pressure gradient) and the hematocrit φRBC(or H). It is difficult to give a good definition of when a system is equilibrated, because most systems keep fluctuating to some extent (as can also be seen in the plots in appendix

B. By choosing the equilibration time as the time until the system starts to only shows small or regular fluctuations the data in table5 is extracted from the data in appendix

B. One can see that in most cases the fluid (represented in the table as equilibration time of umax) equilibrates first, followed by the RBCs. The platelets are the last constituents to form an equilibrium distribution. As in (4.9), it indeed takes longer for the simulated

(26)

equilibrium position. Moreover, faster flowing blood and blood with higher hematocrit H both take less time to equilibrate. For faster flowing blood the cells are in their equilibrium position earlier, and for blood with higher hematocrit holds that platelets are pushed out of the centre of the vessel faster because there are more RBCs to collide with.

Setup Equilibration times Speed

d × l [µm2] H [%] Re PNWE [s] umax [s] CFL [s] umax [m/s]

20 × 40 15 1 0.2 0 0.5 0.65 20 × 40 30 0.5 0.05 0 0.25 0.5 20 × 40 40 0.5 0.1 0 0.25 0.4 20 × 40 50 0.5 0.1 0.2 0.3 0.3 40 × 80 15 4 0.25 0 0.1 0.17 40 × 80 30 3 0.2 0 0.1 0.13 40 × 80 40 2.5 - 0.1 0.2 0.1 40 × 80 50 2.5 - 0.5 0.75 0.9 100 × 200 15 4 0.6 0.3 0.6 0.65 100 × 200 30 2.5 0.5 0.3 0.4 0.4 100 × 200 40 1.5 - 0.3 - 0.25 100 × 200 50 0.5 - - - 0.12 100 × 200 15 20 0.55 0.2 0.2 0.34 100 × 200 30 14.5 0.25 0.05 0.15 0.25 100 × 200 40 10.5 - 0.05 0.05 0.18 100 × 200 50 7.5 0.13 0.14 0.14 0.12 100 × 200 15 41.5 0.14 0.04 0.06 0.7 100 × 200 30 31.5 0.2 0.05 0.05 0.52 100 × 200 40 24.5 0.12 0.06 0.05 0.42 100 × 200 50 19 0.1 0.12 0.04 0.32 100 × 200 15 85 0.05 0.02 0.04 1.4 100 × 200 30 67 0.07 0.04 0.05 1.1 100 × 200 40 53 0.05 0.04 0.02 0.9 100 × 200 50 43.5 0.06 0.06 0.04 0.8 150 × 300 15 7.5 0.6 0.2 0.45 0.09 150 × 300 30 4.5 0.25 0.1 0.1 0.05 150 × 300 40 2.5 0.02 0.04 - 0.03 150 × 300 50 1 - - - 0.01

Table 5: Equilibration times for the platelets (PNWE), blood plasma (umax) and RBCs (CFL)

and speed of the fluid after equilibration. Also shown are thexs width of the channel d, the length of the channel l, the hematocrit H and Reynolds number Re. ‘-’ denotes a lack of proper data.

If we calculate the averages of all the equilibration times in table 5 we get an average equilibration time for the fluid of 0.11 s, for the RBCs of 0.20 s and for the platelets of

(27)

0.21 s. This is not very useful since the setups of the systems averaged over are different, but it shows that the equilibration order is indeed: (i) blood plasma, (ii) RBCs and (iii) platelets. Note that the platelets aren’t really equilibrated at the time measures given here, since very slowly all platelets will move to the two walls until the PNWE is equal to 100 % (more on this in the next section). The equilibration times thus seems to be dependent on the hematocrit, the vessel width and the velocity of the fluid, as (4.9) proposes, but more measurements are needed to see if the ratio of the dependency (linear, quadratic, etc.) as proposed by Nott and Brady for hard sphere suspensions also holds for blood.

4.3 Speed of margination

Every simulation starts with a homogeneous distribution of both platelets and red blood cells, but over time develops the typical cell distribution of blood. Soon after the start of the simulation the margination of platelets starts to take place. To illustrate this margination of platelets this process is plotted in figure4.3for one of the setups. Start-ing with an average platelet concentration over the whole width, this vessel ends up with platelet concentration excesses up to 20 times the average concentration near the walls.

Figure 4.3: The margination of platelets in a channel of 100×200 µm, with Re ≈ 31.5 and H = 30%. The color bar denotes the normalized platelet concentration as in (4.2).

The hematocrit, platelet concentration, shear rate and fluid velocity over the channel width for two moments in the same simulation as in figure4.3are shown in figure4.4. At t = 0.001 s both the red blood cells and platelets are distributed homogeneously, while at t = 0.22 s the platelets have marginated toward the cell-free layer without RBCs. The RBC-peak in the centre of the vessel is not observed in experiments, but it is visible in

(28)

0.22 s can also be seen by looking at the smaller standard deviations compared to the distribution profiles at t = 0.001 s.

Figure 4.4: Distribution of H, ˆCplt, ˙γ and Uxover the channel width d in a 100×200 µm, with

Re ≈ 30 and H = 30% at t = 0.001 s (left) and t = 0.22 s (right). The color filled lines denote the standard deviation.

The size of the cell-free layer (CFL) for different hematocrits is showed in figure4.6. The CFL is larger in simulations with lower hematocrit. The distribution of platelets in flow without RBCs is shown in figure 2.3. So it looks like platelets want to occupy around 60 % of the width of the channel. In blood flow with high hematocrit the platelets get pushed to the walls, but in lower hematocrit their distribution profile is more like their distribution without RBCs. In figure4.5the margination of platelets is shown for blood with a lower hematocrit, it tends more towards the profile in figure 2.3, because there are less RBCs the RBCs need a smaller volume in the centre of the vessel.

Figure 4.5: Left: The margination of platelets in a channel of 100×200 µm, with Re ≈ 67 and H = 15%. The color bar denotes the normalized platelet concentration as in (4.2). Right: the distribution of H, ˆCplt, ˙γ and Uxover the channel width in this system after ≈ 0.6 seconds.

To measure the speed of the margination and compare the margination rates for different setups two quantities are of special interest. Firstly the width of the CFL and secondly the aforementioned PNWE. In figure4.6the width of the CFL is presented as a function

(29)

Figure 4.6: The cell free layer width in µm over H and over Re for setups with vessel width = 20, 40, 100 and 150 µm. Same-coloured cubes had the same imposed Reynolds numbers.

of hematocrit and as a function of the Reynolds number. A lower hematocrit leads to a bigger cell free layer. The two peaks belonging to the PNWE can be seen in figure

4.4 as the two parts with high concentration of platelets and in figure 4.3 as the two peaks. The evolution of the PNWE, together with the evolution of umax and the CLF is plotted in figure4.7. Three phases in the equilibration are distinguishable: very soon the fluid equilibrates, then the CFL gets a constant value and finally the PNWE stays constant. Once the platelets are marginated towards the walls, to won’t return to the center of the vessel.

Figure 4.7: Equilibration of the fluid velocity umax, the cell-free layer (CFL) and the particle

near wall excess (PNWE) for a simulated vessel with channel width d = 100 µm, hematocrit H = 30% and Re ≈ 30.

By looking at the percentage of platelets near the wall as a function of time the speed of the margination is quantified, the results are plotted in figure4.8.

(30)

Figure 4.8: The margination of platelets for different hematocrit and Reynolds number in vessels with 100 × 200µm. The platelet near wall excess (PNWE) or platelet in margins percentage is defined as the volume percentage of platelets who ‘live’ inside an area of 5 µm around the platelet peaks at the walls. The error bars are 1σ.

4.4 Measuring viscosity and diffusivity in the micro-model

The viscosity of a fluid tells how good the fluid resists applied forces. The kinematic viscosity ν is calculated

ν = µ ρ =

σ

ρ ˙γ, (4.10)

with µ the dynamic viscosity, ρ the density of the fluid, σ the stress and ˙γ the shear rate.

To calculate the particle stress σ in order to calculate the viscosity in the simulation Batchelor’s method is used (30). For every LSP the stresses in the x and y direction and the positions are known. The average stress per particle LSP is

hσxyi = 1 V

X

(Fx· xy), (4.11)

where V is the volume of the particle, Fx is the force experienced by a LSP in the direction of the flow and xy is the y coordinate of the LSP. After calculation of the stresses per LSP the stress data is mapped onto a 1µm × the width of the vessel array.

(31)

The stress in Poiseuille flow can also be calculated theoretically. The stress in the centre of the flow should be equal to 0 and the wall stress is

σwall= gL

2 , (4.12)

where g is the body force applied to the fluid and L the width of the channel. The stress profile can be calculated by extrapolating between the wall stress and the stress at the centre of the channel. An example of such a theoretical stress profile and a measured stress profile for the simulated blood constituents is shown in figure 4.9. The shape of the theoretical and measured stress profiles is roughly the same, only at the edges they look different. The stress of the RBCs goes to 0 because there is a cell-free layer near the walls of the vessel and the change in steepness for the fluid stress near the wall is also due to the cell-free layer. In figure 4.10 then the summed stresses are noted. The irregular shape of the graph is due to the irregular presence of platelets. The red and the blue graphs are the same for the eye, because the blood plasma (fluid) stress is an order of magnitude smaller than the cell stresses and so the contribution of the plasma stress tot the total stress is negligible.

Figure 4.9: Left: Example of a theoretical stress profile in a channel with a width of 100 µm. Right: Example of a measured stress profile in a channel with a width of 100 µm. In red the stress on the RBCs, in yellow the stress on the platelets and in blue the stress on the blood plasma (fluid). The stress is still in lattice Boltzmann units, because the conversion to SI units is done later in the calculation for practical reasons.

Concerning suspensions it is interesting to look at the viscosity of the suspension relative to the viscosity of the fluid. The relative viscosity is defined as

νrel= νsuspension νfluid ≈ νcells+ νfluid νfluid . (4.13)

The local viscosities of the RBCs and platelets (calculated with Batchelor’s method) and for the fluid are plotted in figure4.11. There is an order of magnitude difference in the viscosities of the fluid and of the cells, which is probably due to cells. It is strange that the viscosities of the RBCs and the platelets are so similar. One other problem

(32)

Figure 4.10: The sum of the measured stresses from the right image in figure 4.9. In red the sum of the RBC and platelet stress, and in blue the sum of the RBC, platelet and blood plasma (fluid) stress.

with the viscosities is the big gap in the center. Trying to average and fit the quantities unfortunately didn’t fix the problem. As can be seen in figure 4.12 the profile gets smoother, but the gap stays.

Figure 4.11: Left: the local viscosity of the fluid in blue, of the platelets in yellow and of the RBCs in red. The gaps in the yellow line are due to an absence of platelets in those areas. Right: the relative viscosity (4.13). Both for a channel with a width of 100 µm, H = 40% and Re ≈ 1.5.

The diffusion of platelets is measured in the micro model by looking at the root mean square displacement of the platelets over time. It is not sufficient to look at the instan-taneous velocities of the platelets, because they are subject to a collision-free period just after the start of a measurement. Only after a time t all the particles have had a colli-sion and the measurement becomes significant. If one starts measuring immediately the resulting plot will show a unrealistic bend. Therefore it is best to track the displacement of platelets over time and calculate the diffusivity only in retrospect, for different boxes of particle paths. Figure 4.13 shows an example of the displacement in the direction perpendicular to the flow direction for the platelets. Displacement is shown here with all platelets starting at y = 0 while they have different y-coordinates. Next step will be

(33)

Figure 4.12: Left: viscosity calculated with theoretical stress profile (as in figure4.9, instead of measured stress profile) from one time step. Middle: viscosity calculated with theoretical stress profile, averaged over 10 time steps. Right: viscosity calculated with theoretical stress profile, averaged over 10 time steps, where all profiles where fitted to a 4th degree polynomial for all time steps.

to calculate the diffusivity D for every platelet and link it to their real y-coordinate to be able to plot a diffusion profile for the blood flow.

Figure 4.13: Typical displacement profile for platelets in the y-direction perpendicular to the flow direction in simulated blood.

(34)

Discussion

In this work the movement of platelets in blood is studied with a numerical model. Attempts have been made to answer a physical question by running different simulations and comparing the results with experimental results. Unfortunately numerical models always have certain drawbacks. It is not always possible to measure the same quantities in simulation results and in experimental results. For example local quantities such as velocity and stress are relatively easy to measure in a simulation, but harder to measure in an experiment. It is therefore not always possible to compare results.

At first we looked into the different equilibration times of blood. Work on this subject has been done for hard sphere suspensions, but not so much for soft sphere suspensions. It is found that the order of equilibration is: (i) blood plasma, (ii) RBCs and (iii) platelets. For platelets it can take orders of magnitude longer in comparison with the plasma and the RBCs to equilibrate. A good measure for equilibration time is difficult because ‘equilibrated’ is hard to define. Often a system keeps fluctuating or starts evolving slowly after a rapid change. With more and precise measurements it should be possible to define a function similar to (4.9) as proposed by Nott and Brady for blood or soft sphere suspensions in general, because we observed a clear difference in equilibration times when changing channel width, hematocrit or fluid velocity. Unfortunately we weren’t able to make a proper fit through our data.

As is shown in figure 4.8 the margination of platelets can be measured if one uses a slightly arbitrary definition, but the results are nonetheless interesting. Little work has been done on measuring the speed of the margination in the past but in this work a start has been made by looking at the dependencies of platelet margination. From the data it is clear that platelet margination depends on fluid velocity and hematocrit. The next step in this research would be trying to fit functions to the margination speed.

To find the driving force behind the margination of platelets an attempt has been made to measure the viscosity and the diffusivity in the micro-model. The viscosity profiles derived from the simulation data unfortunately show an (until now) unexplainable gap, which is possibly due to a numerical problem. As can be seen in section 4.4 the mea-surement of viscosities produced an order of magnitude difference in the values for the viscosities of the cells and the viscosity of the fluid. In (14) another blood flow model is discussed (which also makes use of the lattice Boltzmann method coupled to the im-mersed boundary method), and the same approach for calculating particle stresses as we use seems to work there, so probably we made a mistake in one of the calculations. A start has been made to measure the diffusivity of the platelets in the suspension but time didn’t permit to finish the measurements. Important to not for future researchers is that it is necessary to follow particles over time instead of using instantaneous velocities. This is because after the start of the measurement there is a collision-free time for all platelet, which would change the value of the diffusion D.

(35)

points representing different RBCs has not been taken into account. So the volumes used in this work deviate a little bit from the real volume of RBCs.

The lattice Boltzmann method is a numerical method to simulate fluid dynamics. There are numerous other techniques to resemble the Navier-Stokes equation such as finite difference or finite volume methods. It is interesting to compare different techniques and their outcomes. From a numerical point of view all methods have their drawbacks, for examples in implementing or in scaling the model. Different models have different ranges of validity (e.g. only for high or low Reynolds numbers), stability of accuracy. Because different methods are used for slightly different approaches it is difficult to compare articles and the work of research groups. A certain method for simulating blood may be accurate to describe blood flow in aneurysms, while another method is better for describing pulsatile flow.

A problem encountered with this thesis was both processor time and storage space. The simulations presented in this thesis used circa 4000 processor hours and fill up more than 600 GB. The aim was to run 8 simulations with different initial conditions for every used setups, because more simulations to average over, means better statistics (and thus smaller error bars). Starting with a study of storage space and runtime used by SuMO would help to better plan the simulations upfront.

Since studying the transport properties of platelets is interesting from a physiological point of view both more realistic setups and more realistic platelet properties should be used. The flow of platelets in aneurysms is already studied (19). Allowing platelets to stick together would allow them to form a thrombus. This way aneurysms can be simulated in more detail, but also various other diseases associated with platelets. One obvious extension of the model is to go from two to three dimensions. Some groups already managed to do this, and the group in Amsterdam is almost there.

(36)

Conclusion

With a 2D model for blood flow the physics of platelet movement in blood is studied, to get a better understanding of the speed of the margination and the driving forces behind the movement of the platelets. The well-known distribution of platelets and red blood cells in normal blood is resolved and the equilibration time for these distributions as well as for the flow velocities is measured. It is shown that the margination of the platelets is by far the last process in blood to equilibrate. The speed of the margination of platelets is compared for different model setups and is found to depend on the fluid velocity and the hematocrit of the blood. Lastly the local viscosity and diffusion of the blood have been measured to investigate the driving process of the platelet margination.

The lattice Boltzmann method coupled with the immersed boundary method seems to be a good way to simulate the flow of blood. The currently used 2D model should be expanded to 3D and biochemical effects in cell interactions should be taken in consid-eration or a more accurate simulation of blood. Eventually it will be possible to model the blood flow around disease-like structures as aneurysms.

Acknowledgements

I would like to thank dr. Alfons Hoekstra for setting me up with this project, prof. dr. Daniel Bonn and prof. dr. Morton Denn for their curiosity and the interesting discussions on suspension physics and dr. Rudolf Sprik and prof. dr. Ton van Leeuwen for their interest in my project and their willingness to serve as second and third readers. I thank SURFsara (www.surfsara.nl) for the support in using the Lisa Compute Cluster and everyone in the computational science research group for their kindness.

I’d like to thank dr. Eric Lorenz for his explanation of SuMO, his critical questions, helpful suggestions and his supervision in general. I am most grateful to Lampros Moun-trakis MSc. While working on this project my physical reasoning and computing skills improved a great deal because of him. I got a peak into how it is to actually do science and I also just had a great time discussing al kinds of subjects together. I not only ended up with a finished project, but also with a new friend.

(37)

References

1. I. Vaartjes, C. Koopman, I. van Dis, F. L. J. Visseren, and M. L. Bots. Hart- en vaatziekten in Nederland 2013, cijfers over leefstijl, risicofactoren, ziekte en sterfte. Hartstichting, 2013.

2. ´E. Guazzelli and J. F. Morris. A Physical Introduction to Suspension Dynamics. Cambridge University Press, 2011.

3. D. A. Fedosov, B. Caswell, A. S. Popel, and G. E. Karniadakis. Blood Flow and Cell-Free Layer in Microvessels. Microcirculation, 17(8):615–628, 2010.

4. H. Lei, D. A. Fedosov, B. Caswell, and G. E. Karniadakis. Blood flow in small tubes: quantifying the transition to the non-continuum regime. Journal of Fluid Mechanics, 722:214–239, 2013.

5. M. O. Bernabeu, R. W. Nash, D. Groen, H. B. Carver, J. Hetherington, T. Kru¨uger, and P.V. Coveney. Impact of blood rheology on wall shear stress in a model of the middle cerebral artery. Interface Focus, 3:20120094, 2013.

6. K. Vahidkhah, S. L. Diamond, and P. Bagchi. Platelet dynamics in three-dimensional simulation of whole blood. Biophysical Journal, 106(11):2529–2540, 2014.

7. L. Mountrakis, E. Lorenz, and A. G. Hoekstra. Validation of an efficient two-dimensional model for dense suspensions of red blood cells. Int. J. Mod. Phys. C, 25(11):1441005, 2014.

8. L. D. Landau and E. M. Lifschitz. Fluid Mechanics, volume 6 of Course of Theoretical Physics. Butterworth-Heinemann, 2nd edition, 1987.

9. F. M. White. Fluid Mechanics. McGraw-Hill, 4th edition, 1999.

10. P. K. Kundu, I. M. Cohen, and D. R. Dowling. Fluid Mechanics. Academic Press, 5th edition, 2012.

11. Y. C. Fung. Biomechanics: Mechanical Properties of Living Tissues. Springer, 2nd edition, 1993.

12. M. V. Kameneva, M. J. Watach, and H. S. Borovetz. Gender difference in theologic properties of blood and risk of cardiovascular diseases. Clinical Hemorheology and Microcirculation, 21(3):357–363, 1999.

13. A. M. Robertson, A. Sequeira, and M. V. Kameneva. Hemorheology. Hemodynamical Flows, 37:63–120, 2008.

14. T. Kr¨uger. Computer simulation study of collective phenomena in dense suspensions of red blood cells under shear. PhD thesis, Ruhr-Universit¨at Bochum, 2011.

(38)

16. R. F˚ahræus and T. Lindqvist. The viscosity of the blood in narrow capillary tubes. Am. J. Phys., 96(3):562–568, 1931.

17. P. A. Aarts, S. A. van den Broek, G. W. Prins, G. D. Kuiken, J. J. Sixma, and R. M. Heethaar. Blood platelets are concentrated near the wall and red blood cells, in the center in flowing blood. Arteriosclerosis, Thrombosis, and Vascular Biology, 8(6):819–824, 1988.

18. G. Segr´e and A. Silberberg. Behaviour of macroscopic rigid spheres in poiseuille flow. ii. experimental results and interpretation. J. Fluid Mech., 14:136–157, 1962. 19. L. Mountrakis, E. Lorenz, and A. G. Hoekstra. Where do the platelets go? A

simulation study of fully resolved blood flow through aneurysmal vessels. Interface Focus, 3(2), 2013.

20. X. Grandchamp, G Coupier, A. Srivastav, C. Minetti, and T. Podgorski. Lift and down-gradient shear-induced diffusion in red blood cell suspensions. Phys. Rev. Lett., 110:108101, 2013.

21. U. Frisch, B. Hasslacher, and Y. Pommeau. Lattice-Gas Automata for the Navier-Stokes Equation. Physical Review Letters, 56(14):1505+, 1986.

22. S. Succi. The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Numer-ical Mathematics and Scientific Computation). Oxford University Press, 2001. 23. A. G. Hoekstra. The lattice boltzmann method: Theory and embedding in micro-,

meso- and macroscopic theories. 2003.

24. T. Kr¨uger. Introduction to the immersed boundary method. 2011. http://www. timm-krueger.de/downloads/Krueger_Edmonton_IBM.pdf.

25. C. S. Peskin. The immersed boundary method. Acta Numerica, 11(-1):479–517, 2002.

26. P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collision processes in gases. i. small amplitude processes in charged and neutral one-component systems. Physical Review, 94(3):511?525, 1954.

27. P. R. Nott and J. F. Brady. Pressure-driven flow of suspensions: simulation and theory. Journal of Fluid Mechanics, 275:157–199, 1994.

28. A. Tokarev, G. Panasenko, and F. Ataullakhanov. Segregation of flowing blood: Mathematical description. Math. Model. Nat. Phenom., 6(5):281–319, 2011.

29. H. Zhao and E. S. G. Shaqfeh. Shear-induced platelet margination in a microchannel. Phys. Rev. E, 83:061924, 2011.

30. G. Batchelor. The stress system in a suspension of force-free particles. J. Fluid Mech., 41(3):545–570, 1970.

Shear-Induced Transport of Platelets in a Red Blood Cell Suspension