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

Inand 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 ( > 1 s) 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. © 2018 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license. ( http://creativecommons.org/licenses/by/4.0/ ) b c a c s f f Y p d v p


Introduction
Cell resolved blood flows are studied with various computer models, see e.g. Fedosov et al. [1] or Mountrakis et al. [2] . These models explicitly represent the red blood cells and/or platelets in their model [3] to mimic and explain various effects (e.g. the Fåhraeus-Lindqvist effect [4] or thrombus formations [5] ). In these studies mostly periodic boundary conditions (PBC) are used for the flow domain. Sometimes, however, it is not possible to use periodic boundaries (e.g. single inlet, multiple outlet flows) and another solution has to be used [6] . 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] .
These boundaries are important since the distribution of cells in a blood vessel is non-trivial. This is because of phenomena such as the cell-free layer, Fåhraeus-Lindqvist effect or margination of platelets. While periodic boundary conditions might circumvent this problem, they themselves can impose further complications. In a system with PBC the outlet can influence the distribution of particles in the inlet boundary condition (e.g. imagine a thrombus where all platelets get stuck, then this directly causes a depletion of active platelets). * Corresponding author.
E-mail address: v.w.azizitarksalooyeh@uva.nl (V.W. Azizi Tarksalooyeh). Therefore, we propose another type of inlet and outlet inspired by Lykov et al. [6] . The in-and 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 adjusted for IB-LBM blood suspension simulations. Furthermore Lykov's approach is only validated for very narrow domains ( < 10 μm). Our method is also validated for wider domains ( ∼ 200 μm). A similar approach is proposed by Ye et al. [9] . We extend their method by providing a detailed explanation 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 conditions and show their application in a selected use case.

Material model for the particles
The material model developed by Mountrakis et al. [10] is used for simulating red blood cells and platelets. This material model is also used for shear experiments [11] . We have tuned this model to increase numerical stability (see Appendix A ). The cells are initialized as disks that relax towards their final shape (biconcave for RBCs, platelets stay ellipse shaped) in the first millisecond of a simulation. In a computational context cells are referred to as particles, while the Lagrangian structure points that make up the particles are referred to as lsps.

In-and outflow boundary conditions
A tubular domain with periodic boundary conditions for both the particles and solvent (pre-inlet) is placed in front of the inlet (see Fig. 1 ). This domain cannot be influenced by the rest of the system in any way. This property implies a one-way communication restriction between the pre-inlet domain and the flow domain. The restriction is enforced by only copying information from the pre-inlet domain to the flow domain. This means that the fluid field and immersed boundary particles must be reconstructed in the flow domain. This is done in the inflow region.
At the end of the flow domain particles have to be removed without influencing the flow domain. This is not possible with a one-way coupling as with the pre-inlet domain and flow domain. Therefore an outflow region is defined in which particles are deleted. The resulting disturbances have no influence on the flow domain before the outlet boundary.
The flow domain is encompassed by the inlet and outlet. Between these boundaries the distribution of particles and fluid is well defined. In the in-and outflow regions the distribution of particles and the fluid flow might not be correct.
It is possible to use different time or space discretisation in the pre-inlet domain because of the one-way communication restriction between the pre-inlet domain and the flow domain. It is even possible to use an entirely different model, as long as the results of that can be interpolated to the inflow region of the flow domain.

Inflow region
At the beginning of the inflow region there is a fluid row which has no neighbors (as it is decoupled from the pre-inlet). This is a classical fluid boundary problem that is well understood for lattice Boltzmann methods [12,13] . We use the macroscopic velocity of a row from the pre-inlet domain to solve the distribution function for this boundary with the method proposed by Zou and He [14] . The immersed boundary method only couples through the macroscopic fluid force and velocity, thus the particles are not affected by the discontinuities in the distribution functions.
The inflow of particles is more intricate. When a lsp of a particle crosses the row from which the macroscopic velocity is copied for the fluid, it is copied to the inflow region. However, all lsps of a particle have to interact to keep the particle in its correct shape. When a particle is transferring from the pre-inlet domain into the inflow region the lsp in the inflow region cannot interact with the lsp in the pre-inlet domain. Therefore, this can cause numerical instability of the particle in the inflow region. Because of this a more complex system for the inlet for particles is proposed.
When a single lsp enters the inflow region it first becomes a ghost lsp (see Fig. 2 ). A ghost lsp is used to calculate repulsion for other lsps in the inflow region. This is to ensure no overlap occurs. A ghost has no other interaction with lsps and it follows the location of the lsp in the pre-inlet domain it was copied from. To keep Fig. 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.
the velocity field of the fluid from deviating from the one in the pre-inlet domain, which would give rise to discontinuities 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.
When the full ghost particle has crossed the Zou-He boundary at the beginning of the inflow region, all the ghost lsps become normal lsps in the inflow region. This happens immediately after the last ghost lsp of the ghost particle enters the inflow region. This solves the problem of numerical instability of a particle while transferring into the flow domain. The transformation of all ghost lsps into normal lsps can be seen as "dropping a particle into the fluid".
The size of the inflow region is dependent on the experiment, but should be at least as long as the length of the longest possible elongated particle in the experiment. With blood flows this is the red blood cell with a maximum length of about 15 μm.

Outflow region
At the end of the outflow region there must be a boundary condition. The same Zou-He boundary as in the inflow region is used, but density is prescribed instead of velocity to calculate the unknown distributions. The density that is used for this is the initial density of the flow domain ρ 0 , therefore ρ 0 = ρ out . This boundary condition does not enforce a parabolic profile as expected in a tubular domain. Therefore, a flattening of the velocity profile in the outflow region is to be expected.
A particle is removed fully when one of its lsps crosses the last row of the outflow region (and consequently behavior of the particle becomes undefined). The pressure gradient from the fluid in the particle as opposed to the fluid outside it is approximately zero. Because of this there are no shockwaves resulting of the deletion of a particle. However, the deletion of a particle does influence the particle distribution and flow profile. For this reason the removal is done in the outflow region.
The size of the outflow region is dependent on the experiment and chosen boundary condition. It should be at least as long as the length of the longest possible elongated particle. When the boundary condition introduces artifacts in the outflow region the size of the outflow region should be at least as large to encompass all these artifacts.

Macroscopic fluid model
The fluid is simulated using the lattice Boltzmann method with the BGK collision operator. The lattice spacing is x = 1 μm . The fluid that is simulated is blood plasma at 37 • C , therefore, the kinematic viscosity is 1 . 28 × 10 −6 m 2 / s [15,16] and the density is 1025 kg/m 3 . The timestep is set to t = 1 . 176 × 10 −7 s . The relaxation time ( τ ) of the fluid is 1.1. With the particles present this gives an observed locally maximum compressibility error of 1% (as observed from our simulations). This means that the quasi- In the pre-inlet domain the fluid is driven by an external forcing term, which is implemented as proposed by Guo et al. [17] . In the beginning ( t = 0 s ) of a simulation, the body force is increased quadratically from zero to the final magnitude over one millisecond to allow the red blood cells to relax from their initial disk shape to their biconcave shape. A quadratic increase is chosen for numerical stability.

Immersed boundary implementation
The fluid-particle coupling is modeled through an immersed boundary method [7] . In this immersed boundary method the velocity of the fluid is interpolated to the particles at every timestep. The same is true for the force exercised by the material model of the particles on the fluid. This interpolation is done to the closest Eulerian points X with a discrete Dirac delta function as constructed by Mountrakis et al. [10] which effectively gives a fourpoint stencil with weights to the closest points. No repulsion particles are used at the fluid boundaries.

Pre-inlet
Firstly, five differently sized pre-inlet domains are tested to ensure that a correct distribution of particles enters the system. To this end five simulations 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 geometries ( L d = 200 μm ). The geometries only differ in length in the direction of the periodic boundary condition ( L p ), giving them different aspect ratios defined as the length of a channel divided by its width ( L p / L d ). The aspect ratios range from 0.5 to 2.5 . The rest of the parameters are the same ( Table 1 ).
All pre-inlets retrieve similar velocity profiles. Therefore, only the pre-inlet domain with aspect ratio 2.5 is plotted in Fig. 3 . Due to the shear thinning behavior of blood the velocity profile is more plug-like than the parabolic profile for a fluid without particle that   has the same average velocity. These finding are supported by Carboni et al. [18] .
In flowing whole blood the platelets will marginate to the cellfree layer. This is used as a measure of correctness for the preinlets. 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 of at most 50% of the mean hematocrit. This definition for the cellfree layer is used in the rest of this paper. In all the simulations of the pre-inlet domains the cell-free layer is 5 μm according to this definition (see Fig. 4 ). A platelet is considered marginated when it is present within this layer. Fig. 5 shows the percentage of marginated platelets over time. The initial margination happens quickly and then asymptotically tends to the final margination percentage, this finding is supported by the work of Pleunis et al. [19] and Zhao et al. [20] . All preinlets 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 finite size system effects. Therefore we have used the pre-inlet with the largest aspect ratio (2.5) in the following experiment.

Plane Poiseuille flow
To test the in-and outlet boundary conditions a plane Poiseuille flow is simulated with a pre-inlet domain of the same length as the flow domain (500 μm). See Table 2 for the parameters that are used. At the start ( t = 0 s ) of the simulation only the pre-inlet is filled with particles. The flow domain only consists of fluid. Fig. 6 shows the domain being filled up after turning on the body force   in the pre-inlet domain. After 0.01 s the flow domain equilibrates around the same hematocrit as the pre-inlet domain. In Fig. 7 the velocity magnitude of the fluid of the pre-inlet domain coupled with the flow domain is shown. The Zou-He velocity inlet attaches to the 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 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 profile towards the end of the outflow region as expected.
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 enough outflow region for this simulation. Fig. 8 shows the pre-inlet domain and flow domain with particles present. The red particles are red blood cells, while the black dots represent the platelets. The graphs on the bottom of Fig. 8 show the averaged hematocrit and platelet concentration for various locations in the flow domain and pre-inlet domain. Through the transition of the pre-inlet domain to the flow domain no noticeable 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 condition artifacts.

Same geometry with pre-inlet and periodic conditions
To show the necessity of our in-and outlet conditions we conduct two experiments with a single geometry. The geometry is such that it can be either run with periodic boundaries or a preinlet. 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 s with a velocity of 20 cm/s. In Fig. 10 we can see clear differences arising between the same geometry with different boundary conditions. The aneurysm with pre-inlet is enriched with more platelets than the aneurysm with periodic boundaries. This is due to depletion of platelets in the tube with the periodic conditions. The aneurysm with pre-inlet even has more platelets than are available in the whole periodic system after 0.23 s. Furthermore, the density of particles in the tube with periodic boundary conditions is lower in the aneurysm. We account this effect due to the non-changing number of particles in the whole periodic geometry. Both effects are unwanted and noticeable when running long simulations of such a geometry with periodic boundary conditions. The geometry that is simulated with our new in-and outlet conditions does not suffer from these effects.

Application to a curved aneurysm
As an example of applying the in-and outlet conditions for an application where periodic boundary conditions cannot be used we simulate the flow of blood within a curved vessel that has an attached aneurysm. Again, the pre-inlet domain that is attached to the curved aneurysm has a width of 200 μm and a length of 500 μm.
The domain is also randomly initialized with the same hematocrit and platelet percentage as the pre-inlet domain. Otherwise the flow domain (aneurysm) would take in the order of seconds to  Fig. 11 shows the aneurysm after one second of simulation. In this figure the location of a platelet is represented by a black dot. The red blood cells are omitted for clarity. On the background the velocity of the fluid and the streamlines are plotted. We can see that platelets start to accumulate in the aneurysm as predicted by Závodszky et al. [21] and observed in earlier work by Mountrakis et al. in a straight vessel geometry with periodic boundary conditions [22] .

Discussion and conclusion
We have shown that Lykov's et al. [6] method can be adapted for two dimensional suspensions simulated with IB-LBM. Furthermore, we show that these two dimensional suspensions can scale to larger inflow boundaries of 200 μm. The tests with the plane Poiseuille flow show that an accurate distribution of particles within the flow domain and a well behaving blood flow is retrieved with our in-and outlet boundaries. The simulation with the curved vessel and attached aneurysm demonstrates the importance of such in-and outlet boundary conditions. Since platelets are accumulating in the aneurysm it is necessary to have an inflow that is independent of the behaviour of particles in the system. These in-and outflow boundary conditions can be extended into three dimensions as a logical next step. As an extension to the curved aneurysm showcase it is possible to do particle residence time analysis [21] 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 detail, also taking into account the pulsatility of blood flow.
The material model used by Mountrakis et al. [10] suffered from instabilities in our simulations. Therefore, some constants and material model equations are adjusted.
In the simulations the red blood cells appeared to be too stiff. This is solved by squaring the angle used in calculating the bending resistance force ( f trs ). The bending constant C trs is decreased a hundredfold to 1 × 10 −9 N / rad .
We noticed numerical instability caused by large forces from the repulsion force ( F rep ). This seemed to be happening because of a exponential term for the distance. Therefore this term is removed and made linear. This dampens the force impact of collisions that with the quadratic term could create such large forces that the simulation became numerically unstable.
As a final step some parameters for the material model are slightly adjusted (see Table 3 ).
In the simulations it became apparent that the repulsion force causes the particles to keep a distance of h cutoff from each other at all times. Because of this the particles act as if they have a diameter that is 0 . 5 × h cutoff = 0 . 3 μm larger. Therefore this extra area is added to the calculation for the hematocrit. This is shown as the apparent area in Table 4 .