Incorporating concentration fields governed by advection-diffusion into Palabos-based CFD models in three dimensions using finite differences

(1)

Incorporating concentration fields governed by

advection-diffusion into Palabos-based CFD models in three

dimensions using finite differences

Report Bachelor Project Physics and Astronomy

Hugo Hoving January 31, 2020

Universiteit van Amsterdam

Faculteit der Natuurwetenschappen, Wiskunde en Informatica

Supervisor dr. G. Zavodsky

Examiner prof. dr. E.T. van Bavel

Institute UvA Informatics Institute

Project timeframe from 18-06-2019 to 20-12-2019

Study load 15 EC

(2)

This project concerns incorporating a concentration field governed by advection-diffusion into Palabos-based CFD models in three dimensions using finite differences. Data processors have been written in Palabos to approximate the derivatives in the advection-diffusion equation with central scheme finite differences in both 2nd and 4th order accuracy in the spatial reso-lution ∆x. Whenever central scheme would involve probing concentrations at locations inside or beyond the vessel wall, a one-sided finite difference is used in the direction opposite to the wall, accurate up to 1st order and 2nd order for the two data processors respectively.

Testing the model against analytical solutions for trivial cases of pure advection and pure diffusion indicates that the model reproduces the analytical solution qualitatively for most cases, but that numerical errors are generated. These errors seem to be larger for the advective part than for the diffusive part of the equation, for the values of the parameters chosen here, which correspond to the biologically relevant values. For pure advection of relatively steep jumps in concentration, central scheme leads to oscillations in the numerical solution. These seem to be mitigated effectively when upwind scheme is applied. Advection of smoother con-centration distributions seems to be captured by central scheme more accurately. However, contradictory to the finding for steeper concentration jumps, the distribution tested here is captured more accurately with central scheme than with upwind scheme.

Samenvatting voor de niet ingevoerde lezer

Dit project betreft het implementeren van een concentratieveld dat wordt beheerst door advectie en diffusie in een model van computational fluid dynamics (CFD) dat is geschreven met behulp van de vrij beschikbare verzameling algoritmes genaamd Palabos, in drie dimensies met behulp van de eindige differentiemethode. Het vakgebied van CFD is gericht op het modelleren van (stromingen van) vloeistoffen met behulp van (super)computers. Palabos is een specifieke vorm van CFD die gebruik maakt van de Lattice-Boltzmann Methode (LBM), welke methode wat geavanceerder is dan de methodes waarvan in de meeste commeri¨ele CFD-software gebruik wordt gemaakt. Hoewel de resultaten van dit project voor veel verschillende toepassingen van CFD relevant kunnen zijn, wordt met dit project de toepassing van CFD in de geneeskunde beoogd. Het gaat dan specifiek om het modelleren van bloed in de menselijke bloedsomloop met als doel om bijvoorbeeld hersenbloedingen en herseninfarcten te kunnen voorspellen. De eindige differentiemethode wordt in dit project gebruikt om de verschillende afgeleides die in de advectie-diffusievergelijking voorkomen te benaderen, zodat zij door de computer kunnen worden gehanteerd. De resultaten die in dit stuk gepresenteerd worden, laten zien dat het model in de meeste gevallen kwalitatief juist is, maar wel numerieke fouten genereert. Het model zal verder onderzocht moeten worden om te bepalen hoe deze fouten zo effectief mogelijk geminimaliseerd kunnen worden.

(5)

1 Introduction

1.1 Clinical decision making in the treatment of aneurysms

1.1.1 The current practise

Whenever someone is diagnosed with an aneurysm, a local dilatation of a blood vessel, the medical specialist in question is faced with the decision whether or not to proceed with surgery on the patient. Aneurysms pose an increased risk of dissection and rupture of the local blood vessel wall, potentially resulting in the disability or death of the patient. At the same time, surgery poses risks for many different complications, sometimes with equally severe conse-quences. If the dissection and rupture are sufficiently unlikely, conservative treatment or watchful waiting may suffice, rendering surgery unnecessary. Careful choice of treatment of aneurysms, therefore, has the potential to significantly decrease morbidity and mortality.

Both the current guidelines for treatment of aortic aneurysms and intracranial aneurysms base indications for surgery primarily on the diameter of the aneurysm, in accordance with the finding that greater diameter increases the risk of dissection and rupture. Secondary to that comes the assessment of other known risk factors, such as blood pressure, genetic disor-ders, family history of similar disease and patient age.1,2

As in all fields of medicine, this process of clinical decision making primarily consists of the systematic evaluation of a predefined set of risk factors and contraindications, resulting in a default indication for treatment from which the medical specialist will only rarely diverge. The existence and the relative significance that is to be attributed to the different determi-nants that are relevant to an indication are derived from medical research. The findings of most traditional medical research are statistic in nature. Whereas the laws of physics and chemistry are generally believed to be absolute, their interplay in the living human body is so complex that absolute laws of medicine can rarely be defined. Medical research, therefore, is traditionally limited to studying the physiology of large groups of individuals in relation to certain circumstances that are thought to effect it and drawing conclusions based on the relative occurrence of observed physiological effects within the group. Relying on statisti-cal conclusions like these, constantly adjusting and revising them in accordance with new findings, makes for the current state of medicine, which is undoubtedly impressive. For the

1_{Erbel, R. et al. 2014 ESC Guidelines on the diagnosis and treatment of aortic diseases. European Heart}

Journal 35, 2873-2916 (2014).

2

Gregory Thompson, B. et al. Guidelines for the Management of Patients With Unruptured Intracranial Aneurysms. Stroke 46, 2368-2400 (2015).

(6)

individual patient, however, these conclusions still often prove unreliable. It is in those cases that techniques capable of providing more patient-specific information can make the differ-ence. In general; the more relevant factors are mapped and monitored, the less is left to chance.

1.1.2 The promise of CFD

Computational fluid dynamics (CFD) is the field of computational science that concerns the modelling of fluids and flows. In a clinical context, CFD can be used to model the blood flow trough an aneurysm with computers and make estimates of certain parameters relating to the blood flow that can be relevant for the decision at hand. Most popular amongst those parameters in the field of CFD is the wall shear stress: the stress induced in the vessel wall parallel to the direction of flow, proportional to the velocity gradient perpendicular to the direction of flow. Both high and low values of shear stress are believed to contribute to unfa-vorable remodelling of the vessel wall. That way, for example, the decision not to surgically intervene in an aneurysm based on the observed aneurysm diameter can be adjusted by the finding that hemodynamics are unfavorable in the specific case.

Another such parameter can be the local concentration of platelets in the blood. Platelets are primarily responsible for the coagulation of blood. Whilst a desired effect for the healing of wounds, they can also induce thrombosis, when present under sufficiently high concentrations. Thrombosis can be fatal when it results in infarction of an artery and subsequent ischeamia of any vital organ. Aneurysms are often treated with endovascular stents or flow diverters. These devices are surgically inserted inside the vessel segment comprising the aneurysm to reinforce the wall or to favorably remodel local hemodynamics. The occurrence of a certain degree of thrombosis is then actually needed for the device to settle in place and conjoin with the adjacent vessel walls. However, the stability of this proces is fragile and the risk for unfavourable developments is always present. With methods of CFD that specifically ac-count for platelet concentration, predictions about local thromogenic conditions can be made, both before and after surgical intervention. These predications can be relevant for decisions regarding the placement of the device itself or for the postoperative prescription of medicines.

(7)

1.2 The present work: modelling platelet concentration

The present work concerns incorporating a concentration field governed by advection-diffusion into a Palabos-based model of blood flow in three dimensions. This concentration field, a three-dimensional scalar field, can be used to model the concentration of platelets in the blood as a function of space and time and make estimates of the local tromobogenic conditions. Additional concentration fields can be added analogously to account for different species of blood cells or any other constituent of the blood that may be of interest. The model into which this field is incorporated, was written by dr. G´abor Z´avodszky and is named Hemoflow. Hemoflow is designed to receive a three-dimensional model of a blood vessel, expose it to the boundary conditions representing human blood flow and resolve the physical parameters describing that flow throughout the model. Whereas Hemoflow treats the blood as a continuum fluid,3 like most CFD models, the addition of concentration fields draws information about individual constituents into view. This allows to monitor more information about the flow, but has the potential of increasing the accuracy of the modelled flow field as well. As the concentration of blood constituents influences the viscosity of the blood, the concentration field could be coupled back to the flow field, thereby increasing its accuracy.

3

With the side note that it concerns an LBM model here, which carries information on the particle distri-bution function as well.

(8)

2 Theoretical framework

2.1 The advection-diffusion equation

The advection-diffusion equation is a second order partial differential equation that expresses the change in concentration of a certain species of particles with respect to time at a given point in space in terms of two distinguishable processes capable of causing such change: ad-vection and diffusion - the two distinguishable modes of particle transport that are most common in nature.

In three dimensions, the advection-diffusion equation can be written compactly in the follow-ing way:

∂c

∂t = − ~∇ · (~vc) + ~∇ · (D ~∇c) (1)

Each of the quantities involved depend both on position (x, y, z) and time (t) - that excludes the derivative operators. Note that concentration c and the diffusivity D are scalars while the velocity ~v is a vector. Note that the inner product in the second term works between the two nablas, but that the derivatives in the first nabla must operate on D as well, so that D must reside in between the two nablas. Only if D is assumed constant, it can be taken in front of the second term.

Expanding (1) gives: = − ∂ ∂xx +ˆ ∂ ∂yy +ˆ ∂ ∂zzˆ · (vxx + vˆ yy + vˆ zz) c(x, y, z, t)ˆ + ∂ ∂xx +ˆ ∂ ∂yy +ˆ ∂ ∂zzˆ · D(x, y, z, t) ∂ ∂xx +ˆ ∂ ∂yy +ˆ ∂ ∂zzˆ c(x, y, z, t) (2)

Now dropping the dependencies and taking the inner products:

= − ∂ ∂x vxc − ∂ ∂y vyc − ∂ ∂z vzc + ∂ ∂x D ∂ ∂xc + ∂ ∂y D ∂ ∂yc + ∂ ∂z D ∂ ∂zc (3) = − vx ∂ ∂xc − c ∂ ∂xvx− vy ∂ ∂yc − c ∂ ∂yvy− vz ∂ ∂zc − c ∂ ∂zvz+ ∂ ∂xD ∂ ∂xc + D ∂2 ∂x2c + ∂ ∂yD ∂ ∂yc + D ∂2 ∂y2c + ∂ ∂zD ∂ ∂zc + D ∂2 ∂z2c (4)

On the right-hand side of (1), the first term represents advection and the second term repre-sents diffusion. Equivalently, in (2), (3) and (4), the first line reprerepre-sents advection and the second line represents diffusion.

(9)

Since the advection-diffusion equation is simply the sum of three identical equations for the three spatial dimensions, there is no loss of generality when confining to one dimension only. In the notation of (3), confining to the x-dimension gives the following comprehensible form of the advection-diffusion equation:

∂c ∂t = − ∂ ∂x vxc + ∂ ∂x D ∂ ∂xc (5) In the next sections, the terms in both parts of the equation are derived and discussed.

2.1.1 Fick’s Law of Diffusion

If the only mechanism responsible for the transport of particles is the random movement of those particles, the change in particle concentration is governed by Fick’s Law of Diffusion. Fick’s Law of Diffusion says that the amount of substance passing through a unit area with normal vector in a certain direction is proportional to minus the gradient of the concentration in that direction, with a constant of proportionality called the diffusivity D.

J = −D∂c

∂x (6)

Since the measure of amount of penetration per unit area per unit time is called flux, J is called the diffusion flux.

The law simply reflects that random movement of particles will result in a net displacement of particles from regions of higher concentration to regions of lower concentration at a rate proportional to how steeply the concentration differs between those regions. In the absence of other processes, diffusion will keep going until the concentration has become constant throughout the volume to which the particles are confined. Once the concentration is constant, there is no concentration gradient left to drive diffusion flux.

Note that any process of random movement of particles will by itself lead to diffusion in accordance to Fick’s Law. Diffusion does not require the presence of any forces that specifically push the particles from regions of higher concentration to regions of lower concentration. One can understand this by imagining a cloud of particles released into a finite volume and considering that the random movement of those particles will by itself result in the equilibrium situation of constant concentration of those particles throughout the entire volume. In specific situations, mechanisms causing particles to move in specific directions may be present in addition to the mechanism(s) causing random movement - the latter being only thermal diffusion in the simplest of cases. It is important, however, to strictly distinguish between those two types of mechanisms when accounting for them in models, as the next point will further illustrate.

(10)

Note that the diffusivity D is no more than a measure of ‘the tendency of particles to move with respect to each other’. The formulation of Fick’s Law requires the diffusivity to always be a positive number and does not foresee the preference of particles to move in a specific direc-tion. That is, the diffusivity may be decomposed into, for example, x-, y- and z-components, but the resulting diffusivities will still not allow to discriminate between positive and negative directions. Intuitively, such discrimination would be achieved by simply allowing the diffu-sivity to attain negative values, but doing so would precisely invert the description of physics provided by Fick’s Law. Negative diffusivity would imply that particles move from regions of lower concentration to regions of higher concentration, thereby increasing the concentration gradient that is already there, which would violate the second law of thermodynamics. It is therefore important to realise that any type of movement of particles that constitutes an absolute tendency of the particles to move with respect to each other can be absorbed into the diffusivity and forced into the workings of Fick’s Law, but that all direction specific types of diffusion have to be treated separately. These types of diffusion are also known as non-Fickian. These can be referred to as drifts.

The diffusion term in the advection diffusion equation follows from Fick’s Law in the fol-lowing way. Consider a cubic volume of finite size with two sides perpendicular to the x-axis, positioned within a cloud of particles whose only mode of transport is diffusion in the x-dimension. Mass conservation now requires that a change in concentration of particles inside the volume is caused by a difference between the amount of particles flowing into the vol-ume and the amount of particles flowing out of the volvol-ume. Since transport is only in the x-direction, both can be determined by only considering the sides of the cube perpendicular to the x-axis. If we define the diffusion flux along the x-axis to be positive in the positive x-direction, the rate of change of the average concentration ¯c of particles inside the volume equals the diffusion flux through the cube’s side situated at lowest x minus the diffusion flux through the cube’s side situated at highest x, divided by the distance between the two sides (see Appendix 2A for proof of the last step).

∂¯c ∂t = J(x1)_{− J}(x2) x2− x1 (7) Substituting (6) gives: ∂¯c ∂t = 1 x2− x1 −D(x1, t) ∂c ∂x _x=x 1 + D(x2, t) ∂c ∂x _x=x 2 (8) Now, in the limit of x2− x1 → 0, the average concentration on the left-hand side becomes the

ordinary concentration and the finite difference on the right-hand side becomes a derivative (see Appendix 2B for proof of the latter), so that we obtain:

∂c ∂t = ∂ ∂x D∂c ∂x (9) Which is the diffusive part of the advection-diffusion equation.

(11)

2.1.2 The advection term

If the only mechanism responsible for the transport of particles is the bulk movement of a medium in which the particles are immersed, the change in particle concentration is governed by advection - for which no specific law can be mentioned.

In stead of diffusion flux, we are now consider the advection flux:

J = ~vc (10)

This equation reflects that the bulk movement of a medium in which particles are immersed causes the particles to move along with the flow at the same velocity, so that the collection of particles is stationary when observed from a reference frame moving along at that velocity. Note that advection as represented by (10) will generally not hold in reality. For one thing, we know from Newton’s Second Law that a given force will provide a particle of smaller mass with higher acceleration than it will an otherwise identical particle of larger mass. Therefore, for example, solid particles of material of high mass density suspended in a fluid of material of low mass density, will, in principle, attain a slower velocity ~v0 than the velocity ~v of the fluid (10). This effect may in turn be cancelled, partially or entirely, due to the viscosity of the fluid. For example, if a fluid is very thick, meaning it has high viscosity, it may be capable of dragging along the solid particles suspended in it, even if its mass density is lower than that of the solid. Usually, multiple effects like these are at play to result in a net advective ~v0 6= ~v of the particles in question, which can be expressed in terms the relevant quantities involved and then substituted for ~v in the general expression of (10). In general, however, the bulk flow of the medium ~v will nonetheless be the most dominant amongst those quantities.

The advection term of the equation can be derived in analogy to the diffusion term. In stead of the diffusion flux, we now substitute the advection flux into (7):

∂¯c ∂t = 1 x2− x1 ~ v(x1, t)c(x1, t) − ~v(x2, t)c(x2, t) (11) Letting x2− x1→ 0: ∂c ∂t = − ∂ ∂x vxc (12) Which is the advective part of the advection-diffusion equation.

(12)

2.2 Solving the equation: the method of finite differences 2.2.1 The method of finite differences

The derivative of quantity f (x) with respect to quantity x is a measure for the amount of change of f (x) as a result of an amount of change of x. Since it is thus a measure of the difference between two quantities in a property shared by those quantities (a measure of the difference between f (x) and x in their tendency to change in relation to each other), it is natural to express it as a ratio. Following the first sentence, it should be the ratio of the change of f (x) to the change of x.

∆f (x)

∆x =

f (x + ∆x) − f (x)

∆x (13)

In the limit where the change of x approaches zero, by definition, the ratio becomes a deriva-tive: lim ∆x→0 f (x + ∆x) − f (x) ∆x ≡ df (x) dx (14)

Outside of this limit, the ratio, or specifically its numerator, is known as a finite difference: f (x + ∆x) − f (x)

∆x ≡

∆f (x)

∆x (15)

The derivative is obtained analytically, by applying predefined methods of differential calcu-lus. The derivative of a function with respect to a certain variable is often itself a function of that variable and can be evaluated over its entire domain by substituting the variable. The finite difference is obtained numerically, by observing the corresponding problem under consideration, reading off f (x), f (x + ∆x) and ∆x and substituting them into (15). The finite difference is a real scalar value and only refers to the part of the domain that is spanned by the difference, so that evaluating the entire domain requires repeating the same computational procedure for each of the differences into which the domain is chosen to be divided.

Note that the finite difference is essentially an average derivative over a certain piece of do-main and that, in that sense, the finite difference is a non-exact alternative to a genuine derivative.

The purpose of the present work is to incorporate a concentration field governed by an advection-diffusion equation into an existing model of computational fluid dynamics. This requires solving the advection diffusion equation for ∂c_∂t and that in turn requires computing multiple derivatives. The remarks above should suffice to make clear that derivatives are preferred over a finite differences. However, analytically computing the derivatives requires that ~v(x, y, z, t), c(x, y, z, t) and D(x, y, z, t) are known functions that we can manipulate with

(13)

differential calculus. This, unfortunately, is not the case here, meaning that we have no choice but to resort to finite differences. Finite differences will enable us to solve for ∂c_∂t, since we have the ability to access the values in the fields of ~v(x, y, z, t), c(x, y, z, t) and D(x, y, z, t) throughout the entire domain.

2.2.2 Central scheme and one-sided scheme

In practice, using the finite difference as expressed in (15) will give rise to a problem. f (x) will generally be confined to a finite domain in x, meaning that there will be a maximum value xmax, above which f (x) is not defined:

f (xmax+ ∆x) − f (xmax)

∆x =(undefined) (16)

The concentration field which is the subject of this paper is confined to a finite domain in space, since it models constituents of blood that are confined to a blood vessel. Probing, for example, c(xmax+ ∆x, y, z, t) in the process of computing ∂c(x,y,z,t)_∂x will not give the desired

result, because it amounts to examining the concentration at a location that lies either inside the wall or beyond the wall of the blood vessel, where there is no concentration to be exam-ined in the first place.

This problem is easily surpassed by taking finite difference in the opposite direction. This introduces the backward difference:

∇f (x)

∆x ≡

f (x) − f (x − ∆x)

∆x (17)

The nabla4∇ distinguishes it from (15), which can now be referred to as the forward difference: ∆f (x)

∆x ≡

f (x + ∆x) − f (x)

∆x (18)

This distinction leaves room for a third type of difference, the central difference: δf (x)

∆x ≡

f (x + ∆x) − f (x − ∆x)

2∆x (19)

In the present work, the central difference is used as the standard for approximating the derivatives and the forward and backward differences are used only when necessary due to the local presence of boundaries. The motivation for this choice will follow from the next subsection.

(14)

2.2.3 Accuracy considerations

An important difference between the forward and backward differences on the one hand and the central difference on the other, is that the error in approximating the corresponding derivative is proportional to ∆x when applying the forward and backward schemes, but pro-portional (∆x)2 when applying the central scheme. This means that the reduction of the error in the approximation of the derivative that results from a reduction ∆x, i.e. from an increase in the spatial resolution of the model, is greater for central scheme than for the one sided schemes. These relations follow from examining the Taylor expansions of the terms in the numerators of each of the three schemes. To see this, one first needs to understand Taylor the concept of expansions.

The Taylor expansion of a function f (x) about a certain point a within its domain is an infinite power series in (x − a) that equals f (x). The Taylor series of f (x) is:

f (x) = ∞ X n=1 (x − an) n! dnf (x) dxn x=a = f (a) + (x − a) 1! df (x) dx x=a+ (x − a)2 2! d2f (x) dx2 x=a +(x − a) 3 3! d3f (x) dx3 _x=a+ (x − a)4 4! d4f (x) dx4 _x=a+ (x − a)5 5! d5f (x) dx5 _x=a+ ... (20)

A function can be approximated by considering a finite amount of terms of its Taylor series. The approximation becomes more accurate as more terms are included. For example:

f (x) ≈ f (a) +(x − a) 1! df (x) dx _x=a+ (x − a)2 2! d2f (x) dx2 _x=a+ (x − a)3 3! d3f (x) dx3 _x=a (21)

It is common to refer to the accuracy of a Taylor expansion by stating the order of the first term that is not included in the series. The order of a term is simply the value of the exponent of (x − a) in that term. Since that term and all of the higher order terms to the left of it are not included in the expansion, their sum equals the error in the given approximation of f (x). Since the first of those terms will be largest, provided that the series is not evaluated too far away from a, the first term gives an indication of the size of the error of the expansion. As an example, (21) is a Taylor expansion of f (x) about a of fourth order.

We can now derive the difference in accuracy between the one-sided schemes and the central scheme. Consider the in the right-sided difference in (18) evaluated at x = a:

∆f (x) ∆x _x=a= f (a + ∆x) − f (a) ∆x (22)

(15)

Following (21), evaluating the Taylor expansion of f (x) about a at a yields: f (a) ≈ f (a) + (a − a) 1! df (x) dx x=a+ (a − a)2 2! d2f (x) dx2 x=a+ (a − a)3 3! d3f (x) dx3 x=a+ ... = f (a) (23)

Equivalently, evaluating the Taylor expansion of f (x) about a at a + ∆x yields: f (a + ∆x) ≈ f (a) +(a + ∆x − a) 1! df (x) dx x=a+ (a + ∆x − a)2 2! d2f (x) dx2 x=a+ (a + ∆x − a)3 3! d3f (x) dx3 x=a+ ... = f (a) + ∆xdf (x) dx x=a + (∆x)2 2 d2f (x) dx2 x=a + (∆x)3 6 d3f (x) dx3 x=a+ ... (24) Now substituting the results of (23) and (24) into (22):

∆f (x) ∆x x=a= f (a) + ∆xdf (x)_dx _x=a+ (∆x)2 2 d2_{f (x)} dx2 _x=a+ (∆x)3 6 d3_{f (x)} dx3 _x=a+ ... − f (a) ∆x = df (x) dx _x=a+ ∆x 2 d2f (x) dx2 _x=a+ (∆x)2 6 d3f (x) dx3 _x=a+ ... (25)

Dropping a in order to generalise to all x: ∆f (x) ∆x = df (x) dx + ∆x 2 d2_{f (x)} dx2 + (∆x)2 6 d3_{f (x)} dx3 + ... (26)

This says that the forward difference of f (x) equals the derivative of f (x) plus some power series in ∆x. This power series thus equals the error in the approximation of the derivative that is given by the forward difference. Since the first term of that series is of order one, we say that the forward difference in the form of (18) and (22) is first order accurate.

Performing the same procedure for the left-sided difference in (17) results in: ∇f (x) ∆x = df (x) dx − ∆x 2 d2f (x) dx2 + (∆x)2 6 d3f (x) dx3 + ... (27)

Which is first order accurate as well. However, for the central scheme in (19) we find: ∇f (x) δx = df (x) dx + (∆x)2 6 d3f (x) dx3 + ... (28)

So the central difference in the form of (19) is second order accurate. This motivates the choice to use central scheme whenever possible, and use one-sided schemes only whenever central scheme isn’t an option, due to the local presence of boundaries.

(16)

2.3 Computational fluid dynamics

Computational fluid dynamics (CFD) is the field of computational science that concerns the modelling of fluids and flows. It seeks to develop numerical solving algorithms for computers that can be used to apply the existing physics of fluid dynamics to flows that are too com-plex to study analytically. Since computers can be many orders of magnitude faster than humans at executing numerical tasks and are entirely flawless in doing so, the computer is an extremely powerful tool in this context and forms an indispensable part of CFD.

One of the most ubiquitous equations in (computational) fluid dynamics is the Navier-Stokes equation (NSE): ρ ∂~v ∂t + ~v · ∇~v = −∇p + µ∇2~v + 1 3µ∇(∇ · ~v) + ~F (29)

With ρ the density of the fluid, ~v its flow velocity, p the pressure in the fluid and F the external body forces acting on the fluid. This form of the equation appertains to the flow of a compressible medium and is thus said to govern ’compressible flow’. Besides, it assumes constant (shear) viscosity (or dynamic viscosity) µ and zero volume viscosity ζ. The expres-sion between parentheses on the left-hand side equals the material derivative of ~v and may be expressed more compactly as D~_Dtv. The NSE follows from applying the law of conserva-tion of momentum (Newton’s classical mechanics) to a volume of fluid. A complete derivaconserva-tion of the equation will not be given here, as this does not serve the main purpose of this project.5

Problems of fluid dynamics typically involve finding the flow velocity ~v of the fluid under consideration at different points inside the fluid. Since this velocity is governed by the NSE, this task amounts to solving the NSE. However, analytic solutions to the NSE can only be found for the simplest types of flow. In most cases, therefore, numerical approaches are needed to approximate ~v. This is where CFD comes into play.

Traditional forms of CFD rely on the method of finite volumes. Roughly speaking, this can be regarded as a three-dimensional analog to the method finite differences discussed in the previous section. The strength of the method of finite volumes is that it is relatively simple in principle. Amongst its weaknesses in the context of solving the NSE, however, is that the approximation schemes can become rather lengthy and cumbersome. The gradient of the velocity ∇~v, for example, is a second order tensor that contains nine derivatives: each of the three spatial components of ~v has to be differentiated with respect to each of the three

5

(17)

dimensions. Taking the inner product with ~v subsequently requires nine multiplications. The matter is complicated even more by the circumstance that the mapping of ~v to ~v · ∇~v is not linear.

2.3.1 The Lattice Boltzmann method

The Lattice Boltzmann method (LBM) is a method of CFD that does not obtain ~v by directly solving the NSE, but by utilising of the particle distribution function f (x, y, z, ξx, ξy, ξz, t) of

the flow. The particle distribution function equals the number of particles at position x, y, z that have velocity ~ξ = ξxx + ξˆ yy + ξˆ zz at time t. The function thus extends the concept ofˆ

density to velocity space, so that density becomes velocity-specific.

The integral of f over velocity space equals the density ρ: ρ(x, y, z, t) =

Z

f (x, y, z, ξx, ξy, ξz, t) d3ξ (30)

The integral of ~ξf over velocity space equals the product of the flow velocity and the density, ~vρ, known as the momentum density:

~v(x, y, z, t)ρ(x, y, z, t) = Z

~

ξf (x, y, z, ξx, ξy, ξz, t) d3ξ (31)

ρ and ~vρ are also known as the zeroth and first order moments of f , respectively. In order to calculate ~v in this way, which is what we are after, we thus first need to obtain the particle distribution function. The distribution function is governed by the Boltzmann equation:

∂f ∂t + (~ξ · ∇)f + ~ F ρ · ∂f ∂~ξ = Ω(f ) (32)

With Ω(f ) the so called collision operator. The distribution function is obtained by solving that equation.

The problem with this is that the Boltzmann equation is, as such, not easier to solve than the NSE. For one thing, f depends on no less than seven variables: x, y, z, ξx, ξy, ξz and

t. The trick is that there exists a numerical approach to solving the Boltzmann equation, which, although a bit more involved in principle, turns out relatively simple in practise. This approach relies on the use of so called Hermite polynomials. As the entire derivation is too lengthy to set out in this text and would not serve its purpose, I refer to chapter 3 of The Lattice Boltzmann Method, Graduate Texts in Physics by T. Kr¨uger et al.. The result which is of relevance here, is the lattice Boltzmann equation which is obtained in that way:

(18)

The particle distribution function f has now become the discrete-velocity particle distribution function fi. Both time and space have been discretised. This means that the velocities are now

also constrained to a domain of certain discrete values. They are denoted ci = (cix, ciy, ciz),

with i the index labeling the different allowed velocities. At each point in discretised space, the distribution function now consists of a collection of values for fi, one for each of the

discrete velocities i, that each equal the density of particles with that velocity i at that point. Those values of fi are also referred to as the particle populations at that point.

The moments of fi equal sums of the function evaluated at the different discrete points of its

domain, in stead of integrals over a continuous domain as was the case with f : ρ(x, y, z, t) =X i f (x, y, z, ξx, ξy, ξz, t)ξ (34) ~v(x, y, z, t)ρ(x, y, z, t) =X i ~ ξf (x, y, z, ξx, ξy, ξz, t) (35)

When we substitute the BGK collision operator Ωi(f ) = −1_τ(fi− f_ieq)∆t, the most commonly

used collision operator in LBM, we obtain the lattice BGK equation: fi(x + ci∆t, t + ∆t) = fi(x, t) −

∆t

τ (fi(x, t) − f

eq

i (x, t)) (36)

Now, the LBM algorithm consists of two steps: collision and streaming. First, at each of the discrete points of the spatial domain, the distribution function is altered as prescribed by collision operator, locally at that point:

f_i∗(x, t) = fi(x, t) −

∆t

τ (fi(x, t) − f

eq

i (x, t)) (37)

The collision operator redistributes the particles at each point amongst the different particle populations, modelling the redistribution of velocities amongst particles that occurs during physical collisions, and ensures that this redistribution is in accordance with the physical requirement that the system progressively moves towards an equilibrium state over time. Second, each of the new populations f_i∗ is substituted for a population at a neighbouring point, as determined by the direction and magnitude of the corresponding velocity ci (the

new distribution is said to be ’streamed’ to the neighbouring points):

fi(x + ci∆t, t + ∆t) = fi∗(x, t) (38)

The net effect of the two steps (37) and (38) now equals one time step in accordance with the lattice Boltzmann equation (36).

In reality, collision and streaming are, of course, practically indistinguishable parts of the same continuous process. That they appear separately in the algorithm, is just the conse-quence of the choice to compute the effect of the collisions on the velocities of the particles

(19)

on the one hand and on the positions of the particles on the other, consecutively in stead of simultaneously. The collision step requires computing feq and therefore requires comput-ing both ρ and ~u. Performing the collision step locally means that the computer does not need to communicate with neighbouring parts of the domain to perform those computations. This allows for more computational efficiency and specifically enables the algorithm to be adequately parallelised. The streaming step contains all the non-local behaviour at once, but that one is a lot less computationally intensive than collision, as it is just the plain exchange of one species of data between different points in the domain.

2.3.2 The Palabos library

Palabos is a freely available library (a collection of pre-written algorithms) for computational fluid dynamics based on the Lattice Boltzmann method. Palabos is written in the program-ming language C++. Palabos is an acronym for ‘Parallel Lattice Boltzmann Solver’.

This section briefly discusses the characteristic aspects of Palabos that are relevant to the present work. For a more complete and detailed description, I refer to the Palabos User Guide.6

2.3.2.1 Parallel computing

Parallel computing refers to the process of dividing a computational task into smaller com-putations and performing those comcom-putations simultaneously. Parallel computing is distin-guished from serial computing, in which a computational task is divided into smaller com-putations that are performed consecutively. A CFD model is, like any other model, usually developed with the aim to be as accurate as possible. One might specifically imagine that the results of CFD simulations will have to meet certain accuracy requirements if they are ever to be used in the context of medical diagnostics - as is the aim of the code corresponding to the present paper. Attaining higher accuracy typically requires performing more intensive (or extensive) computations. At the same time, the computing capacities of all computational cores are finite, both in terms of the quantity of computations they can handle at a time and the speed at which those computations can be performed. Performing a typical CFD simu-lation in serial at a reasonable accuracy will, depending on the precise size of the simulated domain, quickly take up unreasonably large amounts of time. Because the serial approach involves performing all the computations consecutively on one core, it can only exploit the capacities of one computational core at a time. The parallel approach, to the contrary, divides

6

(20)

the problem over many different cores and lets all of them work on a piece of the problem simultaneously, so that the parallel approach can employ the combined capacity of as many cores as there are available. This means that parallel computing has the potential to decrease the runtime of a program dramatically.

Working in parallel does, of course, require that the code is written in a way that allows for this division and distribution of the problem. This logistic task gives rise to some specific computational structures and concepts. In Palabos, these are primarily the BlockXD data structures and the concept of data processors.

2.3.2.2 The BlockXD data structures

In Palabos, all data that is relevant to the user in the context of performing simulations, is stored in data structures that, most generally, resemble arrays or matrices. The most fun-damental types of these arrays are called Block2D and Block3D, for 2D and 3D simulations respectively. For convenience, we refer to them as BlockXD in general. Specific blocks are given more specific names depending on the type of data they are designed to store. The particle populations fi, for example, are stored in blocks named BlockLatticeXD. Regular

scalar valued quantities, such as the concentration modeled in this paper, are stored in data structures named ScalarFieldXD.

Palabos discriminates between two differents species of blocks: the AtomicBlock and the MultiBlock. The AtomicBlock is just the corresponding array in its general form. The MultiBlock is essentially a collection of multiple AtomicBlocks grouped together. Accord-ingly, BlockLatticeXD and ScalarFieldXD have counterparts named MultiBlockLatticeXD and MultiScalarFieldXD. Palabos is programmed to automatically divide MultiBlocks into suitable smaller atomic blocks to cope with irregularities in the domain and to ensure that the problem can be adequately parallelised. Therefore, for most applications of Palabos, the user will use only MultiBlocks.

2.3.2.3 The concept of data processors

One of the most basic concepts and most common tools in any programming language is the for loop. With the for loop, the instructions of an algorithm are repeatedly applied to each element of a data structure, such as a list or an array. This way, the for loop ensures that the algorithm is applied to the entire structure, in a systematic way. The for loop is a most typical example of a serial procedure, as it treats each element in the data structure

(21)

consecutively and separately. This also means that the for loop is, as such, incompatible with a parallel computing approach. In practice, when an algorithm is applied to, for example, a MultiScalarField3D by directly looping that algorithm over the scalar field, Palabos will either crash and report an error or proceed and produce erroneous output. The instructions implied by the for loop disagree with the requirement that the code be parallelised. The for loop could itself be parallelised by rewriting it into smaller for loops tailored to each of the AtomicBlocks into which the MultiBlock is divided and making sure those for loops are assigned to each of those respective AtomicBlocks, but Palabos is not equipped to do this automatically. In stead, the user is required to incorporate the for loop and the algorithm in question into a specific type of structure called a data processor. The data processor contains all sorts of additional instructions to the computer besides the for loop and the algorithm. When these instructions are provided and configured correctly, the data processor can be applied directly to a MultiBlock or an AtomicBlock and should give the desired result.

(22)

3 Methods

3.1 Incorporating the concentration field

3.1.1 Using pointers to the fields

The advection-diffusion algorithm does not manipulate the concentration field directly. In stead, a pointer to the concentration field is used, so that the concentration field is ma-nipulated by the advection-diffusion algorithm after dereferencing that pointer. Declaring a pointer to the field generates a single address for the concentration field which can be com-municated with each of the different computational cores onto which the advection-diffusion algorithm is allocated. As a consequence, the different computational cores can access the concentration field simultaneously and independently, so that the code will work when ran in parallel.

3.1.2 Using different fields for reading and writing

It might surprise that the data processor should receive two ScalarField3D data structures. One of them should surely be the concentration field, but that leaves the other. In fact, one of the scalar fields represents the concentration field at the previous time step, and is used to read from, and the other scalar field represents the concentration field at the second previous time step, and is used to write to. Thus - if we agree that the data processor oper-ates in between the previous and the next time step - the values of the concentration needed for computing the concentration corresponding to the next time step are read off from the concentration field corresponding to the previous time step and the values of the concentra-tion corresponding to the next time step are, once computed, written over the values of the concentration field corresponding to the second previous time step. At that moment, the concentration field corresponding to the second previous time step is no longer relevant.

Using separate fields for reading and writing is required, because, for each subspaces of the spatial domain that result from parallelisation, the advection-diffusion algorithm is swept through the domain in a serial procedure with the use of a for loop. The finite difference computations performed in the course of updating the concentration field at each lattice node in accordance with advection-diffusion involve reading off the values of the concentra-tion field in neighbouring lattice nodes. If one and the same concentraconcentra-tion field were used for reading and writing, the computations at each lattice node, apart from the very first one, would read off concentrations of neighbouring nodes that have, in part, already been updated. This would result some distorted variant of the intended advection-diffusion behaviour. In

(23)

reality, of course, the combined process of advection and diffusion occurs continuously and equally simultaneously throughout the entire domain. Therefore, one time step of modelled advection-diffusion should impose upon all lattice nodes a change relative to the previous time step, based on the previous time step.

(24)

3.2 Incorporating the advection-diffusion equation 3.2.1 The advection-diffusion algorithm

For each spatial dimension, we have: ∂c ∂t = − ∂ ∂x vxc + ∂ ∂x D ∂ ∂xc = − vx ∂ ∂xc − c ∂ ∂xvx+ ∂ ∂xD ∂ ∂xc + D ∂2 ∂x2c (39)

In general, we use the approximation:

cnew = cold + ∂c ∂t old∆t c(xi, ti+ ∆t) = c(xi, ti) + ∂c ∂t (xi,ti) ∆t (40)

So this gives, for each dimension:

c(xi, ti+ ∆t) = c(xi, ti) + − v_x ∂ ∂xc − c ∂ ∂xvx | {z } x part ad + ∂ ∂xD ∂ ∂xc + D ∂2 ∂x2c | {z } x part dif _(x i,ti) ∆t (41)

Approximating the derivatives with central finite differences of second order gives, for each dimension, at each time step:

x part ad = − vx(xi) c(xi+ ∆x) − c(xi− ∆x) 2∆x − c(xi) vx(xi+ ∆x) − vx(xi− ∆x) 2∆x x part dif =D(xi+ ∆x) − D(xi− ∆x) 2∆x c(xi+ ∆x) − c(xi− ∆x) 2∆x + D(xi)(c(xi+ ∆x) − 2c(xi) + c(xi− ∆x)) (42)

In the code, these become:

T x_part_ad = −v_x ∗ (c_xr − c_xl) ∗ .5 − c ∗ (v_x_xr − v_x_xl) ∗ .5

T x_part_dif = (d_x_xr − d_x_xl) ∗ (c_xr − c_xl) ∗ .25 + d_x ∗ (c_xr − 2. ∗ c + c_xl) (43) In the code, then, (40) simply becomes:

c_new = c + (x_part_ad + y_part_ad + z_part_ad + x_part_dif + y_part_dif + z_part_dif)

(44)

Which sums up one step of advection-diffusion - the core of the algorithm contained in the advection-diffusion data processor.

(25)

Four aspects of the result of (43) may benefit from some explanation.

First, the subscripts xr and xl denote the right and the left neighbouring lattice nodes respec-tively. For example, c_xl denotes the concentration at one lattice node to the left, v_x_xr denotes the x component of the velocity at one lattice note to the right and d_x_xr denotes the ‘x-component of the diffusivity’ at one lattice note to the right. For the sake of readability of the code, the corresponding directions in the y- and z-dimensions can be labeled with r and l as well, defining r and l as the positive and negative directions respectively.

Second, the diffusivity may now appear to be decomposed into distinct components along the different dimensions. This may seem in contradiction to the earlier presentation of D as a scalar. In the present work, D is indeed treated as a scalar, but its computation is made direction specific, so that, at each point, it attains a different value for each of the three spatial dimensions. This does not mean that D is treated as a vector, for that would suggest that the net magnitude and direction of the diffusivity are known and that d_x, d_y and d_z represent its three perpendicular vector components. This is not the case. The diffusivity is simply computed in the x-, y- and z-dimensions and the resulting values depend on the shear rate measured in that direction. It is unlikely that these three components will span the vector representing the net diffusivity. The discussion contains further comments on this issue.

Third, when using the LBM, it is common to work in so-called lattice units. Converting from any set of units to lattice units is done by setting both the temporal resolution ∆t and the spatial resolution ∆x equal to 1 and scaling all units related to time and space accordingly. In SI units, both the spatial and the temporal resolution will typically be very small numbers, sometimes with many digits after the decimal point. Multiplication with and division by 1 poses less computational load on computers than multiplication with and division by any other number. Since ∆x and ∆t are involved in the computations at each lattice node, at each spatial step, at each temporal step, replacing them by 1 results in significant gain in computational efficiency. The entire code can be executed in this way, and one can simply convert back to SI units before saving the output. Section 3.3 contains a brief explanation of how this is done in the present work.

Fourth, current Intel processors are slightly quicker at multiplication than at division. There-fore, multiplying with 0.5 is preferred over dividing by 2.

The required velocities in (43) can be easily obtained in Palabos using predefined functions that apply to the BlockLattice. The required values of the concentration field are queried easily as well - we manually put them there in the first place. This leaves the diffusivity, which is defined and computed manually as well, but locally inside the data-processor.

(26)

3.2.2 Modelling diffusivity

The present work uses the following model for the diffusivity suggested by dr. G. Z´avodszky:7

D = C0σ( ˙γ)| ˙γ|Φ (45)

This model agrees with current standards in literature for the diffusivity in shear flows. It differs from the version given in the cited article by including the absolute value of the shear rate, in stead of the shear rate as such. This is to ensure that D attains only positive values, as is required by definition.

C0 is a parameter related to the shape of the particles subject to the diffusion, ˙γ is the shear rate of the flow, σ( ˙γ) is the collision cross-section of the particles subject to the diffusion, based on mutual collisions or collisions with other species of particles present in the suspen-sion, and Φ is the volume fraction of the species of particles responsible for the diffususpen-sion, relative to the total volume of the suspension.8

In the present work, Φ is the volume fraction of red blood cells - also called the hemat-ocrit. In reality, hematocrit is a function of position in the plane perpendicular to the flow. In laminar flow, for example, hematocrit will be a strong function of the perpendicular dis-tance r to the longitudinal axis of the blood vessel, with maximum concentration at r = 0. Accounting for this dependency of Φ on position would require incorporating a second con-centration field to keep track of red blood cell concon-centration. In the present work, Φ is simply taken to be 0.35 (= 35%), which is a reasonable value of the average hematocrit.

The collision cross-section σ1,2 of two hard spheres or radii r1 and r2 is defined as:

σ1,2 = π(r1+ r2)2 (46)

In the present work, σ depends on the shear rate ˙γ, because cells tend to be stretched when exposed flows of nonzero shear.9 Therefore, in the present work, we have:

σ( ˙γ) = π h r1( ˙γ) + r2( ˙γ) i2 (47) 7

Z´avodszky, G. et al. Red blood cell and platelet diffusivity and margination in the presence of cross-stream

gradients in blood flows. Phys. Fluids 31, 031903 (2019).

8

The species of particles subject to the diffusion may differ from the species of particles responsible for the diffusion. In the present work, for example, the particles subject to the diffusion are platelets, but since their diffusivity is caused mainly by their collisions with red blood cells, as red blood cells are more than an order of magnitude more occurrent in the blood than the platelets themselves, the particles responsible for the diffusion are taken to be the red blood cells.

(27)

The present work concerns the diffusivity of platelets (PLTs) in blood. Red bloods cells (RBCs) are the most dominant species of cells in the blood. Therefore, the present work uses for the diffusivity solely the collision cross-section corresponding to PLT-RBC collisions. The radius of platelets is assumed independent of shear rate and set equal to 1.2 µm. The radius of RBCs is assumed dependent of shear rate. Following the work of G. Z´avodszky et al., we use:10,11 r1 = rPLT= 1.2 µm , r2 = rRBC ˙γ = 3.92 µm 1 1.6e − 0.0013 | ˙γ|₊0.6 1.6 (48) Which gives: σ( ˙γ) = π " 1.2 + 3.92 1 1.6e − 0.0013 | ˙γ|₊0.6 1.6 #2 (49) The parameter C0 is taken to be, again following G. Z´avodszky et al.,12 in accordance with aPLT= 1.02 µm2:

C_PLT0 = aPLT σ( ˙γ = 0) =

1.02

π(1.2 + 3.92)2 (50)

This leaves the shear rate ˙γ. The shear rate is defined as the derivative of the velocity with respect to position in a direction perpendicular to the velocity. The shear rate can be modeled in different ways. In the present work, the shear rate is computed separately for each of the three spatial dimensions and is set equal to sum of the two derivatives of the velocity with respect to position evaluated in the two directions perpendicular to that dimension:

˙γ(x)= ∂vy ∂x + ∂vz ∂x , ˙γ (y) ₌ ∂vx ∂y + ∂vz ∂y , ˙γ (z)₌ ∂vx ∂z + ∂vy ∂z (51)

The notation ˙γ(x), for example, indicates that the shear rate is evaluated in the x-direction. Note that this is not the same as ˙γx, the x-component of the (net) shear rate. The discussion

contains further remarks on alternative calculations of the shear rate in this context.

10_Ibid.

11_{This again differs from the cited work in including absolute value signs around the shear rate.}

12_{Ibid. The cited article states that a}

(28)

In view of (43), (45) and (51), the values of ˙γ required for computation of the diffusivity become in the code, approximating the derivatives with finite differences of second order, for each dimension:

s_x = abs((v_y_xr − v_y_xl) ∗ .5) + abs((v_z_xr − v_z_xl) ∗ .5) s_x_xl = abs((v_y − v_y_xll) ∗ .5) + abs((v_z − v_z_xll) ∗ .5) s_x_xr = abs((v_y_xrr − v_y) ∗ .5) + abs((v_z_xrr − v_z) ∗ .5)

(52)

Note that the computations of the shear rate introduce the need for probing not only the neighbouring lattice nodes, but the next neighbouring lattice nodes as well. For example, v_z_xll is the z-component of the velocity at two lattice nodes to the left (i.e., at −2∆x in x relative to the current node).

In the end, the computations of the diffusivity for the x-dimension become in the code: T &d_x = C_PLTRBC() ∗ cc_PLTRBC(s_x) ∗ s_x ∗ Phi_RBC;

T &d_x_xl = C_PLTRBC() ∗ cc_PLTRBC(s_x_xl) ∗ s_x_xl ∗ Phi_RBC; T &d_x_xr = C_PLTRBC() ∗ cc_PLTRBC(s_x_xr) ∗ s_x_xr ∗ Phi_RBC;

(53)

Where C_PLTRBC() is a function that computes C0, cc_PLTRBC() is a function that computes σ( ˙γ) and Phi_RBC is simply Φ.

3.2.3 Treating walls

The concentration field is, of course, confined to the same spatial domain as the blood flow under consideration, as it represents the concentration of a constituent of the blood. There-fore, the advection-diffusion data processor should loop the advection-diffusion algorithm only over the part of the spatial domain that is occupied by the fluid. Any node that lies out-side the fluid domain is called a boundary node, to be distinguished from the fluid nodes. Whether a node is a boundary node or not can be easily queried in Palabos with a predefined function. The above requirement is therefore easily met by including a line at the very top of the algorithm that checks whether the node in question is a boundary node. If so, it skips to the next node in line before the rest of the algorithm is executed.

Computing finite differences requires reading values from neighbouring lattice nodes and for fluid nodes adjacent the vessel wall some of those will be boundary nodes as well. The con-centration and velocity at the boundary nodes will equal 0 and using them in the prescribed computations will lead to erroneous results. For example, when observing the upper line of (43) for a lattice node - let’s call it A - directly adjacent to a wall on its left side, the zero

(29)

concentration and velocity at lattice node B directly to the left of it will imply a flow towards B away from A. As a consequence, the concentration will decrease at A, whilst no concentra-tion increase can ever happen at B, since B lies outside of the fluid domain. To avoid errors like these, the code contains a procedure to check whether any of the neighbouring nodes used by the central scheme is a boundary node. If so, the central scheme is replaced with a one-sided scheme of first order, taken in the direction opposite to that in which the boundary is encountered.

For each dimension, fully executing the core of the algorithm in central scheme requires, in view of the previous paragraph, that the present node is a fluid node and that on each side both the neighbour and next neighbour are fluid nodes. If either the right neighbour or right next neighbour is a boundary node and both the left neighbour and left next neighbour are fluid nodes, the finite differences in x_part_ad and x_part_dif are replaced by left-sided schemes of first order. Also, since d_x still occurs in those schemes and d_x by itself requires, through its dependence on s_x, probing the right neighbour as well, the central scheme used for the calculation of s_x in d_x is replaced with a left-sided scheme as well:

if((lattice.get(iX + 1, iY, iZ).getDynamics().isBoundary() || lattice.get(iX + 2, iY, iZ).getDynamics().isBoundary()) &&

lattice.get(iX − 1, iY, iZ).getDynamics().isBoundary() == false && lattice.get(iX − 2, iY, iZ).getDynamics().isBoundary() == false){ T &c_xll = concentrationField.get(iX − 2, iY, iZ);

s_x = abs(v_y − v_y_xl) + abs(v_z − v_z_xl); d_x = C_PLTRBC() ∗ cc_PLTRBC(s_x) ∗ s_x ∗ Phi_RBC; x_part_ad = −v_x ∗ (c − c_xl) − c ∗ (v_x − v_x_xl);

x_part_dif = (d_x − d_x_xl) ∗ (c − c_xl) + d_x ∗ (c − 2. ∗ c_xl + c_xll); }

(54)

If either the left neighbour or left next neighbour is a boundary node and both the right neighbour and right next neighbour are fluid nodes, the finite differences in x_part_ad and x_part_dif are replaced by right-sided schemes of first order. This again includes replacing the finite difference in d_x:

(30)

if((lattice.get(iX − 1, iY, iZ).getDynamics().isBoundary() || lattice.get(iX − 2, iY, iZ).getDynamics().isBoundary()) &&

lattice.get(iX + 1, iY, iZ).getDynamics().isBoundary() == false && lattice.get(iX + 2, iY, iZ).getDynamics().isBoundary() == false){ T &c_xrr = concentrationField.get(iX + 2, iY, iZ);

s_x = abs(v_y_xr − v_y) + abs(v_z_xr − v_z); d_x = C_PLTRBC() ∗ cc_PLTRBC(s_x) ∗ s_x ∗ Phi_RBC; x_part_ad = −v_x ∗ (c_xr − c) − c ∗ (v_x_xr − v_x);

x_part_dif = (d_x_xr − d_x) ∗ (c_xr − c) + d_x ∗ (c_xrr − 2. ∗ c_xr + c); }

(55)

Finally, if either the right neighbour or the right next neighbour and either the left neighbour or the left next neighbour are boundary nodes, both x_part_ad and x_part_dif are set equal to zero:

if((lattice.get(iX + 1, iY, iZ).getDynamics().isBoundary() || lattice.get(iX + 2, iY, iZ).getDynamics().isBoundary()) && (lattice.get(iX − 1, iY, iZ).getDynamics().isBoundary()||

lattice.get(iX − 2, iY, iZ).getDynamics().isBoundary())){ x_part_ad= 0.;

x_part_dif= 0.;}

(56)

3.2.4 Two data processors: 2nd and 4th order accurate

Besides this 2nd order accurate data processor, a 4th order accurate advection-diffusion data processor is built as well. This way, the results of the two can be compared, so that conclusions may be drawn on the effects of the order of the approximations on the performance of the model. The 4th order data processor is structured analogously to the 2nd order one, but now all central schemes are 4th order accurate and all one-sided schemes are 2nd order accurate. A relevant difference between the two data processors is that the 2nd order one requires 5 consecutive fluid nodes for central scheme and 3 for one-sided scheme, whilst the 4th order one requires 9 consecutive fluid nodes for central scheme and 5 for one-sided scheme.

(31)

3.3 Unit conversion

The operations required to convert a quantity from SI units to lattice units follow from the dimensional analysis of that quantity. For the concentration c, the diffusivity D and the velocity v this entails:

[c]SI= 1/m3 −→ cLBM= (∆x)3 cSI

[D]SI= m2/s −→ DLBM= (∆t/(∆x)2) DSI

[v]SI= m/s −→ vLBM= (∆t/∆x) vSI

(57)

The concentration field c is spawned by the code at the beginning of the simulation and its initial values can be configured in any desired way. If the choice to configure the concentration field is made with reference to certain standards expressed in SI units, the values loaded into the concentration field are to be converted to LBM units before loading them into the initial concentration field. They can then be converted back to SI units before saving, using the inverse operation.

The diffusivity D is computed locally inside the data processor based on C0, σ( ˙γ), ˙γ and Φ. The parameters C0 and Φ are dimensionless numbers and are therefore untouched by unit conversion. The parameter σ( ˙γ) is computed as a function of rPLT and rRBC( ˙γ). The radii

are converted to lattice units via rLBM= (1/∆x) rSI. Finally, ˙γ, the last of the parameters

required to compute D, is purely a function of different values of the velocity v.

The velocity v is obtained from the BlockLattice with a predefined Palabos function. Veloc-ities received in this way are in lattice units by default. This means that no unit conversion of v is required. Consequently, since ˙γ is solely a function of v obtained in this way, no conversion of ˙γ is required either.

Do note, however, that the shear rate in the definition of the RBC radius must remain in SI units, as it functions purely as a parameter in the exponential to scale this radius with a factor between 0 and 1, this scaling having been derived based on the shear rate in SI units. Therefore, it is converted from the default LBM units to SI units using ˙γSI= (1/∆t) ˙γLBM

3.4 Lisa: the SURFsara cluster computer

The simulations done for the present work are executed on a cluster computer named Lisa, managed by the company SURFsara. A word of thanks to SURFsara for granting me per-mission to make use of Lisa for the present work.

(32)

4 Results

4.1 Test cases with analytical solutions

The accuracy of the advection-diffusion algorithm can be examined by comparing the con-centration field it produces to a corresponding analytical solution for the concon-centration field. Because analytical solutions are only available for simple cases, such comparisons will typi-cally require that simplifications are made with respect to the original, complete model. In the next two paragraphs, the results for pure advection and pure diffusion are both sepa-rately compared to analytical solutions. These two test cases correspond to the diffusion of an instantaneous point source with constant diffusivity and the advection of a concentration distribution with constant velocity, both in one dimension.

4.1.1 Diffusion of an instantaneous point source

The concentration due to an amount of mass M of a substance with diffusivity D released into a volume at t = 0 and at x = 0 with initial condition c(x) = M δ(x) equals: 13

c(x, t) = M

A√4πDte

−x2_/4Dt

(58) Where A is the cross-sectional area of the volume perpendicular to x. Because this solution presumes that all of the mass is released at x = 0 precisely, it cannot be fully adequately com-pared to its numerical counterpart. In the model, the spatial resolution of the concentration field is finite and equal to the spatial resolution of the lattice. Therefore, the concentration release at t = 0 and x = 0 cannot be any narrower than ∆x in width.

For this test case, diffusion along a line of 21 lattice nodes in x is considered. The con-centration is set equal to 10 on the middle lattice node and equal to 9 at all other lattice nodes. This means the analytical solution equals the solution of (58) plus 9. The analytical solution is computed in lattice units, so that M/A = 1. Page 33 shows the comparison of the concentrations computed by the model with the analytical concentrations for both the 2nd order accurate data processor (left hand side) and the fourth order accurate data processor (right hand side) for this test case. For each of the data processors, three figures are included showing the development of the solutions in time. The diffusivity is set equal to 5 · 10−10, based on (45) with an average value of the shear rate probed from the modelled flow field itself. Page 34 shows results the same test case, but with diffusivity equal to 10 · 10−10 and the concentration jump increased from 10 vs. 9 to 10 vs. 8.

(33)

-10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 Distance (lattice nodes) (× 37.5506 μμ)

8.8 8.9 9.0 9.1 9.2 9.3 9.4 9.5 9.6 9.7 9.8 9.9 10.0 10.1 Conc ent rat ion Time (s) = 0.823968

Diffusion: finite difference 2nd order vs. analytical solution

Finite difference 2nd order Analytical solution

Figure 1a: Comparison of the 2nd order accurate numerical solution produced by the model with the analytical solution for diffusion of the initial concentration configuration: c = 10 for x = 0 and c = 9 for x 6= 0. Simulation parameters:

Diffusion: finite difference 4th order vs. analytical solution

Finite difference 4th order Analytical solution

Figure 1b: Comparison of the 4th order accurate numerical solution produced by the model with the analytical solution for diffusion of the initial concentration configuration: c = 10 for x = 0 and c = 9 for x 6= 0. Simulation parameters:

(34)

7.50 7.75 8.00 8.25 8.50 8.75 9.00 9.25 9.50 9.75 10.00 10.25 Conc ent rat ion Time (s) = 0.823968

Figure 2a: Comparison of the 2nd order accurate numerical solution produced by the model with the analytical solution for diffusion of the initial concentration configuration: c =

Figure 2b: Comparison of the 4th order accurate numerical solution produced by the model with the analytical solution for diffusion of the initial concentration configuration: c =

Incorporating concentration fields governed by advection-diffusion into Palabos-based CFD models in three dimensions using finite differences