Inflow and outflow boundary conditions for 2D suspension simulations with the immersed boundary lattice Boltzmann method

Inflow and outflow boundary conditions for 2D suspension simulations with the immersed boundary lattice Boltzmann method

Accepted Manuscript Inflow and outflow boundary conditions for 2D suspension simulations with the immersed boundary lattice Boltzmann method Victor W...

2MB Sizes 0 Downloads 48 Views

Accepted Manuscript

Inflow and outflow boundary conditions for 2D suspension simulations with the immersed boundary lattice Boltzmann method Victor W. Azizi Tarksalooyeh, Gabor Zavodszky, Britt J. M. van Rooij, ´ ´ Alfons G. Hoekstra PII: DOI: Reference:

S0045-7930(18)30221-4 10.1016/j.compfluid.2018.04.025 CAF 3870

To appear in:

Computers and Fluids

Received date: Revised date: Accepted date:

19 October 2017 9 March 2018 17 April 2018

Please cite this article as: Victor W. Azizi Tarksalooyeh, Gabor Zavodszky, Britt J. M. van Rooij, ´ ´ Alfons G. Hoekstra, Inflow and outflow boundary conditions for 2D suspension simulations with the immersed boundary lattice Boltzmann method, Computers and Fluids (2018), doi: 10.1016/j.compfluid.2018.04.025

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • New in- and outflow boundary methods for suspension simulations are proposed. • The in- and outflow boundary methods provide

CR IP T

an undisturbed inflow of particles. • The in- and outflow boundary conditions en-

AC

CE

PT

ED

M

AN US

able non-periodic suspension simulations.

1

ACCEPTED MANUSCRIPT

Inflow and outflow boundary conditions for 2D suspension simulations with the immersed boundary lattice Boltzmann method Victor W. Azizi Tarksalooyeha , G´abor Z´avodszkya , Britt J. M. van Rooija , Alfons G. Hoekstraa,b Science Lab, Institute of Informatics, University of Amsterdam, The Netherlands b ITMO University, Saint-Petersburg, Russian Federation

CR IP T

a Computational

Abstract

In- and outflow boundary conditions for 2D immersed boundary lattice Boltzmann suspension simulations, applied to cell based blood flow models, are presented. The inlet is constructed with an one-way coupling

AN US

to a periodic domain containing a correct distribution of suspended particles. This provides an inflow of particles that has a correct distribution and is decoupled from any phenomena in the flow domain. An outflow boundary for the particles that does not influence the distribution of particles in the flow domain is also constructed. With this a method to run long (> 1s) cell based blood flow simulations within any type of domain is provided. These boundary conditions are then used for a simulation of blood flow in a curved

M

vessel with an aneurysm.

ED

Keywords: suspensions, simulations, boundary conditions, hemodynamics, biological fluid dynamics 1. Introduction

PT

Cell resolved blood flows are studied with various computer models, see e.g. Fedosov et al. [1]

CE

or Mountrakis et al. [2]. These models explicitly

In this paper inflow and outflow flow boundaries for complex flows are constructed for suspensions that are modeled with an immersed boundary method [7] on top of a lattice Boltzmann fluid solver [8].

their model [3] to mimic and explain various effects

These boundaries are important since the dis-

AC

represent the red blood cells and/or platelets in

outlet flows) and another solution has to be used [6].

(e.g. the F˚ ahræus-Lindqvist effect [4] or thrombus

tribution of cells in a blood vessel is non-trivial.

formations [5]). In these studies mostly periodic

This is because of phenomena such as the cell-

boundary conditions (PBC) are used for the flow

free layer, F˚ ahræus-Lindqvist effect or margination

domain. Sometimes, however, it is not possible to

of platelets. While periodic boundary conditions

use periodic boundaries (e.g single inlet, multiple

might circumvent this problem, they themselves can impose further complications. In a system with

Email address: [email protected] (Victor W. Azizi Tarksalooyeh) Preprint submitted to Elsevier

PBC the outlet can influence the distribution of April 20, 2018

ACCEPTED MANUSCRIPT

particles in the inlet boundary condition (eg. imag-

In a computational context cells are referred to

ine a thrombus where all platelets get stuck, then

as particles, while the Lagrangian structure points

this directly causes a depletion of active platelets).

that make up the particles are referred to as lsps.

Therefore, we propose another type of inlet and 2.2. In- and outflow boundary conditions

outlet inspired by Lykov et al. [6]. The in- and

Periodic boundary Inflow region

Outlet

CR IP T

outlets are decoupled (the outlet cannot influence the inlet) and a separate tubular domain provides a correct distribution of particles and fluid for the inlet boundary condition. Lykov’s approach is ad-

Pre-inlet domain

justed for IB-LBM blood suspension simulations.

Flow domain

AN US

Furthermore Lykov’s approach is only validated for

Inlet

Outflow region

Figure 1: The various elements necessary for defining the inand outlet boundary conditions.

very narrow domains (< 10µm). Our method is also validated for wider domains (∼ 200 µm). A similar

A tubular domain with periodic boundary condi-

approach is proposed by Ye et al. [9]. We extend

tions for both the particles and solvent (pre-inlet)

their method by providing a detailed explanation

is placed in front of the inlet (see Figure 1). This

M

on how to construct such in- and outlet boundary conditions and how key parameters can be chosen.

Finally, we also validate the in- and outlet boundary

ED

main and the flow domain. The restriction is en-

PT

forced by only copying information from the preinlet domain to the flow domain. This means that the fluid field and immersed boundary particles

CE

2. Methods

tem in any way. This property implies a one-way communication restriction between the pre-inlet do-

conditions and show their application in a selected use case.

domain cannot be influenced by the rest of the sys-

2.1. Material model for the particles

must be reconstructed in the flow domain. This

The material model developed by Mountrakis et

AC

is done in the inflow region.

al. [10] is used for simulating red blood cells and

At the end of the flow domain particles have to

platelets. This material model is also used for shear

be removed without influencing the flow domain.

experiments [11]. We have tuned this model to in-

This is not possible with a one-way coupling as with

crease numerical stability (see appendix A). The

the pre-inlet domain and flow domain. Therefore

cells are initialized as disks that relax towards their

an outflow region is defined in which particles are

final shape (biconcave for RBCs, platelets stay el-

deleted. The resulting disturbances have no influ-

lipse shaped) in the first millisecond of a simulation.

ence on the flow domain before the outlet boundary. 3

ACCEPTED MANUSCRIPT

The flow domain is encompassed by the inlet and

the pre-inlet domain into the inflow region the lsp

outlet. Between these boundaries the distribution

in the inflow region cannot interact with the lsp

of particles and fluid is well defined. In the in- and

in the pre-inlet domain. Therefore, this can cause

outflow regions the distribution of particles and the

numerical instability of the particle in the inflow

fluid flow might not be correct.

region. Because of this a more complex system for the inlet for particles is proposed.

CR IP T

It is possible to use different time or space

When a single lsp enters the inflow region it first

the one-way communication restriction between the

becomes a ghost lsp (see figure 2). A ghost lsp is

pre-inlet domain and the flow domain. It is even

used to calculate repulsion for other lsps in the in-

possible to use an entirely different model, as long

flow region. This is to ensure no overlap occurs.

as the results of that can be interpolated to the

A ghost has no other interaction with lsps and it

AN US

discretisation in the pre-inlet domain because of

inflow region of the flow domain.

follows the location of the lsp in the pre-inlet domain it was copied from. To keep the velocity field

2.2.1. Inflow region

of the fluid from deviating from the one in the pre-

At the beginning of the inflow region there is

inlet domain, which would give rise to discontinu-

a fluid row which has no neighbours (as it is de-

This is a classical

M

coupled from the pre-inlet).

fluid boundary problem that is well understood for

ED

lattice Boltzmann methods [12, 13]. We use the

ities when the particle is finally ”dropped into” the inflow region, the force from a ghost lsp is communicated through the immersed boundary method to the fluid field in the inflow region as well.

macroscopic velocity of a row from the pre-inlet

When the full ghost particle has crossed the Zou-

PT

domain to solve the distribution function for this

He boundary at the beginning of the inflow region,

[14]. The immersed boundary method only couples

all the ghost lsps become normal lsps in the in-

CE

boundary with the method proposed by Zou and He

flow region. This happens immediately after the

thus the particles are not affected by the disconti-

last ghost lsp of the ghost particle enters the inflow

nuities in the distribution functions.

region. This solves the problem of numerical insta-

AC

through the macroscopic fluid force and velocity,

The inflow of particles is more intricate. When

bility of a particle while transferring into the flow

a lsp of a particle crosses the row from which the

domain. The transformation of all ghost lsps into

macroscopic velocity is copied for the fluid, it is

normal lsps can be seen as “dropping a particle into

copied to the inflow region. However, all lsps of a

the fluid”.

particle have to interact to keep the particle in its

The size of the inflow region is dependent on the

correct shape. When a particle is transferring from

experiment, but should be at least as long as the 4

ACCEPTED MANUSCRIPT

Pre-inlet domain A

Inflow region Flow domain

son the removal is done in the outflow region.

A

the experiment and chosen boundary condition. It

The size of the outflow region is dependent on

B

should be at least as long as the length of the longest

Inlet

possible elongated particle. When the boundary

Figure 2: Particle A is transferring into the inflow region. All the lsps of particle A in the inflow region are ghosts. Particle B has no more ghost particles because it fully entered the inflow region and therefore has been dropped into the fluid. Particle B cannot overlap with particle A because it is repelled from the ghost lsps of particle A.

CR IP T

condition introduces artifacts in the outflow region the size of the outflow region should be at least as large to encompass all these artifacts.

length of the longest possible elongated particle in

2.3. Macroscopic fluid model

the experiment. With blood flows this is the red

The fluid is simulated using the lattice Boltz-

blood cell with a maximum length of about 15 µm.

AN US

mann method with the BGK collision operator.

2.2.2. Outflow region

The lattice spacing is ∆x = 1 µm. The fluid that is simulated is blood plasma at 37 ◦C, therefore, the

a boundary condition. The same Zou-He boundary

kinematic viscosity is 1.28 × 10−6 m2 /s [15, 16] and

as in the inflow region is used, but density is pre-

the density is 1025 kg/m3 . The timestep is set to

scribed instead of velocity to calculate the unknown

∆t = 1.176 × 10−7 s. The relaxation time (τ ) of the

distributions. The density that is used for this is

fluid is 1.1. With the particles present this gives an

the initial density of the flow domain ρ0 , therefore

observed locally maximum compressibility error of

ρ0 = ρout . This boundary condition does not en-

1% (as observed from our simulations). This means

force a parabolic profile as expected in a tubular

that the quasi-incompressibility constraint of the

PT

ED

M

At the end of the outflow region there must be

domain. Therefore, a flattening of the velocity pro-

lattice Boltzmann method is conserved.

file in the outflow region is to be expected. A particle is removed fully when one of its lsps

an external forcing term, which is implemented

crosses the last row of the outflow region (and con-

as proposed by Guo et al. [17]. The magnitude

AC

CE

In the pre-inlet domain the fluid is driven by

sequently behaviour of the particle becomes unde-

of this forcing term is calculated from the Hagen-

fined). The pressure gradient from the fluid in the

Poiseuille equation adjusted for flow between par-

particle as opposed to the fluid outside it is approx-

allel plates and a given Reynolds number. In our

imately zero. Because of this there are no shock-

experiments we choose Replasma = 40. This results

waves resulting of the deletion of a particle. How-

in v¯plasma = 17 cm/s for a domain with a width

ever, the deletion of a particle does influence the

of 200 µm. In our experiments the cells add ex-

particle distribution and flow profile. For this rea-

tra viscosity to the fluid. Since our driving force is 5

ACCEPTED MANUSCRIPT

3. Results

constant this results in a lower mean velocity and lower recovered Reynolds number. We measured

3.1. Pre-Inlet

the resulting mean velocity of the suspension which Firstly, five differently sized pre-inlet domains are

is v¯susp = 7.5 cm/s . Since our flow exhibits lami-

tested to ensure that a correct distribution of par-

nair behaviour we assume that

Resusp ≈

CR IP T

ticles enters the system. To this end five simuv¯susp Replasma v¯plasma

which gives us Resusp ≈ 18.

lations of each of the five different geometries for the pre-inlet domain are performed.

The width

of the pre-inlet domain is 200µm for all geome-

Reynolds number varies slightly with the local

tries (Ld = 200µm). The geometries only differ

hematocrit of the system since the presence of cells

in length in the direction of the periodic boundary

AN US

Furthermore this

in the fluid changes the viscosity.

condition (Lp ), giving them different aspect ratios

In the beginning (t = 0 s) of a simulation, the

defined as the length of a channel divided by its

body force is increased quadratically from zero to

width (Lp /Ld ). The aspect ratios range from 0.5

the final magnitude over one millisecond to allow

to 2.5 . The rest of the parameters are the same.

M

the red blood cells to relax from their initial disk shape to their biconcave shape. A quadratic in-

ED

crease is chosen for numerical stability.

2.3.1. Immersed boundary implementation

parameter Diameter, Ld Length, Lp Hematocrit, H Platelet Ratio, P Reynolds Number, Re Simulation Time, s

value 200 µm 100 − 500 µm 42 % 1 20

18 1.18 s

PT

Table 1: Dimensions and constants for the pre-inlet plane Poiseuille flow

The fluid-particle coupling is modeled through an All pre-inlets retrieve similar velocity profiles.

boundary method the velocity of the fluid is in-

Therefore, only the pre-inlet domain with aspect

terpolated to the particles at every timestep. The

ratio 2.5 is plotted in Figure 3. Due to the shear

AC

CE

immersed boundary method [7]. In this immersed

same is true for the force exercised by the material

thinning behaviour of blood the velocity profile is

model of the particles on the fluid. This interpola-

more plug-like than the parabolic profile for a fluid

tion is done to the closest Eulerian points X with a

without particle that has the same average velocity.

discrete Dirac delta function as constructed by by

These finding are supported by Carboni et al. [18].

Mountrakis et al. [10] which effectively gives a fourpoint stencil with weights to the closest points. No

In flowing whole blood the platelets will

repulsion particles are used at the fluid boundaries.

marginate to the cell-free layer. 6

This is used

ACCEPTED MANUSCRIPT

et al. [20]. All pre-inlets achieve the same average margination. All pre-inlets seem to recover a correct distribution, however, the system with the largest aspect ratio inherently has the smallest fiFigure 3: The flow velocity profile of the pre-inlet with aspect ratio 2.5 perpendicular to the flow direction. The black dotted line is the parabolic profile for a fluid flow without particles present in a plane Poiseuille flow with the same average velocity.

nite size system effects. Therefore we have used the

CR IP T

pre-inlet with the largest aspect ratio (2.5) in the following experiment.

as a measure of correctness for the pre-inlets. The margination is measured as the percentage of platelets in the cell-free layer. In relevant literature

AN US

there is no clear definition of the size of the cell-free layer. Because of this we define it as the maximum

length from the wall where there is a hematocrit

Figure 5: A plot of the percentage of marginated platelets for the different pre-inlet domains over time.

of at most 50% of the mean hematocrit. This definition for the cell-free layer is used in the rest of

3.2. Plane Poiseuille Flow

M

this paper. In all the simulations of the pre-inlet

domains the cell-free layer is 5 µm according to this

ED

definition (see Figure 4). A platelet is considered

CE

PT

marginated when it is present within this layer.

Figure 4: A plot of the hematocrit (solid red line) and platelet concentration (dotted black line) perpendicular to the flow. Fifty samples from five experiments from the largest pre-inlet (500 µm) are used. The samples are taken after a simulation time of 1.0 s. The location of the cell-free layer is plotted on the right.

parameter Height, d System Length, L Pre-inlet Length, Lp Hematocrit, H Platelet Ratio, P Reynolds Number, Re Sample Step Size Simulation Time Inflow region Outflow region

value 200 µm 500 µm 500 µm 42% 1 20

18 0.00118 s 1.18 s 50 µm 50 µm

AC

Table 2: Dimensions and constants for the plane Poiseuille flow experiment

Figure 5 shows the percentage of marginated

a plane Poiseuille flow is simulated with a pre-

platelets over time. The initial margination hap-

inlet domain of the same length as the flow do-

pens quickly and then asymptotically tends to the

main (500 µm). See table 2 for the parameters that

final margination percentage, this finding is sup-

are used. At the start (t = 0 s) of the simulation

ported by the work of Pleunis et al. [19] and Zhao

only the pre-inlet is filled with particles. The flow

To test the in- and outlet boundary conditions

7

ACCEPTED MANUSCRIPT

domain only consists of fluid. Figure 6 shows the

Pre-inlet domain

domain being filled up after turning on the body force in the pre-inlet domain. After 0.01s the flow domain equilibrates around the same hematocrit as

Inflow region

Flow domain

Outflow region

CR IP T

the pre-inlet domain.

Figure 7: Velocity magnitude of the fluid in the pre-inlet and the system after 1 second, particles are omitted for visibility.

main with particles present. The red particles are red blood cells, while the black dots represent the

AN US

Figure 6: A plot of the hematocrit over time for the first 0.07 seconds. The hematocrit is measured and averaged over the full domain (including in- and outflow regions in the flow domain).

platelets. The graphs on the bottom of Figure 8 show the averaged hematocrit and platelet concen-

In Figure 7 the velocity magnitude of the fluid of

tration for various locations in the flow domain and

the pre-inlet domain coupled with the flow domain

pre-inlet domain. Through the transition of the

is shown. The Zou-He velocity inlet attaches to the

pre-inlet domain to the flow domain no noticeable

M

pre-inlet with no visible artifacts. At the end of the

outflow region the Zou-He density outlet boundary

ED

treats all unknown lattice Boltzmann nodes as fluid

change in this profile is visible. At the end of the outflow region a slight disturbance is visible due to deletion of particles and Zou-He boundary condi-

nodes. Therefore the outlet boundary does not act

tion artifacts.

PT

as an infinite elongation of the plane Poiseuille flow. This results in an more uniform velocity magnitude Inflow Outflow region Pre-Inlet domain Flow domain region

pected.

Platelet concentration

The transformation from the parabolic velocity

AC

magnitude profile towards a more uniform profile of the fluid is present in the last 50 µm of the domain. Deviations from the correct particle distributions are not present anymore at that distance from the outlet. Therefore we assume 50 µm to be a large

Location (µm)

Hematocrit (%)

CE

profile towards the end of the outflow region as ex-

Figure 8: A comparison of the hematocrit and platelet concentration at different places in the system. The graphs are averaged over 10 samples taken after 1 second of simulation.

enough outflow region for this simulation. Figure 8 shows the pre-inlet domain and flow do8

ACCEPTED MANUSCRIPT

Flow direction

CR IP T

Figure 9: The profiles from figure 8 shown again in a single picture. The solid red line is the hematocrit of the pre-inlet. The solid green lines are the hematocrit of the other three places from figure 8. The dotted red and green lines are the standard deviations for these profiles respectively.

3.3. Same geometry with pre-inlet and periodic con-

(a) Periodic

ditions

AN US

To show the necessity of our in- and outlet con-

(b) Pre-inlet

ditions we conduct two experiments with a single geometry. The geometry is such that it can be ei-

ther run with periodic boundaries or a pre-inlet. The geometry is a straight pipe with diameter of

(c)

M

0.1 cm and an attached aneurysm. We let blood flow through it for 0.4 seconds with a velocity of

ED

20 cm/s.

In figure 10 we can see clear differences arising

unwanted and noticeable when running long simu-

between the same geometry with different bound-

lations of such a geometry with periodic boundary

PT

ary conditions.

Figure 10: The same geometry run with periodic (a) and pre-inlet (b) boundaries after 0.4 s. The red cells represent red blood cells and the black cells represent platelets. In (c) the number of platelets in the aneurysm is plotted over time.

The aneurysm with pre-inlet is

conditions. The geometry that is simulated with

enriched with more platelets than the aneurysm

CE

our new in- and outlet conditions does not suffer

with periodic boundaries. This is due to deple-

from these effects.

AC

tion of platelets in the tube with the periodic condi3.4. Application to a curved aneurysm

platelets than are available in the whole periodic

As an example of applying the in- and outlet

system after 0.23 s. Furthermore, the density of

conditions for an application where periodic bound-

particles in the tube with periodic boundary con-

ary conditions cannot be used we simulate the flow

ditions is lower in the aneurysm. We account this

of blood within a curved vessel that has an at-

effect due to the non-changing number of particles

tached aneurysm. Again, the pre-inlet domain that

in the whole periodic geometry. Both effects are

is attached to the curved aneurysm has a width of

tions. The aneurysm with pre-inlet even has more

9

ACCEPTED MANUSCRIPT

200 µm and a length of 500 µm.

4. Discussion and Conclusion We have shown that Lykov’s et al. [6] method

the same hematocrit and platelet percentage as

can be adapted for two dimensional suspensions

the pre-inlet domain. Otherwise the flow domain

simulated with IB-LBM. Furthermore, we show

(aneurysm) would take in the order of seconds

that these two dimensional suspensions can scale

to reach the correct hematocrit.

In the plane

to larger inflow boundaries of 200µm. The tests

Poiseuille flow experiment all these particles would

with the plane Poiseuille flow show that an accu-

be replaced with particles from the pre-inlet domain

rate distribution of particles within the flow do-

within 0.1 s. Thus giving approximately the same

main and a well behaving blood flow is retrieved

simulation results after this time.

with our in- and outlet boundaries. The simula-

AN US

CR IP T

The domain is also randomly initialized with

tion with the curved vessel and attached aneurysm

simulation. In this figure the location of a platelet is

demonstrates the importance of such in- and outlet

represented by a black dot. The red blood cells are

boundary conditions. Since platelets are accumu-

omitted for clarity. On the background the velocity

lating in the aneurysm it is necessary to have an

of the fluid and the streamlines are plotted. We

inflow that is independent of the behaviour of par-

M

Figure 11 shows the aneurysm after one second of

can see that platelets start to accumulate in the

aneurysm as predicted by Z´avodszky et al. [21] and

ticles in the system. These in- and outflow boundary conditions can be extended into three dimensions as a logical next

a straight vessel geometry with periodic boundary

step. As an extension to the curved aneurysm show-

conditions [22].

case it is possible to do particle residence time anal-

PT

ED

observed in earlier work by Mountrakis et al. in

ysis [21] and other analyses that were previously

CE

only available to pure fluid solvers within hemodynamics. In future work we intend to investigate the behaviour of platelets in aneurysms in more de-

AC

500 µm

200 µm

tail, also taking into account the pulsatility of blood flow.

Pre-Inlet domain 600 µm

Acknowledgements

Figure 11: Curved Aneurysm with platelets (black dots), velocity magnitude of the fluid (color bar) and streamlines (white lines). The fluid is moving from the left to the top

This work was supported by the European Union Horizon 2020 research and innovation programme 10

ACCEPTED MANUSCRIPT

under grant agreement no.

merically unstable.

675451, the Comp-

BioMed project and grant agreement no. 671564,

As a final step some parameters for the material

the ComPat project.

model are slightly adjusted (see table 3).

This work was sponsored by NWO Exacte

Parameter Spring constant, Cspr Cell-cell constant, Crep Cutoff distance, hcutoff

Wetenschappen (Physical Sciences) for the use

CR IP T

of supercomputer facilities, with financial support

Value 1 × 10−3 N/m 2 × 10−22 N · m2 0.6 µm

from the Nederlandse Organisatie voor Weten-

Table 3: Adjusted constants

schappelijk Onderzoek (Netherlands Organization

In the simulations it became apparent that the re-

for Science Research, NWO).

pulsion force causes the particles to keep a distance of hcutoff from each other at all times. Because of

AN US

Appendices

this the particles act as if they have a diameter that is 0.5 × hcutoff = 0.3µm larger. Therefore this extra

A. Variables and Formulae for the Material model

area is added to the calculation for the hematocrit. This is shown as the apparent area in Table 4.

M

The material model used by Mountrakis et al. [10] suffered from instabilities in our simulations.

equations are adjusted.

ED

Therefore, some constants and material model

Parameter Area Apparent area Circumference Nlsp

Red blood cells 10.00 µm2 15.30 µm2 16.71 µm 26

Platelets 1.81 µm2 5.03 µm 8

Table 4: Cell properties

PT

In the simulations the red blood cells appeared to be too stiff. This is solved by squaring the anReferences

CE

gle used in calculating the bending resistance force (ftrs ). The bending constant Ctrs is decreased a

[1] D. A. Fedosov, M. Dao, G. E. Karniadakis, S. Suresh,

hundredfold to 1 × 10−9 N/rad.

Computational biorheology of human blood flow in health and disease, Annals of biomedical engineering

AC

We noticed numerical instability caused by large

42 (2) (2014) 368–387.

forces from the repulsion force (Frep ). This seemed

[2] L. Mountrakis, Transport of blood cells studied with

to be happening because of a exponential term for

fully resolved models, Ph.D. thesis, Faculty of Science,

the distance. Therefore this term is removed and

University of Amsterdam (2015). [3] E. Moeendarbary, T. Ng, M. Zangeneh, Dissipative par-

made linear. This dampens the force impact of col-

ticle dynamics: introduction, methodology and complex

lisions that with the quadratic term could create

fluid applicationsa review, International Journal of Ap-

such large forces that the simulation became nu-

plied Mechanics 1 (04) (2009) 737–763.

11

ACCEPTED MANUSCRIPT

[4] R. F˚ ahræus, T. Lindqvist, The viscosity of the

cia, P. Arratia, C. Wagner, Rheology of human blood

blood in narrow capillary tubes, American Journal of

plasma: Viscoelastic versus newtonian behavior, Phys-

Physiology–Legacy Content 96 (3) (1931) 562–568.

ical review letters 110 (7) (2013) 078305. [16] S. A. Peters, M. Woodward, A. Rumley, H. D.

amond, T. Sinno, L. F. Brass, A systems approach

Tunstall-Pedoe, G. D. Lowe, Plasma and blood vis-

to hemostasis: 2. computational analysis of molecu-

cosity in the prediction of cardiovascular disease and

lar transport in the thrombus microenvironment, Blood

mortality in the scottish heart health extended cohort

124 (11) (2014) 1816–1823.

study, European journal of preventive cardiology (2016)

CR IP T

[5] M. Tomaiuolo, T. J. Stalker, J. D. Welsh, S. L. Di-

[6] K. Lykov, X. Li, H. Lei, I. V. Pivkin, G. E. Karniadakis,

2047487316672004.

Inflow/outflow boundary conditions for particle-based

[17] Z. Guo, C. Zheng, B. Shi, Discrete lattice effects on the

blood flow simulations: Application to arterial bifur-

forcing term in the lattice boltzmann method, Physical

cations and trees, PLoS Comput Biol 11 (8) (2015)

Review E 65 (4) (2002) 046308.

e1004410.

[18] E. J. Carboni, B. H. Bognet, G. M. Bouchillon, A. L.

[7] C. S. Peskin, The immersed boundary method, Acta

AN US

Kadilak, L. M. Shor, M. D. Ward, A. W. Ma, Direct

numerica 11 (2002) 479–517.

tracking of particles and quantification of margination

[8] S. Chen, G. D. Doolen, Lattice boltzmann method for

in blood flow, Biophysical journal 111 (7) (2016) 1487–

fluid flows, Annual review of fluid mechanics 30 (1) (1998) 329–364.

1495.

[19] Z. Pleunis, E. Lorenz, L. Mountrakis, R. Sprik, T. van

[9] S. S. Ye, M. Ju, S. Kim, Recovery of cell-free layer and

Leeuwen, Shear-induced transport of platelets in a red

wall shear stress profile symmetry downstream of an ar-

M

blood cell suspension, Personal Communication.

teriolar bifurcation, Microvascular research 106 (2016) 14–23.

[10] L. Mountrakis, E. Lorenz, A. G. Hoekstra, Validation

[20] R. Zhao, M. V. Kameneva, J. F. Antaki, Investigation of platelet margination phenomena at elevated shear

ED

stress, Biorheology 44 (3) (2007) 161–177. [21] G. Z´ avodszky, G. K´ arolyi, I. Szikora, G. Pa´ al, Frac-

sions of red blood cells, International Journal of Modern

tals and chaos in the hemodynamics of intracranial

Physics C 25 (12) (2014) 1441005.

aneurysms, in: The Fractal Geometry of the Brain,

PT

of an efficient two-dimensional model for dense suspen-

[11] L. Mountrakis, E. Lorenz, A. Hoekstra, Scaling of shear-

Springer, 2016, pp. 263–277. [22] L. Mountrakis, E. Lorenz, A. Hoekstra, Where do the

sion, EPL (Europhysics Letters) 114 (1) (2016) 14002.

platelets go? a simulation study of fully resolved blood

[12] J. Latt, B. Chopard, O. Malaspinas, M. Deville,

flow through aneurysmal vessels, Interface focus 3 (2)

CE

induced diffusion and clustering in a blood-like suspen-

A. Michler, Straight velocity boundaries in the lattice

(2013) 20120089.

AC

boltzmann method, Physical Review E 77 (5) (2008) 056703.

[13] S. Succi, The lattice Boltzmann equation: for fluid dynamics and beyond, Oxford university press, 2001.

[14] Q. Zou, X. He, On pressure and velocity boundary conditions for the lattice boltzmann bgk model, Physics of Fluids (1994-present) 9 (6) (1997) 1591–1598. [15] M. Brust, C. Schaefer, R. Doerr, L. Pan, M. Gar-

12