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.
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-
able non-periodic suspension simulations.
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
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
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
vessel with an aneurysm.
Keywords: suspensions, simulations, boundary conditions, hemodynamics, biological fluid dynamics 1. Introduction
Cell resolved blood flows are studied with various computer models, see e.g. Fedosov et al. 
or Mountrakis et al. . 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  on top of a lattice Boltzmann fluid solver .
their model  to mimic and explain various effects
These boundaries are important since the dis-
represent the red blood cells and/or platelets in
outlet flows) and another solution has to be used .
(e.g. the F˚ ahræus-Lindqvist effect  or thrombus
tribution of cells in a blood vessel is non-trivial.
formations ). 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
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. . The in- and
Periodic boundary Inflow region
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-
justed for IB-LBM blood suspension simulations.
Furthermore Lykov’s approach is only validated for
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. . 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
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
main and the flow domain. The restriction is en-
forced by only copying information from the preinlet domain to the flow domain. This means that the fluid field and immersed boundary particles
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
is done in the inflow region.
al.  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 . 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
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
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
coupled from the pre-inlet).
fluid boundary problem that is well understood for
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-
domain to solve the distribution function for this
He boundary at the beginning of the inflow region,
. The immersed boundary method only couples
all the ghost lsps become normal lsps in the in-
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-
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
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
Pre-inlet domain A
Inflow region Flow domain
son the removal is done in the outflow region.
the experiment and chosen boundary condition. It
The size of the outflow region is dependent on
should be at least as long as the length of the longest
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.
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
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. . The magnitude
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
constant this results in a lower mean velocity and lower recovered Reynolds number. We measured
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
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.
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
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.
the red blood cells to relax from their initial disk shape to their biconcave shape. A quadratic in-
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
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
immersed boundary method . 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. .
Mountrakis et al.  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
et al. . 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
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
this paper. In all the simulations of the pre-inlet
domains the cell-free layer is 5 µm according to this
definition (see Figure 4). A platelet is considered
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
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.  and Zhao
only the pre-inlet is filled with particles. The flow
To test the in- and outlet boundary conditions
domain only consists of fluid. Figure 6 shows the
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
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
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
pre-inlet with no visible artifacts. At the end of the
outflow region the Zou-He density outlet boundary
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
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
The transformation from the parabolic velocity
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
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
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-
To show the necessity of our in- and outlet con-
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
0.1 cm and an attached aneurysm. We let blood flow through it for 0.4 seconds with a velocity of
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
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
our new in- and outlet conditions does not suffer
with periodic boundaries. This is due to deple-
from these effects.
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
200 µm and a length of 500 µm.
4. Discussion and Conclusion We have shown that Lykov’s et al.  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-
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-
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.  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-
case it is possible to do particle residence time anal-
observed in earlier work by Mountrakis et al. in
ysis  and other analyses that were previously
only available to pure fluid solvers within hemodynamics. In future work we intend to investigate the behaviour of platelets in aneurysms in more de-
tail, also taking into account the pulsatility of blood flow.
Pre-Inlet domain 600 µm
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
under grant agreement no.
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
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.
The material model used by Mountrakis et al.  suffered from instabilities in our simulations.
equations are adjusted.
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
In the simulations the red blood cells appeared to be too stiff. This is solved by squaring the anReferences
gle used in calculating the bending resistance force (ftrs ). The bending constant Ctrs is decreased a
 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
We noticed numerical instability caused by large
42 (2) (2014) 368–387.
forces from the repulsion force (Frep ). This seemed
 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).  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.
 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.  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
 M. Tomaiuolo, T. J. Stalker, J. D. Welsh, S. L. Di-
 K. Lykov, X. Li, H. Lei, I. V. Pivkin, G. E. Karniadakis,
Inflow/outflow boundary conditions for particle-based
 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.
 E. J. Carboni, B. H. Bognet, G. M. Bouchillon, A. L.
 C. S. Peskin, The immersed boundary method, Acta
Kadilak, L. M. Shor, M. D. Ward, A. W. Ma, Direct
numerica 11 (2002) 479–517.
tracking of particles and quantification of margination
 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.
 Z. Pleunis, E. Lorenz, L. Mountrakis, R. Sprik, T. van
 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-
blood cell suspension, Personal Communication.
teriolar bifurcation, Microvascular research 106 (2016) 14–23.
 L. Mountrakis, E. Lorenz, A. G. Hoekstra, Validation
 R. Zhao, M. V. Kameneva, J. F. Antaki, Investigation of platelet margination phenomena at elevated shear
stress, Biorheology 44 (3) (2007) 161–177.  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,
of an efficient two-dimensional model for dense suspen-
 L. Mountrakis, E. Lorenz, A. Hoekstra, Scaling of shear-
Springer, 2016, pp. 263–277.  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
 J. Latt, B. Chopard, O. Malaspinas, M. Deville,
flow through aneurysmal vessels, Interface focus 3 (2)
induced diffusion and clustering in a blood-like suspen-
A. Michler, Straight velocity boundaries in the lattice
boltzmann method, Physical Review E 77 (5) (2008) 056703.
 S. Succi, The lattice Boltzmann equation: for fluid dynamics and beyond, Oxford university press, 2001.
 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.  M. Brust, C. Schaefer, R. Doerr, L. Pan, M. Gar-