# A voxelized immersed boundary (VIB) finite element method for accurate and
efficient blood flow simulation

G. C. Bourantas<EMAIL_ADDRESS>B. F. Zwick  D. S. Lampropoulos  V. C.
Loukopoulos  K. Katsanos  A. A. Dimas  V. N. Burganos  A. Wittek  K. Miller

###### Abstract

We present an efficient and accurate immersed boundary (IB) finite element
(FE) method for internal flow problems with complex geometries (e.g., blood
flow in the vascular system). In this study, we use a voxelized flow domain
(discretized with hexahedral and tetrahedral elements) instead of a box
domain, which is frequently used in IB methods. The proposed method utilizes
the well-established incremental pressure correction scheme (IPCS) FE solver,
and the boundary condition-enforced IB (BCE-IB) method to numerically solve
the transient, incompressible Navier–Stokes flow equations. We verify the
accuracy of our numerical method using the analytical solution for the
Poiseuille flow in a cylinder, and the available experimental data (laser
Doppler velocimetry) for the flow in a three-dimensional 90° angle tube bend.
We further examine the accuracy and applicability of the proposed method by
considering flow within complex geometries, such as blood flow in aneurysmal
vessels and the aorta, flow configurations that would otherwise be difficult
to solve by most IB methods. Our method offers high accuracy, as demonstrated
by the verification examples, and high applicability, as demonstrated through
the solution of blood flow within complex geometry. The proposed method is
efficient, since it is as fast as the traditional finite element method used
to solve the Navier–Stokes flow equations, with a small overhead (not more
than 5$\%$) due to the numerical solution of a linear system formulated for
the IB method.

###### keywords:

Transient incompressible Navier–Stokes, Incremental pressure correction scheme
(IPCS), Immersed Boundary (IB), Internal flows, Voxelized flow domain

††journal: Computers & Fluids

##  1 Introduction

Computational fluid dynamics (CFD) simulations are becoming a reliable tool
for analysis of complex flow regimes in healthy and pathological blood
vessels. To accurately model such flows, often a scale-resolving turbulence
model is required, and high-order methods are used for discretizing and
numerically solving the governing flow equations. However, high-order
numerical methods often necessitates the use of structured grids, which may
not be able to handle a variety of complex geometries that arise in arterial
flow domains. Additionally, even for laminar flow regimes, where the flow
regime is less complex, computing derived quantities such as wall shear stress
(WSS), time average wall shear stress (TAWSS) and oscillatory shear index
(OSI), requires a high-quality mesh, which is extremely difficult to obtain.

The immersed boundary (IB) methods emerged as an attractive approach because
of its ability to efficiently handle complex moving and rotating geometries on
structured Cartesian grids, and to simulate flow past rigid and
moving/deforming bodies with irregular/complex shapes [1, 2]. IB methods, in
contrast to mesh-based methods, avoid tedious mesh generation, since they rely
on Cartesian grids to solve the governing flow equations, and on a discrete
set of points (often referred to as Lagrangian points) for the imposition of
prescribed velocity boundary conditions. The IB methods are particularly
efficient for moving objects and deforming geometries—they are attractive with
moving boundaries as they avoid re-meshing—where it does not require the
generation of a new mesh at each time step, but only the updated position of
the points describing the immersed boundary.

Immersed boundary method-based simulations are frequently used in external
aerodynamics problems [3, 4, 5]. However, their application to internal fluid
flow in complex geometries, such as blood flow in arteries, remains
challenging. Yokoi et al.[6] simulated blood flow in a cerebral artery with
multiple aneurysms using an immersed boundary method (they discretized the
flow domain using a Cartesian grid). Although Yokoi et al. didn’t mention the
percentage of nodes outside the flow domain over the total number of grid
nodes (they used 0.6 million Cartesian grid nodes in total), it was apparent
that a large portion of grid nodes were in the exterior of the fluid domain.
Delorme et al. [7] performed Large Eddy Simulation (LES) studies of powered
Fontan hemodynamics using a relatively short vena cava and long pulmonary
arteries to reduce the percentage of Cartesian grid nodes located in the
exterior of the flow domain. Despite that, the authors reported that a
significant number of grid nodes were located outside of the flow domain. To
extend IB methods to simulate complex arterial geometries, de Zélicourt et al.
[8] simulated blood flow in a geometry of the total cavo-pulmonary connection
(TCPC) using an unstructured Cartesian grid. However, studies that apply the
IB methods to complex internal flow configurations remain scarce. The key
reason for this is increase in the internal memory and computational efforts
for single block (structured and unstructured) Cartesian grids that are often
used in the Ib methods [8]. Additionally, handling complex geometries on a
structured grid, all the inflow and outflow boundaries of the geometry have to
terminate on the boundary faces of the global bounding box that encloses both
the fluid and solid regions. Even if this requirement could be met in certain
cases by truncating the complex geometry, in many cases it is not doable. To
overcome the aforementioned problems and extend the applicability of the IB
method to simulate blood flow in complex anatomies, authors in [9] proposed a
method that combined multi-block structured grid and IB method on an
inherently parallel framework. The method enabled the simulation of fluid flow
in complex geometries, reducing the percentage of grid nodes located outside
the immersed boundary.

In recent decades, a large class of IB methods based on the finite element
(FE) method has been developed. Löhner et al. [10, 11] applied the kinetic and
kinematic impositions of boundary conditions introduced in immersed boundary
methods [12] to adaptive nodal finite element grids. Glowinski et al. [13, 14]
simulated viscous flow interacting with rigid particles by using a distributed
Lagrange multiplier field to enforce the rigid body motion in each particle
subdomain onto the overlapping fluid field. Zhang and co-workers [15, 16]
proposed the immersed finite element method (IFEM) to use a flexible
Lagrangian solid mesh that moves on top of a background Eulerian fluid mesh,
which spans the entire computational domain. In IFEM the mesh generation is
greatly simplified, and the continuity between the fluid and solid sub-domains
is enforced by interpolating the velocities. Forces are distributed using the
reproducing kernel particle method (RKPM) delta function. In fact, the higher-
ordered RKPM delta function enables non-uniform spatial meshes (on arbitrary
geometries) for the fluid domain. Wang and Liu [17] proposed the extended
immersed boundary method (EIBM), by including the submerged elastic solid
occupying a finite volume in the flow domain. This was done by replacing the
kinematic and dynamic matching of the fluid-solid interface and the effect of
the immersed solid with nodal forces calculated in the context of finite
element formulations. The equations in both the fluid and solid domains are
approximated using the finite element method, while the continuity between the
fluid and solid domains is enforced by interpolating the velocities and the
distribution of forces delta function with the reproducing kernel particle
method (RKPM).

Finite cell method [18] is an example of well-established IB method. In the
finite cell method, the computational domain and the embedded domain are
discretized using a finite number of computational cells (usually having
simple shapes like squares or cubes) that cover the entire domain. The unknown
field variables are smoothly extended beyond the cells, cut by the original
boundary. Special integration techniques are employed to treat cells cut by
the boundary. The finite cell method (FCM) belongs to the fictitious domain
class [14, 19] of embedded domain methods [20, 10, 11]. The main idea of the
fictitious domain methods is that they extend the physical domain of interest
beyond the complex boundaries into a larger embedding domain discretized by a
simple structured grid. In the finite cell method, the computational domain
and the embedded domain are discretized using a finite number of computational
cells. Special integration techniques are employed to treat cells cut by the
boundary. The fictitious domain concept has been used in numerous immersed
boundary type methods, which have become increasingly popular over the last
decade for the simulation of large-scale computational fluid mechanics
problems. Some examples include the extended finite element method [21, 22,
23], isogeometric immersed boundary methods [24, 25, 26, 27, 28], or embedded
domain approaches focusing on Lagrange multipliers [29, 30], Nitsche type
methods [31, 32] and discontinuous Galerkin methods [33, 34].

In the present study, we develop and use the voxelized finite element (FE)
boundary condition-enforced immersed boundary (VBCE-IB) method, which is an
extension of the BCE-IB method introduced in [35], to numerically solve blood
flow in arteries (internal flow). The flow domain is discretized using
voxels/cubes that envelope the immersed object and it is no longer a box
domain. The main advantage of the VBCE-IB method over the BCE-IB method is
that a large percentage of the Eulerian grid nodes are located inside the
immersed boundary Additionally, the VBCE-IB method accurately satisfies both
the governing equations and boundary conditions using velocity and pressure
correction procedures. The velocity correction applies implicitly such that
the velocity on the immersed boundary (Lagrangian points) interpolated from
the corrected velocity values computed on the mesh nodes (Eulerian nodes)
accurately satisfies the prescribed velocity boundary conditions. The
incremental pressure correction scheme (IPCS) is a modified version of the
projection method, which provides improved accuracy at little extra
computational cost. We implement our method using FEniCS [36] (an open-source
software package that solves partial differential equations) to numerically
solve the incompressible Navier–Stokes equations. FEniCS encompasses software
to specify variational forms of PDEs in the high-level Unified Form Language
(UFL) and compiles these descriptions into numerical routines. This software
has been shown to achieve high levels of computational performance, even in
problems with billions of degrees of freedom [37]. The proposed IB scheme
applies to both external and internal flows, but in the present study we are
particularly interested in internal flow problems. Our numerical scheme
increases the computational efficiency of the IB method for internal flow
cases by reducing the number of unused elements through a sophisticated and
efficient mesh refinement inside the immersed boundary. Reducing the
computational cost is of paramount importance for IB methods applied to
internal flows because not all of the advantageous features of IB methods that
have made them popular for external flow simulations apply to internal flow
simulations.

The remainder of the paper is organized as follows. In Section 2, we present
the numerical formulation of the proposed immersed boundary (IB) method. In
Section 3, we solve a benchmark problem (Poiseuille flow) to verify the
accuracy of the proposed IB method. Section 4 presents numerical examples that
demonstrate the accuracy and applicability of the proposed IB method. Section
5 contains conclusions.

##  2 Numerical method

We consider the incompressible viscous flow in a spatial domain $\Omega$,
which contains an immersed boundary in the form of a closed surface $\partial
S$. The immersed boundary is modeled as localized body forces acting on the
surrounding fluid. The incompressible N-S equations—in their primitive
variables (velocity $\bm{u}$ and pressure $p$) formulation—accounting for the
immersed boundary are written as:

$\frac{\partial\bm{u}}{\partial
t}+\left(\bm{u}\cdot\bm{\nabla}\right)\bm{u}=-\bm{\nabla}p+\bm{\nabla}\cdot
2\nu_{f}\bm{\varepsilon}(\bm{u})+\bm{f},$ (1) $\bm{\nabla}\cdot\bm{u}=0,$ (2)

subject to the velocity boundary condition on $\partial S$

$\bm{u}(\bm{X}(s,t),t)=\bm{U}_{B},$ (3)

with $\bm{U}_{B}$ being the prescribed velocity of the immersed boundary, and
$v_{f}$ the kinematic viscosity of the fluid. The term
$\bm{\varepsilon}(\bm{u})$ is the strain-rate tensor defined as:

$\bm{\varepsilon}(\bm{u})=\frac{1}{2}\left(\bm{\nabla}\bm{u}+\left(\bm{\nabla}\bm{u}\right)^{\mathrm{T}}\right).$
(4)

The forcing term $\bm{f}$ is added to the right-hand side (RHS) of the
momentum Eq.(1) to account for the presence of the immersed boundary. The
forcing term $\bm{f}\left(\bm{x},t\right)$ is the local body force density at
the fluid nodes (Eulerian nodes), and it is obtained from the surface force
density $\bm{F}\left(s,t\right)$ ($s$ is a parametrization of the immersed
boundary surface) at the immersed boundary (Lagrangian points—we use
traditional terminology used in [38]). The forcing term is expressed as:

$\bm{f}(\bm{x},t)=\int_{\partial
S}\bm{F}(s,t)\,\delta\left(\bm{x}-\bm{X}(s,t)\right)\mathrm{d}s,$ (5)

where $\bm{x}$ and $\bm{X}\left(s,t\right)$ denote the Eulerian nodes and the
Lagrangian point coordinates that discretize the fluid domain and represent
the immersed boundary, respectively. Eulerian nodes and Lagrangian points
interact through the Dirac delta function
$\delta\left(\bm{x}-\bm{X}\left(s,t\right)\right)$. Therefore, the velocity at
the immersed boundary points is calculated from the velocity at the Eulerian
nodes using the Dirac delta function

$\bm{U}(\bm{X}(s,t))=\int_{\Omega}\bm{u}(\bm{x})\,\delta\left(\bm{x}-\bm{X}(s,t)\right)\mathrm{d}V.$
(6)

###  2.1 Incremental Pressure Correction Scheme (IPCS)

The Incremental Pressure Correction Scheme (IPCS) [39] is an operator
splitting method that decouples the pressure and velocity fields of the N-S
equations. IPCS is a modified version of the fractional step method proposed
by Chorin [40] and Témam [41]. It improves the accuracy of the original scheme
with little extra cost. In our description, we consider no external forces for
the N-S equations, and we explain in detail how the IB forces were
incorporated. The spatial discretization of velocity and pressure was
performed using linear Lagrange $(P_{1}/P_{1})$ elements [42, 43]. Since we
used a splitting approach (IPCS method), the Ladyzhenskaya-Babuska-Brezzi
(LBB) condition is inherently satisfied, and thus stabilization of the spatial
discretization was not required [43].

The IPCS scheme involves three steps. In the first step, we compute a
tentative velocity $\bm{u}^{*}$ by advancing the linear momentum Eq.(1) in
time using the backward Euler difference scheme,
$\dot{\bm{u}}=\frac{\bm{u}^{n+1}-\bm{u}^{n}}{\Delta t}$. We linearize the
nonlinear terms using a semi-implicit method, such that the advection term
$\bm{u}\cdot\bm{\bm{\nabla}}\bm{u}$ becomes
$\bm{u}^{n}\cdot\bm{\bm{\nabla}}\bm{u}^{n+1}$. Using this semi-implicit
approach for linearization, the Courant-Friedrichs-Lewy (CFL) condition, which
limits the time step according to the velocity and spatial discretization to
ensure stability of the solution, becomes less restrictive. Additionally, the
kinematic viscosity is written as $v_{f}=v_{f}\left(\bm{u}^{n}\right)$, which
simplifies the N-S equations that are now written as a linearized set of
equations known as the Oseen equations:

$\bm{u}^{n+1}+\Delta t\bm{u}^{n}\cdot\bm{\nabla}\bm{u}^{n+1}-\Delta
t\bm{\nabla}\cdot 2v_{f}\bm{\varepsilon}(\bm{u}^{n+1})+\Delta
t\bm{\nabla}p^{n+1}=\bm{u}^{n},$ (7) $\bm{\nabla}\cdot\bm{u}^{n+1}=0.$ (8)

In Eq.(7), $p^{n+1}$ is still unknown. Therefore, we will compute the
tentative velocity $\bm{u}^{*}$ ($\bm{u}^{*}\approx\bm{u}^{n+1}$), by
replacing in Eq.(7) $p^{n+1}$ with the known value $p^{n}$. This makes Eq.(7)
easier to solve for the tentative velocity $\bm{u}^{*}$:

$\bm{u}^{*}+\Delta t\,\bm{u}^{n}\cdot\bm{\nabla}\bm{u}^{n+1}-\Delta
t\,\bm{\nabla}\cdot 2v_{f}\bm{\varepsilon}\left(\bm{u}^{*}\right)+\Delta
t\,\bm{\nabla}p^{n}=\bm{u}^{n}.$ (9)

Equation (9) is numerically solved to compute the tentative velocity
$\bm{u}^{*}$ without using the incompressibility constraint Eq.(2).

In the second step, we use the tentative velocity $\bm{u}^{*}$ to compute the
updated velocity $\bm{u}^{n+1}$, which is divergence-free (and should fulfill
the incompressibility constraint). We define a function for the velocity
correction as $\bm{u}^{c}=\bm{u}^{n+1}-\bm{u}^{*}$, and after some algebra
(subtracting Eq.(9) from Eq.(7)) and using that
$\bm{\nabla}\cdot\bm{u}^{n+1}=0\Rightarrow\
\bm{\nabla}\cdot\bm{u}^{c}=-\bm{\nabla}\cdot\bm{u}^{*}$ (incompressibility
constraint) we obtain a new set of flow equations for $\bm{u}^{c}$:

$\bm{u}^{c}+\Delta t\,\bm{u}^{n}\cdot\bm{\nabla}\bm{u}^{c}-\Delta
t\,\bm{\nabla}\cdot 2v_{f}\bm{\varepsilon}\left(\bm{u}^{c}\right)+\Delta
t\,\bm{\nabla}\Phi^{n+1}=0,$ (10)
$\bm{\nabla}\cdot\bm{u}^{c}=-\bm{\nabla}\cdot\bm{u}^{*},$ (11)

where $\Phi^{n+1}=p^{n+1}-p^{n}$. Equation (10) can be further simplified to
$\bm{u}^{c}+\Delta\text{t\ }\bm{\bm{\nabla}}\Phi^{n+1}=\bm{0}$ and the system
of Eqs.(10)-(11) is reduced to an elliptic type (Poisson) partial differential
equation for the pressure difference $\Phi^{n}$

$\nabla^{2}\Phi^{n+1}=\frac{\bm{\nabla}\cdot\bm{u}^{*}}{\Delta t},$ (12)

with boundary conditions being those applied to the original flow problem.

In the third and final step, we compute the updated velocity $\bm{u}^{n+1}$
and pressure $p^{n+1}$ through $p^{n+1}=\Phi^{n+1}+p^{n}$ and
$\bm{u}^{n+1}=\bm{u}^{*}+\Delta t\bm{\nabla}\Phi^{n+1}$, respectively.

The algorithmic procedure for our implementation of the IPCS for solving N-S
equations is summarized below:

  1. 1.

Calculate the tentative velocity $\bm{u}^{*}$ by solving Eq.(9)

  2. 2.

Solve the Poisson equation
$\bm{}\nabla^{2}\Phi^{n+1}\bm{=}\frac{\nabla\cdot\bm{u}^{*}}{\Delta t}$
(Eq.(12)) to obtain $\Phi^{n+1}$

  3. 3.

Calculate the corrected velocity $\bm{u}^{n+1}=\bm{u}^{*}+\Delta
t\bm{\nabla}\Phi^{n+1}$ and pressure $p^{n+1}=\Phi^{n+1}+p^{n}.$

###  2.2 Boundary Condition Enforced Immersed Boundary (BCE-IB) method

In this section, we describe the boundary condition enforced immersed boundary
(BCE-IB) method for the numerical solution of the incompressible 3D
Navier–Stokes equations. The immersed boundary method has two steps:

  * 1.

_Predictor step_ : where we numerically solve the N-S equations to compute the
predicted velocity field $\tilde{\bm{u}}\left(\bm{x}\right)$ by disregarding
the body force terms in Eq.(1) ($\bm{f}=\bm{0}$)

$\frac{\tilde{\bm{u}}-\bm{u}^{n}}{\Delta
t}+\left(\bm{u}^{n}\cdot\bm{\nabla}\bm{u}^{n}\right)=+\Delta
t\bm{\nabla}p^{n}-\Delta t\bm{\nabla}\cdot 2v_{f}\bm{\varepsilon}(\bm{u}^{n})$
(13)

  * 1.

_Corrector step_ : where we account for the effect of body forces and update
the predicted velocity field $\tilde{\bm{u}}$ to the updated one
$\bm{u}^{n+1}$, which satisfies the boundary condition
$\bm{u}\left(\bm{X}\left(s,t\right),t\right)=\bm{U}_{B}$

$\frac{\bm{u}^{n+1}-\tilde{\bm{u}}}{\Delta t}=\bm{f}^{n+1}.$ (14)

In the predictor step, Eq.(13) is used to compute the predicted velocity field
$\tilde{\bm{u}}$ under the incompressibility constraint (mass conservation),
which couples the velocity and pressure. The predictor step accounts for the
first of the three steps in the IPCS method. In this step, we compute the
predicted velocity $\tilde{\bm{u}}$, which in this case is the tentative
velocity $\bm{u}^{*}$ (Step 1 in the IPCS algorithmic procedure). The
corrector step involves the evaluation of the unknown body force
$\bm{f}^{n+1}$ and the update of the predicted velocity $\tilde{\bm{u}}$ to
the physical one $\bm{u}^{n+1}$ (Steps 2 and 3 in the IPCS algorithmic
procedure).

Accurate and efficient computation of the unknown body force $\bm{f}^{n+1}$ is
an important advantage of the Boundary Condition Enforced Immersed Boundary
(BCE-IB) method. In this method, the body force $\bm{f}^{n+1}$ is equivalent
to a velocity correction, which applies implicitly such that the velocity
$\bm{U}\left(\bm{X}\left(s,t\right),t\right)$ on the Lagrangian points,
interpolated from the physical velocity $\bm{u}\left(\bm{x},t\right)$ computed
on the Eulerian nodes, equals to the prescribed boundary velocity $\bm{U}_{B}$
(Eq. (3)). The body force term $\bm{f}^{n+1}$ is computed using the equation

$\bm{U}_{B}^{n+1}\left(\bm{X}^{n+1}\right)=\int_{\Omega}\left(\tilde{\bm{u}}+\Delta
t\frac{\bm{f}^{n+1}}{\rho}\right)\delta\left(\bm{x}-\bm{X}^{n+1}\right)\mathrm{d}V,$
(15)

derived by substituting Eq.(14) to Eq.(6). The force density
$\bm{f}^{n+1}\left(\bm{x},t\right)$, is evaluated on the Eulerian nodes, and
it is computed by distributing (‘spreading’) the boundary force
$\bm{F}\left(\bm{X}\left(s,t\right),t\right)$ on the Lagrangian points through
the Dirac delta function $\delta\left(\bm{x-X}\left(s,t\right)\right)$ (see
Eq.(5)). Eq.(15) is written as

$\bm{U}_{B}^{n+1}\left(\bm{X}^{n+1}\right)=\int\left(\tilde{\bm{u}}+\Delta
t\frac{\int{\bm{F}^{n+1}\left(\bm{X}^{n+1}\right)\delta\left(\bm{x}-\bm{X}^{n+1}\right)\mathrm{d}s}}{\rho}\right)\delta\left(\bm{x}-\bm{X}^{n+1}\right)\mathrm{d}V,$
(16)

and $\bm{U}_{B}^{n+1}$ is no longer correlated with $\bm{f}^{n+1}$ but instead
with $\bm{F}^{n+1}$.

To compute the boundary force $\bm{F}^{n+1}$, we rewrite Eq.(16) in a discrete
form that results in an algebraic system of equations. We represent the
immersed boundary by a set of Lagrangian points
$\bm{X}_{i}=\left(X_{i},Y_{i},Z_{i}\right),\ i=1,2,\ldots,M$ and the flow
domain by the Eulerian nodes $\bm{x}_{j}=\left(x_{j},y_{j},z_{j}\right),\
j=1,2,\ldots,N$. In our analysis, the Eulerian mesh in the vicinity of the
immersed object has a uniform Cartesian grid structure with grid spacing $h$.
Furthermore, the Dirac function
$\delta\left(\bm{x}-\bm{X}\left(s,t\right)\right)$ is approximated by a
continuous kernel distribution $D\left(\bm{x}_{i}-\bm{X}^{j}\right)$ given as

$D_{ij}=D\left(\bm{x}_{i}-\bm{X}^{j}\right)=\left(\frac{1}{h}\delta\left(\frac{x_{i}-X^{j}}{h}\right)\right)\left(\frac{1}{h}\delta\left(\frac{y_{i}-Y^{j}}{h}\right)\right)\left(\frac{1}{h}\delta\left(\frac{z_{i}-Z^{j}}{h}\right)\right),$
(17)

with the kernel $\delta\left(r\right)$ proposed by Lai and Peskin [2]

$\delta(r)=\begin{cases}\begin{aligned}
\frac{1}{8}\left(3-2|r|+\sqrt{1+4|r|-4r^{2}}\right),&\quad&|r|&\leq 1\\\
\frac{1}{8}\left(5-2|r|-\sqrt{-7+12|r|-4r^{2}}\right),&\quad&1<|r|&\leq 2\\\
0,&\quad&|r|&>2\end{aligned}\end{cases}$ (18)

Using Eq.(5) for the force term,

$\bm{f}^{n+1}(\bm{x}_{j})=\sum_{i=1}^{M}{\bm{F}^{n+1}\left(\bm{X}_{i}^{n+1}\right)D_{h}^{ij}\Delta
S_{i}},$ (19)

Eq.(16) is written as

$\bm{U}_{B}^{n+1}(\bm{X}_{i}^{n+1})=\sum_{j=1}^{N}\tilde{\bm{u}}(\bm{x}_{j})D_{h}^{ij}h^{3}+\sum_{j=1}^{N}\sum_{k=1}^{M}\frac{\bm{F}^{n+1}(\bm{X}_{i}^{n+1})\Delta
t}{\rho}D_{h}^{kj}D_{h}^{ij}h^{3},$ (20)

where $\Delta S_{i}$ is the area of the _i th_ surface element (segment).
Eq.(20) forms a well-defined system of equations for the variables
$\bm{F}_{i}^{n+1}\bm{\ }\left(i=1,2,\ldots,M\right)$. Eq.(20) can be written
in matrix notation as

$\bm{A}_{\bm{F}}\bm{F}=\bm{B}_{\bm{F}},$ (21)

with $\bm{A}_{\bm{F}},\bm{F}$ and $\bm{B}_{\bm{F}}$ defined as

$\bm{A}_{\bm{F}}=\frac{\Delta t}{\rho}h^{3}\begin{bmatrix}D_{11}\Delta
S_{1}&D_{12}\Delta S_{1}&\dots&D_{1N}\Delta S_{1}\\\ D_{21}\Delta
S_{2}&D_{22}\Delta S_{2}&\dots&D_{2N}\Delta S_{2}\\\
\vdots&\vdots&\ddots&\vdots\\\ D_{M1}\Delta S_{M}&D_{M2}\Delta
S_{M}&\dots&D_{MN}\Delta
S_{M}\end{bmatrix}\begin{bmatrix}D_{11}&D_{12}&\dots&D_{1M}\\\
D_{21}&D_{22}&\dots&D_{2M}\\\ \vdots&\vdots&\ddots&\vdots\\\
D_{N1}&D_{N2}&\dots&D_{NM}\end{bmatrix},$ (22)
$\bm{B}_{\bm{F}}=\bm{U}-h^{3}\bm{D}^{\mathrm{T}}\tilde{\bm{u}}=\begin{bmatrix}\bm{U}_{1}\\\
\bm{U}_{2}\\\ \vdots\\\
\bm{U}_{M}\end{bmatrix}-h^{3}\begin{bmatrix}D_{11}&D_{21}&\dots&D_{N1}\\\
D_{12}&D_{22}&\dots&D_{N2}\\\ \vdots&\vdots&\ddots&\vdots\\\
D_{1M}&D_{2M}&\dots&D_{NM}\end{bmatrix}\begin{bmatrix}\tilde{\bm{u}}_{1}\\\
\tilde{\bm{u}}_{2}\\\ \vdots\\\ \tilde{\bm{u}}_{N}\end{bmatrix},$ (23)

$\bm{F}=\begin{bmatrix}\bm{F}_{1}\\\ \bm{F}_{2}\\\ \vdots\\\
\bm{F}_{M}\end{bmatrix},$ (24)

where $\bm{U}_{i}\ (i\bm{=}1,2,\ldots,M\bm{)}$, $\bm{F}_{i}\bm{\
(}i\bm{=}1,2,\ldots,M\bm{)}$ and ${\tilde{\bm{u}}}_{j}\bm{\
}\left(j=1,2,\ldots,M\right)$ are the abbreviations for
$\bm{U}_{B}^{n+1}\left(\bm{X}_{\bm{i}}^{n+1}\right)$,
$\bm{F}^{n+1}\left(\bm{X}_{\bm{i}}^{n+1}\right)$ and
$\tilde{\bm{u}}\left(\bm{x}_{j}\right)$, respectively. By solving the system
of equations (Eq.(22)) using a direct solver, we obtain the unknown boundary
force $\bm{F}_{\bm{i}}^{n+1}\bm{\ }\left(i=1,2,\ldots,M\right)$ at all
Lagrangian points. Boundary forces $\bm{F}_{\bm{i}}^{n+1}$ are then
substituted into Eq.(19) and Eq.(16) to calculate the body force
$\bm{f}^{n+1}$ on the Eulerian nodes, and the corrected physical velocity
$\bm{u}^{n+1}$_._

###  2.3 Algorithmic procedure

Our approach combines the IPCS (section 2.1) with the BCE-IB method (section
2.2). The velocity field predicted by IPCS is corrected to account for the
force contribution from the immersed boundary. The proposed solution procedure
can be summarized as follows; to update the solution from time instance $n$ to
$n+1$:

  1. 1.

Calculate the tentative velocity $\tilde{\bm{u}}$ using Eq.(13)

  2. 2.

Compute the matrix $\bm{A}_{\bm{F}}$ using Eq.(22)

  3. 3.

Solve the system (Eq.(21)) to compute the boundary force
$\bm{F}^{n+1}\left(\bm{X}_{i}^{n+1}\right)\bm{,\ (}i\bm{=}1,2,\ldots,M\bm{)}$
on the Lagrangian points. Substitute the computed boundary force
$\bm{F}^{n+1}\left(\bm{X}_{i}^{n+1}\right)$ into Eq.(19) to obtain the body
force $\bm{f}^{n+1}\left(\bm{x}_{j}^{n+1}\right)\bm{,\
(}j\bm{=}1,2,\ldots,N\bm{)}$ on the Eulerian points.

  4. 4.

Use $\tilde{\bm{u}}$ to update the corrected velocity $\bm{u}^{(n+1)}=\Delta
t\bm{f}^{n+1}+\tilde{\bm{u}}$ using Eq.(14)

  5. 5.

Solve the Poisson equation
$\bm{}\nabla^{2}\Phi^{n+1}\bm{=}\frac{\bm{\nabla}\cdot\bm{u}^{n+1}}{\Delta t}$
for $\Phi^{n}$ (Eq.(12))

  6. 6.

Recompute the updated velocity $\bm{u}^{n+1}=\bm{u}^{n}+\Delta
t\bm{\bm{\nabla}}\Phi^{n+1}$ and pressure $p^{n+1}=\Phi^{n+1}+p^{n}.$

  7. 7.

Set physical/updated velocity $\bm{u}^{\left(n+1\right)}$ as
$\bm{u}^{\left(n\right)}$ and repeat steps (1) to (4) until the desired
solution is achieved.

Steps 1, 5, 6 apply to the IPCS method without the presence of the immersed
boundary, while Steps 2, 3, 4 calculate the velocity field (on the Eulerian
nodes) due to the presence of the immersed boundary.

##  3 Algorithm verification

In this section, we demonstrate the accuracy of the proposed method in
internal flow cases. We consider two benchmark flow problems, Hagen-Poiseuille
flow in a straight tube and flow in a U-bend geometry. The numerical results
obtained were compared against the analytical solution for the Hagen-
Poiseuille flow example and against experimental data in the U-bend geometry
flow case.

###  3.1 Poiseuille flow in a cylindrical tube

In the first benchmark flow example, we consider the Hagen-Poiseuille flow in
a straight cylindrical tube. The driving force of the flow is the pressure
difference applied to the inlet and outlet. In our simulations, we set the
density of the fluid to $\rho=1,050~{}\mathrm{kg\,m^{-3}}$ and the dynamic
viscosity to $\mu=0.00345~{}\mathrm{N\,s/m^{2}}$. The tube has length $L=10R$
and radius $R=0.005$ m, and is immersed in the voxelized flow domain, as shown
in Fig. 1. We demonstrate the accuracy of the voxelized BCE-IB method by
comparing our numerical findings with the analytical solution (seec[35]). The
immersed boundary (cylinder) is represented by a set of surface points
(Lagrangian nodes). The voxelized flow domain is discretized first with
hexahedral elements and second with tetrahedral elements (generated by
splitting the hexahedral elements). At the inlet (red circle in Fig. 1b) and
outlet (orange circle in Fig. 1b), we apply pressure boundary conditions,
while at the remaining walls (surface yellow Fig. 1b), we apply no-slip
velocity boundary conditions ($\bm{u}=0$).

Figure 1:  Poiseuille flow cylinder geometry (a) immersed boundary given as
STL file with nodes and facets (b) flow domain for the Poiseuille flow
example, with inlet (circle in red color), outlet (circle in orange color) and
walls (yellow color) boundaries; (c) immersed boundary into the voxelized flow
domain, and (d) voxelized flow domain along with the cross-section of the
mesh.

We do not conduct a mesh-independence analysis in order to ensure a mesh-
independent numerical solution. Instead, we use a dense mesh for the voxelized
flow domain. As previously demonstrated (see [35]), we set the voxelized
Eulerian mesh resolution based on the distribution of the Lagrangian nodes. We
set the Eulerian mesh resolution $h$ using the mean value $\Delta S_{mean}$ of
the surface area $\Delta S$ assigned to each Lagrangian point, such that
$h=\sqrt[]{\Delta S_{mean}}$. Furthermore, we examine how the accuracy of the
numerical results is affected (boundary effects) by the dimensions (radius) of
the voxelized flow domain (the number of hexahedral elements (cubes) located
outside the immersed boundary (see Fig. 2)). In our simulations we use five,
six and seven rows of elements outside the immersed boundary (see Table 1).
The time step was set to $dt=2.5\times 10^{-3}$ (this time step ensures the
accuracy and stability of the numerical solution). The simulation ends when
the normalized root mean square difference

$\text{NRMSE}=\frac{\sqrt{\frac{\sum_{i=1}^{N}\left(u_{i}^{t+dt}-u_{i}^{t}\right)^{2}}{N^{2}}}}{\left(\max\left(u_{i}^{t}\right)-\min\left(u_{i}^{t}\right)\right)},$
(25)

between two successive time instances for all three velocity components is
less than $10^{-6}$ (practically, the flow reaches a steady state). The exact
solution for the velocity field components is $u_{y}=u_{z}=0$ and
$u_{x}\left(r\right)=U_{\max}\left(1-\left(\frac{r}{R}\right)^{2}\right)$,
with $U_{\max}=-\frac{R^{2}}{4\mu}\frac{\text{dp}}{\text{dz}}$ being the
maximum velocity. By setting the pressure difference to
$\frac{\text{dp}}{\text{dz}}=\frac{1}{L}$ (applying $p_{\text{inlet}}=1$ and
$p_{\text{outlet}}=0$) the Reynolds number becomes $\mathrm{Re}=\frac{\rho
U_{\max}L}{\mu}\approx 551$.

We use an immersed boundary consisting of 4,869 nodes and 9,570 facets. The
mean area of the facets is $\Delta S_{mean}=2.5\times 10^{-7}$ $m^{2}$, and
the voxelized Eulerian mesh (flow domain) resolution is $h=\sqrt[]{\Delta
S_{mean}}=0.0005$ m. We discretize the voxelized flow domain using hexahedral
elements. We use five, six and seven rows of elements outside the immersed
boundary (Fig.2), which leads to 57,626 elements (63,761 nodes), 62,853
elements (69,264 nodes), and 67,570 elements (74,245 nodes), respectively. The
numerical results computed using the different voxelized flow domains were
interpolated on a body-fitted mesh (38,756 nodes and 213,913 elements) with
tetrahedral elements (used only for visualization).

Table 1 lists the NRMSE for the different voxelized flow domains (with
hexahedral elements) compared with the analytical solution computed on the
body-fitted tetrahedral mesh. The results show that the dimensions of the
voxelized domain do not affect the accuracy of the numerical solution. To
demonstrate the applicability of the proposed scheme, we used tetrahedral
elements to discretize the flow domain. We use the same geometrical specifics
used in the tetrahedral meshes above, and the hexahedral elements are now
split into tetrahedral elements, generating tetrahedral meshes with 345,756
elements (63,761 nodes), 377,118 elements (69,264 nodes) and 405,420 elements
(74,245 nodes), respectively. Table 2 lists the NRMSE for the different
voxelized flow domains (with tetrahedral elements) compared with the
analytical solution computed on the body-fitted tetrahedral mesh. The results
show that the dimensions of the voxelized domain do not affect the accuracy of
the numerical solution.

Table 1: Normalized root mean square error (NRMSE) of the voxelized immersed boundary method (with hexahedral elements) numerical results against the analytical solution omputed on the body-fitted tetrahedral mesh for the Poiseuille flow example. The three hexahedral meshes correspond to five, six and seven rows of elements outside the immersed boundary.  mesh | $h$ | number of nodes | number of elements |  $u_{x}$ \- NRMSE |  $u_{y}$ \- NRMSE |  $u_{z}$ \- NRMSE
---|---|---|---|---|---|---
$1$ | $5\times 10^{-4}$ | 63,761 | 57,626 | $1.33\times 10^{-2}$ | $5.52\times 10^{-4}$ | $3.77\times 10^{-4}$
$2$ | $5\times 10^{-4}$ | 69,264 | 62,853 | $2.42\times 10^{-2}$ | $4.63\times 10^{-4}$ | $3.56\times 10^{-4}$
$3$ | $5\times 10^{-4}$ | 74,245 | 67,570 | $1.96\times 10^{-2}$ | $3.23\times 10^{-4}$ | $2.18\times 10^{-4}$
Table 2: Normalized root mean square error (NRMSE) of the immersed boundary method (with tetrahedral elements) numerical results against the analytical solution omputed on the body-fitted tetrahedral mesh for the Poiseuille flow examplefor the Poiseuille flow example. The three tetrahedral meshes correspond to five, six and seven rows of elements outside the immersed boundary.  mesh | $h$ | number of nodes | number of elements |  $u_{x}$ \- NRMSE |  $u_{y}$ \- NRMSE |  $u_{z}$ \- NRMSE
---|---|---|---|---|---|---
$1$ | $5\times 10^{-4}$ | 63,761 | 345,756 | $1.29\times 10^{-2}$ | $4.02\times 10^{-4}$ | $3.77\times 10^{-4}$
$2$ | $5\times 10^{-4}$ | 69,264 | 377,118 | $1.31\times 10^{-2}$ | $3.12\times 10^{-4}$ | $3.56\times 10^{-4}$
$3$ | $5\times 10^{-4}$ | 74,245 | 405,420 | $1.15\times 10^{-2}$ | $2.96\times 10^{-4}$ | $2.46\times 10^{-4}$
Figure 2:  Poiseuille flow example; immersed boundary (green cylinder) and a
cross-section of the voxelized mesh (hexahedral elements) with (a) five (b)
six, and (c) seven elements outside the immersed boundary.

Figure 3 shows the velocity streamlines for $\mathrm{Re}=551$ produced on the
tetrahedral body-fitted mesh. The velocity field computed with the IB method
on the Eulerian nodes has been projected on the visualization mesh for
visualization purposes (the visualization mesh is not used in the computation
of the velocity and pressure field).

Figure 3:  Velocity streamlines for $\mathrm{Re}=551$ using a visualization
mesh of 213,913 linear ($P_{1}/P_{1}$) elements for the Poiseuille flow
example. The maximum velocity computed for $\mathrm{Re}=551$ was 0.0361 (the
theoretical value is 0.0362).

###  3.2 3D angle tube bend

In the second benchmark problem, we demonstrate the accuracy of the proposed
method considering the flow in a three-dimensional 90° angle tube bend, as
shown in Fig.4. With this flow example we demonstrate that the proposed
numerical IB method can accurately simulate flows in tortuous tubes where
secondary flows can occur.

For the U-bend flow problem, experimental flow data are available. Van De
Vosse et al. [44] performed laser Doppler velocimetry experiments to obtain
the center-plane axial velocities at a set of different angles around the tube
bend. To reduce flow disturbances due to inflow boundary conditions, we
extended the inlet and outlet by $R_{D}=4$ mm (the length of the inner
diameter $D_{i}$). In our numerical simulations, we set the Reynolds number —
computed based on the tube inner diameter $D_{i}$ and the mean inlet velocity
$\overline{U}$ — to $\mathrm{Re}=300$, to match the experimental conditions.
The inner diameter $D_{i}$ and the curvature radius $R$ of the tube were set
to 4 and 24 mm, respectively.

Figure 4:  U-bend geometry (a) top and (b) side view. The points on the
surface define the Lagrangian points of the immersed boundary $R_{D}=1\ mm$,
$R=24\ mm$, $D_{i}=4\ mm$.  Figure 5:  U-bend flow example (a) Immersed
boundary given as STL file with nodes and facets (b) Flow domain for the
Poiseuille flow example, with inlet (circle in red color), outlet (circle in
blue color) and walls (yellow color) boundaries; (c) Immersed boundary into
the voxelized flow domain, and (d) voxelized flow domain along with the cross-
section of the tetrahedral mesh.

The U-bend geometry (immersed object) is embedded within the voxelized flow
domain (Fig.5a). We use a parabolic velocity profile
($u_{y}(r)=U_{\max}\left(1-\left(\frac{r}{R}\right)^{2}\right)$,
$U_{\max}=0.122625~{}\mathrm{m/s}$) at the inlet (red circle in Fig.5b), while
at the outlet (blue circle in Fig.5b) we apply zero pressure
$\left(p=0\right)$. At the remaining boundaries (shown in yellow in Fig.5b) of
the flow domain, we set no-slip boundary conditions.

We do not apply a mesh-independence analysis study to ensure a mesh-
independent numerical solution. Instead, we use a dense mesh to discretize the
voxelized flow domain. The resolution $h$ of the uniform voxelized mesh
(Eulerian mesh) is based on the average area $\Delta S$ (and consequently on
the average length) of the surface elements (facets) of the immersed boundary.
We use tetrahedral elements to discretize the flow domain. We use an immersed
boundary consisting of 15,810 nodes and 31,462 facets. The mean area of the
facets is $\Delta S_{mean}=2.5\times 10^{-7}$ m2 and therefore the voxelized
Eulerian mesh (flow domain) resolution is $h=\sqrt[]{\Delta S_{mean}}=0.0005$
m. We use six rows of elements outside the immersed boundary. The voxelized
flow domain consists of 3,162,114 linear elements ($P_{1}/P_{1}$) and 558,401
nodes. The numerical results computed using the voxelized flow domain were
interpolated on a body-fitted mesh (used only for visualization having 143,406
nodes and 808,984 tetrahedral elements). We compared the results obtained
using body-fitted mesh and voxelized flow domain. The NRMSE is $2\times
10^{-3}$ (maximum value for all three velocity components), highlighting that
the proposed method with the voxelized flow domain is accurate, and decreases
the computational cost for the IB method for internal flows.

Fig.6 shows a comparison of the axial velocity computed using the IB method,
and the experimental data (red dots) at the tube outlet. The axial velocity is
computed on the center plane at the start of the outlet extension (see
Fig.4a), where excellent agreement between the velocity predicted using the
proposed IB method and experimental results is observed. The numerical results
obtained using the linear elements ($P_{1}/P_{1}$) are also in good agreement
with the experimental data, highlighting the accuracy of the proposed method.

Figure 6:  Axial velocity computed with the proposed IB method (dashed line)
compared with the experimental data (red dots) Van De Vosse et al. [44] at the
tube outlet for the U-bend flow problem.

##  4 Numerical results - Applications to complex flow domains

In this section, we demonstrate the accuracy and applicability of the proposed
method. We highlight the use of a voxelized geometry as the flow domain in
terms of robustness, and we demonstrate the applicability of our method in
internal flow problems with complex geometries. In the numerical examples
considered in this section, we chose to use linear $(P_{1}/P_{1})$ tetrahedral
elements (instead of hexahedral elements) to discretize the voxelized flow
domain. Even if the number of tetrahedral elements is higher than the number
of hexahedral elements (for the same voxelized flow domain), local refinement
using tetrahedral elements is more straightforward (there are no hanging
nodes) compared to hexahedral elements.

###  4.1 Flow in a blood vessel bifurcation

In this example, we simulate the blood flow in a bifurcation (Fig.7a). The
bifurcation geometry (immersed boundary) is represented by a triangle mesh and
given as an STL file. The inlet and two outlets of the vessel geometry were
extruded to ensure that the flow is fully developed. The vessel has length
$L=0.0296$ m, while the inlet has a mean radius of $R_{in}=0.0034$ m and the
two outlets of $R_{out}=0.0023$ m. The voxelized flow domain (Fig.7b)
envelopes the immersed boundary, as shown in Fig.7c.

Figure 7:  Flow in a vessel bifurcation example (a) immersed boundary given as
STL file (b) voxelized flow domain (c) immersed boundary enveloped into the
voxelized flow domain (d) flow domain for the Poiseuille flow example, with
inlet (circle in red color), outlet (circle in blue color) and walls (yellow
color) boundaries (e) voxelized flow domain along with the cross-section of
the hexahedral mesh (f) voxelized flow domain along with the cross-section of
the tetrahedral mesh.

We do not apply a mesh-independence analysis to ensure a mesh-independent
numerical solution. Instead, we use a dense mesh for the voxelized flow
domain. We used 4,601,112 linear $(P_{1}/P_{1})$ tetrahedral elements (816,570
nodes) to discretize the voxelized flow domain (the voxelized mesh was
generated using an in-house Python code). The dense mesh of 4,601,112
tetrahedral elements ensures a grid-independent numerical solution (in a
previous study, a mesh of 627,040 tetrahedral elements ensured a grid-
independent solution). The resolution $h$ of the uniform voxelized mesh
(Eulerian mesh) is based on the average area $\Delta S_{mean}$ (and
consequently on the average length) of the surface elements (facets) of the
immersed boundary. We used a resolution of $h=0.0005$ m, which corresponds to
$\Delta S_{mean}=2.5\times 10^{-7}$ m2, and seven elements outside the
immersed boundary, to eliminate flow instabilities near the flow domain
boundaries due to the presence of the immersed boundary. The mesh generation
is a straightforward and efficient procedure since it takes less than 60
seconds to generate the voxelized mesh with 4,601,112 linear tetrahedral
elements using a laptop computer with an i7 quad-core processor and 16 GB of
internal memory.

We set the total time for the simulation to $T=3.2$ sec (four cardiac cycles),
and the time step to $dt=5\times 10^{-3}$ sec. At the inlet (red area located
at the bottom surface in Fig.7d), we apply the velocity waveform shown in
Fig.8 (we use constant velocity with components only in the $z-$ axis), and
zero pressure boundary conditions at the two outlets (green area on the top
surface in Fig.7d). At the remaining walls (yellow surface in Fig.7d), we
apply no-slip boundary conditions. We use the Newtonian model for blood flow,
with dynamic viscosity $\mu=0.00345Pas$ and density of
$\rho=1,056\frac{kg}{m^{3}}$. To demonstrate the accuracy of our IB method, we
compare the results obtained using the current method with those for a body-
fitted mesh with tetrahedral elements. The body-fitted mesh consists of
154,681 nodes and 860,818 linear tetrahedral elements
$\left(P_{1}/P_{1}\right)$.

Table 3 lists the normalized root mean square error

$\text{NRMSE}=\frac{\sqrt{\frac{\sum_{i=1}^{N}\left(u_{i}^{IB}-u_{i}^{bodyfitted}\right)^{2}}{N^{2}}}}{\left(\max\left(u_{i}^{bodyfitted}\right)-\min\left(u_{i}^{bodyfitted}\right)\right)}$
(26)

at different time instances. The results obtained using our IB method and
body-fitted mesh reported in Table 3 are, for practical purposes,
indistinguishable as in the bioengineering applications differences of less
than 5$\%$ would be overwhelmed by uncertainties in biomechanical properties.
The decisive advantage of our method is the ease of generating patient-
specific computational grids with practically no increase in computational
cost. The only computational overhead of our IB method in comparison to the
approach relying on a body-fitted mesh is due to the numerical solution of the
linear system given by Eq.(21). For the computation of flow in the coronary
artery bifurcation conducted here, this overhead was negligible (less than
1$\%$ of the total computation time). Fig.9 shows the streamlines in the
coronary artery bifurcation at time $t=2.7$ sec and $0.25$ sec (peak systole).

Figure 8:  (a) Pulsatile velocity waveform (shown for three cardiac cycles) and (b) the pulsatile velocity waveform (one cardiac cycle) imposed at the inlet for the flow in a vessel bifurcation example.  Table 3: NRMSE of the voxelized flow domain using tetrahedral and hexahedral elements at different time instances for the flow in a vessel bifurcation example (see Fig.7 and Fig.9).  velocity component | $t=2.5s$ | $t=2.6s$ | $t=2.7s$ | $t=2.8s$ | $t=2.9s$ | $t=3.0s$ | $t=3.1s$ | $t=3.2s$
---|---|---|---|---|---|---|---|---
$u_{x}$ | 2.29 | 3.78 | 1.78 | 3.18 | 4.27 | 2.51 | 2.25 | 2.22
$u_{y}$ | 3.05 | 4.04 | 1.93 | 3.71 | 3.97 | 3.16 | 3.64 | 3.88
$u_{z}$ | 4.53 | 4.18 | 3.60 | 2.82 | 3.78 | 4.27 | 3.33 | 2.82
Figure 9:  Streamlines (left) and pressure (right) at (a) peak systole
pressure and (b) diastolic minimum for the coronary artery bifurcation flow
example.

###  4.2 Flow in the descending aorta with aneurysm

In this example, we simulate blood flow in the descending aorta with an
aneurysm shown in Fig.10a. The immersed boundary (descending aorta) is 25 cm
long in the $z$\- axis and is embedded within a voxelized geometry that
envelopes the immersed boundary (Fig.10c).

Figure 10:  Flow in the descending aorta with aneurysm example (a) immersed
boundary given as STL file (b) voxelized flow domain (yellow color) (c)
immersed boundary along with the cross-section of the tetrahedral mesh (d)
flow domain for the descending aorta with aneurysm flow example, with inlet
(circle in red color), outlet (circle in green color) and walls (yellow color)
boundaries.

The voxelized flow domain is discretized using linear
$\left(P_{1}/P_{1}\right)$ tetrahedral elements, that envelopes the immersed
boundary. The immersed boundary consists of 116,218 nodes and 232,127
triangular surface elements (facets). The mesh resolution $h$ of the voxelized
flow domain is computed based on the average area $\Delta S_{mean}$ of the
facets of the immersed boundary. The mean area of the facets is $\Delta
S_{mean}=1.06\times 10^{-7}$ m2 and therefore, the mesh resolution of the
voxelized flow domain is $h=\sqrt[]{\Delta S_{mean}}=0.00035$ m. We use six
rows of elements outside the immersed boundary. We discretize the flow domain
using 12,718,206 ($P_{1}/P_{1}$) elements and 3,558,401 nodes. This mesh,
which ensures a grid-independent numerical solution, was generated in less
than $30$ sec.

Figure 11:  (a) Pulsatile pressure waveform (shown for four cardiac cycles)
and (b) the pulsatile pressure waveform (one cardiac cycle) imposed at the
inlet for the flow in the descending aorta with aneurysm example.

We simulate four cardiac cycles, and we set the total time of the simulation
to $T=3.44$ sec and the time step to $dt=1\times 10^{-4}$ sec. We apply the
time-dependent pressure profile shown in Fig.11 at the inlet (red area located
at the top surface in Fig. 10d) and zero pressure boundary conditions at the
two outlets (green areas on the top surface in Fig.10d). At the remaining
walls (yellow surface in Fig.10d), we apply no-slip boundary conditions. We
use the Newtonian model for blood flow, with dynamic viscosity
$\mu=0.00345Pas$ and density of $\rho=1,056\frac{kg}{m^{3}}$. We visualize the
results for the voxelized immersed boundary method simulation by projecting
the numerical solution computed on the voxelized mesh to a body-fitted mesh
with tetrahedral elements (this mesh is only for visualization and is not used
in the simulation). We use the visualization body-fitted tetrahedral mesh (in
total 4,018,371 elements) to numerically solve the flow equations using the
same boundary conditions, and we compare the numerical solution (velocity
field components) of the voxelized immersed boundary solver projected on the
body-fitted mesh tetrahedral mesh nodes with the numerical solution obtained
using the body-fitted mesh tetrahedral mesh at different time instances in the
cardiac cycle. Table 4 lists the normalized root mean square error (NRMSE)
between the two solutions (IB projected and body-fitted) on the body-fitted
tetrahedral mesh vertices at different time instances. The results obtained
using our IB method and body-fitted mesh reported in Table 4 are, for
practical purposes, indistinguishable as in the bioengineering applications;
differences of less than 5$\%$ would be overwhelmed by uncertainties in
biomechanical properties.

Table 4: NRMSE of the voxelized flow domain using tetrahedral elements at different time instances for the descending aorta with aneurysm flow example (see Fig.10 and Fig.12).  $NRMSE$ | $t=0.1s$ | $t=0.2s$ | $t=0.3s$ | $t=0.4s$ | $t=0.5s$ | $t=0.6s$ | $t=0.7s$ | $t=0.8s$
---|---|---|---|---|---|---|---|---
$u_{x}$ | 1.53 | 2.38 | 2.03 | 1.94 | 3.91 | 1.88 | 1.84 | 1.74
$u_{y}$ | 2.33 | 3.38 | 2.43 | 2.34 | 2.94 | 2.38 | 3.42 | 3.13
$u_{z}$ | 3.13 | 2.48 | 2.26 | 3.24 | 3.44 | 2.82 | 3.31 | 2.28
Figure 12:  Streamlines (left) and pressure (right) at (a) $t=1.92$ and (b)
$t=2.32$ sec for the flow in the descending aorta with aneurysm example.

Figure 12 shows the pressure contours and the velocity streamlines at time
$t_{1}=1.92$ sec and at $t_{2}=2.32$ sec. The results plotted are the
numerical results computed on the voxelized mesh and then projected on the
body-fitted tetrahedral mesh.

###  4.3 Blood flow in the ascending aorta

In the final example, we simulate blood flow in the ascending aorta (the
brachiocephalic, left common carotid, and left subclavian arteries are
included in the model) shown in Fig. 13a.

Figure 13:  Blood flow in the ascending aorta example (a) immersed boundary
given as STL file (b) voxelized flow domain (yellow color) (c) immersed
boundary along with the cross-section of the tetrahedral mesh (d)-(e) immersed
boundary along with the cross-section of the tetrahedral mesh and (f) flow
domain for the descending aorta with aneurysm flow example, with inlet (circle
in red color), outlets (circle in green color) and walls (yellow color)
boundaries.

The voxelized flow domain (we use linear $\left(P_{1}/P_{1}\right)$
tetrahedral elements) envelopes the immersed boundary. The immersed boundary
consists of 15,938 nodes and 31,684 triangular surface elements (facets). We
compute the mesh resolution $h$ of the voxelized flow domain based on the
average area $\Delta S_{mean}$ (and consequently on the average length) of the
surface elements (facets) of the immersed boundary. The mean area of the
facets is $\Delta S_{mean}=5.6\times 10^{-7}$ m2 and, therefore, the mesh
resolution of the voxelized flow domain is $h=\sqrt[]{\Delta
S_{mean}}=0.00075$ m. We use six rows of elements outside the immersed
boundary. We discretize the flow domain with 3,162,114 ($P_{1}/P_{1}$)
elements and 558,401 nodes. This mesh ensures a grid-independent numerical
solution. The voxelized mesh was generated in less than $30$ sec using a
laptop computer with an i7 quad-core processor and 16 GB internal memory.

We simulate four cardiac cycles, and we set the total time of the simulation
to $T=3.44$ sec and the time step to $dt=10^{-4}$ sec. We apply a time-
dependent pressure profile at the inlet (circle in green color in Fig.13f); in
the five outlets (circles in red color in Fig.13f), we apply zero pressure;
and at the remaining boundaries (yellow color in Fig.13f) we apply no-slip
boundary conditions. We use the Newtonian model for blood flow, with dynamic
viscosity $\mu=0.00345Pa\cdot s$ and density of $\rho=1,056\frac{kg}{m^{3}}$.
To demonstrate the accuracy of our IB method, we compare the results obtained
using the current method—after projecting the numerical solution on a body-
fitted mesh—with those for a body-fitted mesh with tetrahedral
$\left(P_{1}/P_{1}\right)$ elements (92,786 nodes and 504,628 elements). The
comparison applied at different time instances (every $0.1$ s interval) of the
cardiac cycle indicated a normalized root mean square error (NRMSE) of less
than $5\%$ for the three velocity components.

We visualize the results for the voxelized immersed boundary method simulation
by projecting the numerical solution computed on the voxelized mesh to a body-
fitted mesh used to highlight the accuracy of the proposed method.
Additionally, we have conducted a qualitative comparison of the predicted
streamlines and pressure contours at the time instances of $t=1.92$ sec and
$t=2.32$ sec (Fig.14). There are no visually distinguishable differences
between the results obtained using our IB method and approach using body-
fitted mesh.

Figure 14:  Pressure (left) and streamlines (right) at (a) $t=1.92$ and (b)
$t=2.32$ sec for the blood flow in the ascending aorta example.

##  5 Conclusions

In this contribution, we presented a voxelized immersed boundary method for
internal flows with application to blood flow simulations. The proposed scheme
combines the finite element (FE) method to numerically solve the Navier–Stokes
equations and the boundary condition enforced immersed boundary (BCE-IB)
method to account for complex geometries. In this study, we focused on mesh
generation and did not consider image analysis (such as segmentation and
geometry reconstruction) of the blood vessels used in the simulations. The
surface of blood vessels can be automatically generated from computed
tomography angiography (CTA) using algorithms available in existing image
processing tools such as 3D Slicer.

The proposed method increases the computational efficiency of the IB method
for internal flows by decreasing the time to mesh generation and by reducing
the number of unused mesh points through a sophisticated and efficient method
for voxelized mesh generation in the enveloping of the immersed boundary. This
method facilitated the generation of a computational mesh consisting of over
three and a half million elements (computation of flow in an aortic
bifurcation—see Fig. 10) in less than 30 seconds on an off-the-shelf laptop
with a quad-core i7 processor and 16 GB of internal memory. These promising
results make it possible to regard our flow solver as compatible with the time
constraints of clinical workflows.

The proposed method utilizes an efficient and accurate finite element solver
based on the incremental pressure correction scheme (IPCS). The accuracy of
the proposed scheme was successfully verified through comparison with an
analytical solution for the Poiseuille flow in a cylinder with a circular
cross-section. We further validated the accuracy of the proposed method by
solving a flow case where experimental data are available. The applicability
and efficiency of the proposed method was demonstrated through the solution of
flow examples with complex geometry. We solved several problems related to the
IB and FE methods. In the former, for internal flows, we use a flow domain in
which most of the elements are located inside the immersed boundary, while for
the latter, we use a high-quality mesh.

## Acknowledgments

A. Wittek and K. Miller acknowledge the support by the Australian Government
through the Australian Research Council’s Discovery Projects funding scheme
(project DP160100714). The views expressed herein are those of the authors and
are not necessarily those of the Australian Research Council.

## References

  * Peskin [1972] C. S. Peskin,  Flow patterns around heart valves: A numerical method,  Journal of Computational Physics 10 (1972) 252–271.
  * Lai and Peskin [2000] M.-C. Lai, C. S. Peskin,  An Immersed Boundary Method with Formal Second-Order Accuracy and Reduced Numerical Viscosity,  Journal of Computational Physics 160 (2000) 705–719. doi:10.1006/jcph.2000.6483.
  * Choi et al. [2007] J.-I. Choi, R. C. Oberoi, J. R. Edwards, J. A. Rosati,  An immersed boundary method for complex incompressible flows,  Journal of Computational Physics 224 (2007) 757–784.
  * Colonius and Taira [2008] T. Colonius, K. Taira,  A fast immersed boundary method using a nullspace approach and multi-domain far-field boundary conditions,  Computer Methods in Applied Mechanics and Engineering 197 (2008) 2131–2146. Immersed Boundary Method and Its Extensions.
  * Mittal and Iaccarino [2005] R. Mittal, G. Iaccarino,  Immersed boundary methods,  Annual Review of Fluid Mechanics 37 (2005) 239–261.
  * Yokoi et al. [2005] K. Yokoi, F. Xiao, H. Liu, K. Fukasaku,  Three-dimensional numerical simulation of flows with complex geometries in a regular cartesian grid and its application to blood flow in cerebral artery with multiple aneurysms,  Journal of Computational Physics 202 (2005) 1–19.
  * Delorme et al. [2013] Y. Delorme, K. Anupindi, A. Kerlo, D. Shetty, M. Rodefeld, J. Chen, S. Frankel,  Large eddy simulation of powered fontan hemodynamics,  Journal of Biomechanics 46 (2013) 408–422. Special Issue: Biofluid Mechanics.
  * de Zélicourt et al. [2009] D. de Zélicourt, L. Ge, C. Wang, F. Sotiropoulos, A. Gilmanov, A. Yoganathan,  Flow simulations in arbitrarily complex cardiovascular anatomies – an unstructured cartesian grid approach,  Computers & Fluids 38 (2009) 1749–1762.
  * Anupindi et al. [2013] K. Anupindi, Y. Delorme, D. A. Shetty, S. H. Frankel,  A novel multiblock immersed boundary method for large eddy simulation of complex arterial hemodynamics,  Journal of Computational Physics 254 (2013) 200–218. doi:10.1016/j.jcp.2013.07.033.
  * Löhner et al. [2004] R. Löhner, J. D. Baum, E. Mestreau, D. Sharov, C. Charman, D. Pelessone,  Adaptive embedded unstructured grid methods,  International Journal for Numerical Methods in Engineering 60 (2004) 641–660. doi:10.1002/nme.978.
  * Löhner et al. [2008] R. Löhner, J. R. Cebral, F. E. Camelli, S. Appanaboyina, J. D. Baum, E. L. Mestreau, O. A. Soto,  Adaptive embedded and immersed unstructured grid techniques,  Computer Methods in Applied Mechanics and Engineering 197 (2008) 2173–2197. doi:10.1016/j.cma.2007.09.010.
  * Sotiropoulos and Yang [2014] F. Sotiropoulos, X. Yang,  Immersed boundary methods for simulating fluid–structure interaction,  Progress in Aerospace Sciences 65 (2014) 1–21. doi:10.1016/j.paerosci.2013.09.003.
  * Glowinski et al. [2001] R. Glowinski, T. Pan, T. Hesla, D. Joseph, J. Périaux,  A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow,  Journal of Computational Physics 169 (2001) 363–426. doi:10.1006/jcph.2000.6542.
  * Glowinski and Kuznetsov [2007] R. Glowinski, Y. Kuznetsov,  Distributed lagrange multipliers based on fictitious domain method for second order elliptic problems,  Computer Methods in Applied Mechanics and Engineering 196 (2007) 1498–1506. doi:10.1016/j.cma.2006.05.013.
  * Zhang et al. [2004] L. Zhang, A. Gerstenberger, X. Wang, W. K. Liu,  Immersed finite element method,  Computer Methods in Applied Mechanics and Engineering 193 (2004) 2051–2067. doi:10.1016/j.cma.2003.12.044.
  * Liu et al. [2006] W. K. Liu, Y. Liu, D. Farrell, L. Zhang, X. S. Wang, Y. Fukui, N. Patankar, Y. Zhang, C. Bajaj, J. Lee, J. Hong, X. Chen, H. Hsu,  Immersed finite element method and its applications to biological systems,  Computer Methods in Applied Mechanics and Engineering 195 (2006) 1722–1749. doi:10.1016/j.cma.2005.05.049.
  * Wang and Liu [2004] X. Wang, W. K. Liu,  Extended immersed boundary method using fem and rkpm,  Computer Methods in Applied Mechanics and Engineering 193 (2004) 1305–1321.
  * Düster et al. [2008] A. Düster, J. Parvizian, Z. Yang, E. Rank,  The finite cell method for three-dimensional problems of solid mechanics,  Computer Methods in Applied Mechanics and Engineering 197 (2008) 3768–3782. doi:10.1016/j.cma.2008.02.036.
  * Ramière et al. [2007] I. Ramière, P. Angot, M. Belliard,  A fictitious domain approach with spread interface for elliptic problems with general boundary conditions,  Computer Methods in Applied Mechanics and Engineering 196 (2007) 766–781. doi:10.1016/j.cma.2006.05.012.
  * Neittaanmäki and Tiba [1995] P. Neittaanmäki, D. Tiba,  An embedding of domains approach in free boundary problems and optimal design,  SIAM Journal on Control and Optimization 33 (1995) 1587–1602. doi:10.1137/S0363012992231124.
  * Belytschko et al. [2003] T. Belytschko, C. Parimi, N. Moës, N. Sukumar, S. Usui,  Structured extended finite element methods for solids defined by implicit surfaces,  International Journal for Numerical Methods in Engineering 56 (2003) 609–635. doi:10.1002/nme.686.
  * Haslinger and Renard [2009] J. Haslinger, Y. Renard,  A new fictitious domain approach inspired by the extended finite element method,  SIAM Journal on Numerical Analysis 47 (2009) 1474–1499. doi:10.1137/070704435.
  * Legrain et al. [2012] G. Legrain, N. Chevaugeon, K. Dréau,  High order x-fem and levelsets for complex microstructures: Uncoupling geometry and approximation,  Computer Methods in Applied Mechanics and Engineering 241-244 (2012) 172–189. doi:10.1016/j.cma.2012.06.001.
  * Casquero et al. [2015] H. Casquero, C. Bona-Casas, H. Gomez,  A nurbs-based immersed methodology for fluid–structure interaction,  Computer Methods in Applied Mechanics and Engineering 284 (2015) 943–970. doi:https://doi.org/10.1016/j.cma.2014.10.055.
  * Casquero et al. [2021] H. Casquero, C. Bona-Casas, D. Toshniwal, T. J. Hughes, H. Gomez, Y. J. Zhang,  The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics,  Journal of Computational Physics 425 (2021) 109872. doi:https://doi.org/10.1016/j.jcp.2020.109872.
  * Xu et al. [2016] F. Xu, D. Schillinger, D. Kamensky, V. Varduhn, C. Wang, M.-C. Hsu,  The tetrahedral finite cell method for fluids: Immersogeometric analysis of turbulent flow around complex geometries,  Computers & Fluids 141 (2016) 135–154. doi:https://doi.org/10.1016/j.compfluid.2015.08.027.
  * Kamensky et al. [2015] D. Kamensky, M.-C. Hsu, D. Schillinger, J. A. Evans, A. Aggarwal, Y. Bazilevs, M. S. Sacks, T. J. Hughes,  An immersogeometric variational framework for fluid–structure interaction: Application to bioprosthetic heart valves,  Computer Methods in Applied Mechanics and Engineering 284 (2015) 1005–1053. doi:https://doi.org/10.1016/j.cma.2014.10.040.
  * Kamensky et al. [2017] D. Kamensky, M.-C. Hsu, Y. Yu, J. A. Evans, M. S. Sacks, T. J. Hughes,  Immersogeometric cardiovascular fluid–structure interaction analysis with divergence-conforming b-splines,  Computer Methods in Applied Mechanics and Engineering 314 (2017) 408–472. doi:https://doi.org/10.1016/j.cma.2016.07.028.
  * Burman and Hansbo [2010] E. Burman, P. Hansbo,  Fictitious domain finite element methods using cut elements: I. a stabilized lagrange multiplier method,  Computer Methods in Applied Mechanics and Engineering 199 (2010) 2680–2686. doi:10.1016/j.cma.2010.05.011.
  * Gerstenberger and Wall [2010] A. Gerstenberger, W. A. Wall,  An embedded dirichlet formulation for 3d continua,  International Journal for Numerical Methods in Engineering 82 (2010) 537–563. doi:10.1002/nme.2755.
  * Baiges et al. [2012] J. Baiges, R. Codina, F. Henke, S. Shahmiri, W. A. Wall,  A symmetric method for weakly imposing dirichlet boundary conditions in embedded finite element meshes,  International Journal for Numerical Methods in Engineering 90 (2012) 636–658. doi:10.1002/nme.3339.
  * Hautefeuille et al. [2012] M. Hautefeuille, C. Annavarapu, J. E. Dolbow,  Robust imposition of dirichlet boundary conditions on embedded surfaces,  International Journal for Numerical Methods in Engineering 90 (2012) 40–64. doi:10.1002/nme.3306.
  * Bastian and Engwer [2009] P. Bastian, C. Engwer,  An unfitted finite element method using discontinuous galerkin,  International Journal for Numerical Methods in Engineering 79 (2009) 1557–1576. doi:10.1002/nme.2631.
  * Rangarajan et al. [2009] R. Rangarajan, A. Lew, G. C. Buscaglia,  A discontinuous-galerkin-based immersed boundary method with non-homogeneous boundary conditions and its application to elasticity,  Computer Methods in Applied Mechanics and Engineering 198 (2009) 1513–1534. doi:10.1016/j.cma.2009.01.018.
  * Bourantas et al. [2021] G. Bourantas, D. Lampropoulos, B. Zwick, V. Loukopoulos, A. Wittek, K. Miller,  Immersed boundary finite element method for blood flow simulation,  Computers & Fluids 230 (2021) 105162.
  * unk [2019] FEniCS Project, 2019. URL: https://fenicsproject.org.
  * Richardson et al. [2019] C. N. Richardson, N. Sime, G. N. Wells,  Scalable computation of thermomechanical turbomachinery problems,  Finite Elements in Analysis and Design 155 (2019) 32–42.
  * Wu and Shu [2009] J. Wu, C. Shu,  Implicit velocity correction-based immersed boundary-lattice Boltzmann method and its applications,  Journal of Computational Physics 228 (2009) 1963–1979. doi:10.1016/j.jcp.2008.11.019.
  * Goda [1979] K. Goda,  A multistep technique with implicit difference schemes for calculating two- or three-dimensional cavity flows,  Journal of Computational Physics 30 (1979) 76–95. doi:10.1016/0021-9991(79)90088-3.
  * Chorin [1968] A. J. Chorin,  Numerical solution of the Navier-Stokes equations,  Mathematics of Computation 22 (1968) 745–762. doi:10.1090/S0025-5718-1968-0242392-2.
  * Témam [1969] R. Témam,  Sur l’approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires (I),  Archive for Rational Mechanics and Analysis 32 (1969) 135–153. doi:10.1007/BF00247678.
  * Khan et al. [2019] M. O. Khan, K. Valen-Sendstad, D. A. Steinman,  Direct numerical simulation of laminar-turbulent transition in a non-axisymmetric stenosis model for newtonian vs. shear-thinning non-newtonian rheologies,  Flow, Turbulence and Combustion 102 (2019) 43–72. URL: https://doi.org/10.1007/s10494-018-9905-7. doi:10.1007/s10494-018-9905-7.
  * Bathe [2001] K.-J. Bathe,  The inf–sup condition and its evaluation for mixed finite element methods,  Computers & Structures 79 (2001) 243–252. doi:10.1016/S0045-7949(00)00123-1".
  * Van De Vosse et al. [1989] F. N. Van De Vosse, A. A. Van Steenhoven, A. Segal, J. D. Janssen,  A finite element analysis of the steady laminar entrance flow in a 90° curved tube,  International Journal for Numerical Methods in Fluids 9 (1989) 275–287. doi:10.1002/fld.1650090304.