{'fraction_non_alphanumeric': 0.039379900213827514, 'fraction_numerical': 0.01701710620099786, 'mean_word_length': 4.5291871921182265, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Immersed boundary finite element method for blood flow simulation

We present an efficient and accurate immersed boundary (IB) finite element (FE) solver for numerically solving incompressible Navier--Stokes equations. Particular emphasis is given to internal flows with complex geometries (blood flow in the vasculature system). IB methods are computationally costly for internal flows, mainly due to the large percentage of grid points that lie outside the flow domain. In this study, we apply a local refinement strategy, along with a domain reduction approach in order to reduce the grid that covers the flow domain and increase the percentage of the grid nodes that fall inside the flow domain. The proposed method utilizes an efficient and accurate FE solver with the incremental pressure correction scheme (IPCS), along with the boundary condition enforced IB method to numerically solve the transient, incompressible Navier--Stokes flow equations. We verify the accuracy of the numerical method using the analytical solution for Poiseuille flow in a cylinder. 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 which would otherwise be extremely difficult to solve by most IB methods. Our method offers high accuracy, as demonstrated by the verification examples, and high efficiency, as demonstrated through the solution of blood flow within complex geometry on an off-the-shelf laptop computer.


Introduction
Immersed boundary methods have been used to simulate flow past rigid and moving/deforming bodies with irregular/complex shapes using Cartesian grids. 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 discrete set of points (following the literature [1], referred to them as Lagrangian points) for the imposition of prescribed velocity boundary condition. 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 a new mesh at each time step, but only the updated position of the points describing the immersed boundary. However, high computational cost of most currently used IB methods in application to internal flows remains a challenge. [2][3][4].

State-of-the art of the existing IB methods
In IB methods, flow equations are discretized on a Cartesian grid that does not conform with the boundaries of the immersed object, and the boundary conditions are imposed indirectly through modifications of the governing Navier-Stokes (N-S) equations. The modification applies as a forcing function in the governing equations that incorporates the presence of the immersed boundary. According to a forcing function f (x, t) formulation, the existing IB methods [5] are divided in two groups: continuous (or feedback) forcing and discrete (or direct) forcing methods. These two approaches are often called diffuse and sharp interface methods, respectively.

Continuous forcing approach
In the continuous forcing approach, the forcing function f (x, t) is included into the momentum equation (N-S equations). This approach was introduced by Peskin [1,6] to model blood flow through a beating heart. Peskin's method is a mixed Eulerian-Lagrangian finite difference method for computing the flow interaction with a flexible immersed boundary. The fluid flow is governed by the incompressible N-S equations, numerically solved on a stationary Cartesian grid, while the immersed boundary is represented by a set of massless elastic fibres. The locations of the fibres is tracked in a Lagrangian fashion by a collection of massless points that move with the local fluid velocity.
Peskin's method was later applied to rigid bodies by Goldstein et al. [7]. In this method, often called virtual boundary method, the immersed boundary is treated as a fictitious boundary embedded into the fluid. The embedded boundary applies force to the fluid, such that the fluid will be at rest on the surface (no-slip condition). The force F (s, t) on the boundary is computed such that the fluid velocity u (x, t) satisfies the no-slip condition on the boundary. The body force, is not known a priori and is computed through the velocity on the immersed boundary in a feedback forcing manner.
Saiki and Biringen [8] extended the feedback forcing approach, and eliminated the spurious oscillations caused by the applied feedback forcing term at the boundary. The force F (s, t) on the boundary has been modified, and a more accurate interpolation of the fluid velocity at the boundary points has been developed. A similar approach to the virtual boundary method has been introduced by Lai and Peskin [9]. In this approach, to simulate the flow around a rigid boundary, the boundary is not fixed and undergoes small movement. The immersed boundary points X are connected to fixed equilibrium points X eq using a very stiff spring with large stiffness constant. In case the boundary points move away from the desired location, the spring force will pull these boundary points back.
Khadra et al. [10] developed the penalty continuous forcing approach. The main idea, is that the entire flow occurs in a porous medium, governed by the Navier-Stokes-Brinkman equations.
The equations contain an additional term of volume drag, called Darcy drag, which accounts for the action of the porous medium on the flow. Su et al. [11] proposed a new implicit force formulation on the Lagrangian marker to ensure exact imposition of no-slip boundary condition at the immersed boundary. The immersed solid boundary is represented by discrete Lagrangian markers, which apply forces to the Eulerian fluid domain. The Lagrangian markers and the fluid variables on the fixed Eulerian grid exchange information through a simple discrete delta function.
Glowinski et al. [12] proposed the distributed Lagrange multiplier method (DLM). The method uses a finite element (FE) method (projection or predictor-corrector method) and introduces Lagrange multipliers (i.e. body force) on the immersed boundary to satisfy the no-slip condition. Lee and LeVeque [13] developed the immersed interface method (IIM), a method that was initially developed for elastic membranes. In the IIM, the boundary force/force strength F (s, t) is decomposed into tangential and normal components. The interface is explicitly tracked in a Lagrangian manner. The tangential component of the force is included in the momentum equation as an explicit term and the explicit normal boundary force is implemented into the governing equations in terms of a pressure jump condition across the interface [14].

Discrete forcing approach
In the discrete forcing approach, the flow equations are discretized on a Cartesian grid neglecting the immersed boundary. Instead, the discretization at the Cartesian grid points near the immersed boundary is adjusted to account for the presence of the immersed boundary. This discrete approach is is better suited for higher Reynolds numbers, since the velocity boundary conditions at the immersed boundary are imposed without introducing any forcing term.
Mohd-Yusof [15] developed a spectral method, which uses a forcing term that is numerically computed using the difference between the interpolated and the prescribed velocity on the immersed boundary points. This way, the errors between the calculated and the prescribed velocity on the immersed boundary are compensated by the forcing term. Fadlun et al. [16] extended the discrete-time forcing approach introduced by Mohd-Yusof [15] to a three-dimensional finite difference method on a standard marker-and-cell (MAC) staggered grid and showed that the approach was more efficient than feedback-forcing. Balaras [17] proposed a better reconstruction scheme, based on the method of Mohd-Yusof [15] and Fadlun et al. [16], which applies spreading/interpolation along the line normal to the immersed boundary. The algorithm eliminates the ambiguities associated with interpolation along the Cartesian grid lines, but it is limited only to flows with immersed boundaries aligned with one coordinate direction. Gilmanov et al. [18] developed a new reconstruction scheme, which applies to complex, three-dimensional immersed boundaries. The proposed scheme maintains a sharp fluid/immersed boundary interface by discretizing the immersed boundary using an triangular mesh. The numerical solution on the Cartesian grid points close to the immersed boundary points (Lagrangian points) is reconstructed via linear interpolation along the boundary surface normal vector.
Zhang and Zheng [19] developed an improved version of Mohd-Yusof's method by developing a bilinear interpolation/extrapolation function to interpolate force. Herein, both tangential and the normal velocity components (as opposed to Mohd-Yusof's method) are used to interpolate force. Therefore, the accuracy in the vicinity of the immersed boundary is enhanced without using higherorder schemes or using body-fitted grids near the immersed boundary. Choi et al. [20] developed a finite-volume immersed boundary method valid to all Reynolds numbers and is suitable for implementation on arbitrary grid topologies. They introduced the concept of tangency correction by decomposing the velocity into tangential and normal components along the outward normal direction to the immersed surface. The tangential velocity component is expressed as a power-law function (which is rather arbitrary) of the wall normal distance.

Contributions of this study
In the present study we combine the boundary condition-enforced immersed boundary (BCE-IB) method [21,22] with the incremental pressure correction scheme (IPCS) [23] using the finite element (FE) method. The main advantage of the BCE-IB method is that it satisfies accurately both the 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 satises the prescribed velocity boundary conditions. The incremental pressure correction scheme (IPCS) is a modified version of projection method, which provides improved accuracy at little extra computational cost. We implement the method using FEniCS [24] (an open-source software package that solves partial differential equations) to numerically solve the incompressible Navier-Stokes equations. 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, since the advantageous features of IB methods, that have made them popular for external flow simulations, do not entirely apply for internal flows. Users of IB methods are aware of this problem and have come up with various strategies to circumvent this liability [2][3][4].
The remainder of the paper is organized as follows. In Section 2, we present the numerical formulation for the proposed immersed boundary (IB) method. In Section 3, we solve benchmark problems to verify the accuracy of the proposed IB method. Section 4 is concerned with the numerical examples that demonstrate the accuracy and applicability of the proposed IB method. Section 5 contains conclusions.

Governing equations
We consider the incompressible viscous flow in a spatial domain Ω, which contains an immersed boundary in the form of a closed surface ∂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 u and pressure p) formulation -accounting for the immersed boundary are written as: subject to the no-slip boundary condition on ∂S with U B being the prescribed velocity of the immersed boundary, and v the kinematic viscosity of the fluid. The term ε(u) is the strain-rate tensor defined as: The forcing term 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 f (x, t) is the local body force density at the fluid nodes (Eulerian nodes), it is distributed from the surface force density F (s, t) (s is a parametrization of the immersed boundary surface) at the immersed boundary points (Lagrangian points -we use traditional terminology used in [22]). The forcing term is expressed as: where x and X (s, t) denote the Eulerian nodes and Lagrangian point coordinates that discretize the fluid domain and the immersed boundary, respectively. Eulerian nodes and Lagrangian points interact through Dirac delta function δ (x − X (s, t)). Therefore, the velocity at the immersed boundary points is calculated from the velocity at the Eulerian nodes using the Dirac delta function

Incremental Pressure Correction Scheme (IPCS)
The Incremental Pressure Correction Scheme (IPCS) [23] is an operator splitting method, which transforms the nonlinear Navier-Stokes (N-S) equations into algebraic equations by coupling the pressure and velocity field values (N-S equations are difficult to solve due to their inherent non-linearity). IPCS is a modified version of the fractional step method proposed by Chorin [25] and Témam [26]. 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 in Section 2.3 we explain in detail how IB forces were incorporated.
The IPCS scheme involves three steps. In the first step, we compute a tentative velocity u * by advancing the linear momentum Eq.(1) in time using the backward Euler difference scheme, u = u n+1 −u n ∆t . We linearize the nonlinear terms using a semi-implicit method, such that the advection term u · ∇u becomes u n · ∇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 = v (u n ), which simplifies the N-S equations that are now written as a linearized set of equations, known as the Oseen equations: ∇ · u n+1 = 0.
In Eq. (7), p n+1 is still unknown. Therefore, we will compute the tentative velocity u * (u * ≈ 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 u * : Equation (9) is numerically solved to compute the tentative velocity u * without using the incompressibility constraint Eq.(2).
In the second step, we use the tentative velocity u * to compute the updated velocity u n+1 , which is divergence-free (and should fulfill the incompressibility constraint). We define a function for the velocity correction as u c = u n+1 − u * , and after some algebra (subtracting Eq.(9) from Eq. (7)) and using that ∇ · u n+1 = 0 ⇒ ∇ · u c = −∇ · u * (incompressibility constraint) we obtain a new set of flow equations for u c : where Φ n+1 = p n+1 − p n . Equation (10) can be further simplified to u c + ∆t ∇Φ n+1 = 0 and the system of Eqs.(10)-(11) is reduced to an elliptic type (Poisson) partial differential equation for the pressure difference Φ n with boundary conditions being those applied to the original flow problem.

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: • Predictor step: where we numerically solve the N-S equations to compute the predicted velocity fieldũ (x) by disregarding the body force terms in Eq.
• Corrector step: where we account for the effect of body forces and update the predicted velocity fieldũ to the updated one u n+1 , which satisfies the boundary condition u (X (s In the predictor step, Eq.(13) is used to compute the predicted velocity fieldũ 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ũ, which in this case is the tentative velocity u * (Step 1 in the IPCS algorithmic procedure). The corrector step, involves the evaluation of the unknown body force f n+1 and the update of the predicted velocityũ to the physical one u n+1 (Steps 2 and 3 in the IPCS algorithmic procedure).
Accurate and efficient computation of the unknown body force f n+1 is an important advantage of the Boundary Condition Enforced Immersed Boundary (BCE-IB) method. In this method, the body force f n+1 is equivalent to a velocity correction, which applies implicitly such that the velocity U (X (s) , t) on the Lagrangian points, interpolated from the physical velocity u (x, t) computed on the Eulerian nodes, equals to the prescribed boundary velocity U B (Eq. (3)). The body force term f n+1 is computed using the equation derived by substituting Eq. (14) to Eq. (6). The force density f n+1 (x, t), is evaluated on the Eulerian nodes, and it is computed by distributing ('spreading') the boundary force F (X (s) , t) on the Lagrangian points through the Dirac delta function δ (x − X (s, t)) (see Eq. (5)). Eq. (15) is written as and U n+1 B is no longer correlated with f n+1 but instead with F n+1 . To compute the boundary force 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 . . , M and the flow domain by the Eulerian nodes x j = x j , y j , z j , j = 1, 2, . . . , 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 δ (x − X (s, t)) is approximated by a continuous kernel distribution D x i − X j given as with the kernel δ (r) proposed by Lai and Peskin [9] δ(r) = Using Eq. (6) for the force term Eq. (16) is written as where ∆S i is the area of the i th surface element (segment). Eq. (20) forms a well-defined system of equations for the variables F n+1 (20) can be written in matrix notation as with A F , F and B F defined as where i andũ x j , respectively. By solving he system of equations (Eq. (22)) using a direct solver, we obtain the unknown boundary force F n+1 i (i = 1, 2, . . . , M) at all Lagrangian points. Boundary forces F n+1 i are then substituted into Eq. (19) and Eq. (16) to calculate the body force f n+1 on the Eulerian nodes, and the corrected physical velocity u n+1 .

Algorithmic procedure
Our approach combines the IPCS (section 2.2) with the BCE-IB method (section 2.3). 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. Calculate the tentative velocityũ using Eq.(13) 2. Compute the matrix A F using Eq. (22) 3. Solve the system (Eq. (21)) to compute the boundary force F n+1 X n+1 4. Useũ to update the corrected velocity u (n+1) = ∆t f n+1 +ũ ρ using Eq. (14) 5. Solve the Poisson equation ∆t for Φ n (Eq.(12)) 6. Recompute the updated velocity u n+1 = u n + ∆t∇Φ n+1 and pressure p n+1 = Φ n+1 + p n . 7. Set physical/updated velocity u (n+1) as u (n) 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.

Surface area for Lagrangian points
In 3D flow cases, the surface area ∆S i that is assigned to each Lagrangian point is not as straightforward to define as in two dimensions. Immersed boundaries are represented by a set of Lagrangian points. Therefore, the number of these points and the way they are distributed over the immersed boundary directly affect the accuracy of the IB method.
In 2D, the distribution and number of Lagrangian points are always investigated through the ratio of ds h , with ds being spacing of the Lagrangian points and h the Eulerian nodal spacing. High ratio of ds h will cause fluid leakage (spreading from Lagrangian points to Eulerian nodes and vice-versa is not accurate), while low ratio will increase the computational cost (as rows in matrix A F in Eq.(21) will be linearly dependent and the matrix will become ill-conditioned and therefore difficult to invert). In the literature, there are many studies which report on selection of ds h ; Uhlmann [27] suggests that Lagrangian points spacing should be equal to the Eulerian mesh/grid resolution ds h = 1, Su et al. [11] and Kang and Hassan [28] used ds h = 1.5 (Lagrangian points spacing is larger than the Eulerian mesh/grid resolution), while Lima E Silva et al. [29] reported ds h ≤ 0.9 (Lagrangian points spacing is less than the Eulerian nodal resolution). In 3D, the surface area ∆S i appears as a key feature of the IB method. In this study, we propose and use an easy to apply method to compute ∆S i . The area of each surface element ∆s i should be close to (ah) 2 , with 0.1 ≤ a ≤ 2.5 (just like having in 2D 0.3 ≤ ds h ≤ 2.5, with h being the spacing of the Eulerian nodes).
In this study, we uniformly distribute the Lagrangian points to assign the same value ∆S to all Lagrangian points. To obtain a distribution of Lagrangian points with uniform surface element area, we generate surface mesh (the nodes are the Lagrangian points used in the IB method) with element edges having lengths with small standard deviation. The area element ∆S for each Lagrangian point is the computed mean value of the facets area.

Derived quantities-wall shear stress (WSS)
An important aspect of blood flow simulations is the numerical computation of flow parameters that are of clinical importance, such as blood vessel wall shear stress (WSS). While focusing on WSS, we present an accurate meshless point collocation to compute derived quantities from the velocity field obtained from blood flow simulations. Meshless methods (MMs) have been wellestablished for the efficient and accurate representation of 3D complex geometries as point clouds, and for the accurate computation of spatial derivatives. In our study, we use the Discretization Corrected Particle Strength Exchange (DC PSE) [30] meshless method. The DC PSE method, is one of the most accurate meshless point collocation methods to compute spatial derivatives of field functions defined on point clouds.
WSS is defined as the magnitude of the surface traction vector t s = t − (t ·n) ·n, wheren is the unit normal vector pointing outward, and t = 2µε is the shear stress, defined based on the fluid's viscosity µ and the rate of deformation tensor ε. Let u (x, t) be continuously differentiable velocity vector field describing the blood flow at a point in the domain of interest. The velocity field values are computed by interpolating the velocity field computed on the Cartesian flow domain using the IB method, on a body-fitted mesh of the immersed object (the mesh is used only for visualization of the flow field and derived quantities such as WSS).We consider a discretization of the domain of interest (immersed object) using a finite-element mesh (tetrahedral linear elements), and denote the corresponding set of nodes by N. The nodal subset N S with points on the surface of the tetrahedral mesh will be used for the calculation of WSS. Using indicial notation for the components of u, the deformation (strain-rate) tensor components are defined by ε i j = 1 2 u i, j + u j,i for i, j = {x, y, z}, where u i, j = ∂u i ∂x j . The terms ∂u i ∂x j of the deformation tensor ε are computed on the surface points N S using the DC PSE method as with Φ , j being the j th spatial derivative computed, and u i being a column vector of the i th velocity vector components.

Algorithm verification
We demonstrate the accuracy, efficiency and applicability of our method on internal flow cases. The numerical results obtained were compared against analytical solutions (where available) and experimental data. We consider two 3D benchmark flown problems, the Poiseuille flow in a straight tube and the angle tube bend flow problem, where an analytical solution (Poiseuille flow) and experimental data (tube bend) are available, respectively.

Poiseuille flow in a cylindrical tube
We demonstrate the accuracy of the proposed IB method considering the Poiseuille flow in a straight cylindrical tube. The tube has length L = 10R and radius R = 0.005 m, and is immersed in a domain (which is the spatial domain where the flow equations are solved) with dimensions −0.006 m ≤ y, z ≤ 0.006 m and 0 m ≤ x ≤ 0.05 m, as shown in Fig. 1. The immersed cylinder is described by a set of points on its surface (Lagrangian points). The box domain is discretized with a high quality tetrahedral mesh. The vertices of the tetrahedral elements (Eulerian nodes) form a uniform Cartesian grid. To generate the mesh we use the FEniCS build-in function BoxMesh [24].
The driving force of the flow is the pressure difference applied to the inlet and outlet. In our simulations, the density of the fluid was set to ρ = 1, 050 kg m −3 , and the dynamic viscosity to µ = 0.00345 N s/m 2 . At the inlet (red circle in Fig. 1a) and outlet (green circle in Fig. 1a), we apply pressure boundary conditions, while at the remaining walls (surface yellow Fig. 1a), we apply no-slip velocity boundary conditions (u = 0). The time step was set to dt = 2.5 × 10 −3 . The simulation ends when the normalized root mean square difference between two successive time instances for all three velocity components is less than 10 −8 (practically, the flow reaches a steady state). The exact solution for the velocity field components is u x = u y = 0 and u z (r) = U max 1 − r R 2 , with U max = − R 2 4µ dp dz being the maximum velocity. By setting the pressure difference to dp dz = 1 L (applying p inlet = 1 and p outlet = 0) the Reynolds number becomes Re = U max L µ ≈ 551, while for p inlet = 2 Reynolds number becomes Re = U max L µ ≈ 1, 102.  We use successively denser tetrahedral meshes to obtain a mesh independent numerical solution. The meshes are described in Table 1. Table 2 for different mesh resolution, for Re = 551 and Re = 1, 102. For Re = 551 the convergence rates for the L ∞ and L 2 error norms are 1.29 and 2.79, respectively. The convergence rates for L ∞ and L 2 error norms do not change, demonstrating that the relative error does not change. Figure 2 shows the L ∞ and L 2 error norms versus the mesh resolution of the Eulerian mesh. For both Re = 551 and Re = 1, 102, both L ∞ and L 2 error norms, the results suggest that the proposed scheme is accurate. Figure 3 shows the velocity streamlines for Re = 551 and Re = 1, 102. The streamlines were produced on a tetrahedral body-fitted mesh (we refer to as visualization mesh) that consist of 345,600 linear elements. 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).

3D angle tube bend
We further demonstrate the accuracy of the proposed method, considering the flow in a threedimensional 90 angle tube bend, as shown in Fig. 4. We use this flow example to demonstrate the ability of the proposed method to accurately capture flows in tortuous tubes where secondary flows can occur. For U-bend, flow data are available; Van De Vosse et al. [31] performed laser Doppler velocimetry experiments to obtain the center-plane axial velocities at a set of different  angles around the tube bend. To reduce disturbances in the flow field due to inflow boundary conditions, we extended the inlet and outlet by 1 mm. In our numerical simulations, we set the Reynolds number-computed based on the tube inner diameter D i and the mean inlet velocity U-to 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.
The U-bend geometry is embedded within a box domain with dimensions −0.0064 m ≤ x ≤ 0.032 m, 0 m ≤ y ≤ 0.0384 m and −0.0064 m ≤ z ≤ 0.0064 m, as shown in Fig. 5b. We discretize the flow domain using tetrahedral elements. The box domain is discretized with a high quality tetrahedral mesh (the vertices of the tetrahedral elements are the Eulerian nodes), generated using a uniform Cartesian grid (to generate the mesh we used the FEniCS build-in function BoxMesh).
The tetrahedral mesh is locally refined in the vicinity and within the interior of the immersed boundary, using the FEniCS built-in function Refine [24] which uses the algorithm by Plaza and Carey [32]. The local refinement (tetrahedral mesh resolution h close to the Lagrangian points) is based on the surface area ∆S of each Lagrangian point, such that h = a √ ∆S , with 0.1 ≤ a ≤ 2.5 (see Section 2.5). We uniformly distribute the Lagrangian points, and we use as ∆S the mean value of the surface area of the facets of the immersed boundary. We examine the distribution of the Lagrangian points through a histogram plot of the area of the facets of the immersed boundary. Figure 6a shows the histogram plot for the edge length of the immersed boundary facets, while Fig. 6b shows the histogram plot for the area of the immersed boundary facets for the finest Lagrangian point cloud (16,760 nodes and 31,462 facets). For this point cloud, the mean value for the area of the facets is ∆S mean = 4.2848058 × 10 −8 and ∆S std = 8.96531 × 10 −9 . The tetrahedral mesh resolution h in the vicinity and inside the immersed boundary that corresponds to the √ ∆S mean = 2.06997 × 10 −4 ≈ 2 × 10 −4 (see Table 3) is h = 2 × 10 4 , such that ratio of h √ ∆S mean ≈ 1. We use successively denser tetrahedral meshes to obtain a mesh independent numerical solution. The meshes used are listed in Table 3. We use a parabolic velocity profile (u z (r) =  Fig. 5a), while at the outlet (green circle in Fig. 5a) we apply zero pressure (p = 0). At the remaining boundaries (shown in yellow in Fig. 5a) of the flow domain, we set no-slip boundary conditions. The (rectangular) flow domain has rigid walls (yellow faces in Fig. 5a) with no-slip boundary conditions, an inlet (red circle in Fig. 5a), where the fluid enters, and an outlet (green circle in Fig. 5a), where the fluid exits. The numerical solutions obtained using the the three meshes described in Table 3, were projected/interpolated on a body-fitted mesh consisting of 808,984 linear tetrahedral elements and 143,406 nodes. Fig. 7a shows the axial velocity computed with the proposed method using successively denser meshes listed in Table 3. Fig. 7b 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 in excellent agreement    Table 3 with h being the mesh resolution close to the immersed boundary, and (b) axial velocity computed with the proposed IB method (dashed line) compared against the experimental data (red dots) at the tube outlet for the U-bend flow problem.
between the velocity predicted using the proposed IB method and experimental results is observed. The numerical results obtained using the linear (P 1 /P 1 ) elements are also in good agreement with the experimental data, highlighting the accuracy of the proposed method.

Numerical examples: blood flow simulations
In this section, we demonstrate the accuracy and robustness of the proposed scheme through blood flow simulations in complex vascular geometries. In the first flow problem, we simulate blood flow in coronary artery bifurcation (Fig. 8), while in the second we consider blood flow in the ascending aorta (Fig. 11). Finally, we consider the flow in the ascending and descending aorta (Fig. 14).
The flow problems considered demonstrate the applicability of the proposed method for internal flow simulations. Internal flow cases, are practically out of reach for all IB methods. The computational cost of having a tetrahedral mesh (uniform or locally refined) to discretize the flow domain is high, since only a small percentage of the total number of mesh elements, falls into the domain of interest [4]. The proposed IB method allows for cropping the tetrahedral mesh (see example 4.3) that covers the complex (vascular) geometry, and at the same time increases the percentage of elements that are located inside the flow domain of interest. In all flow cases considered, we use a locally refined mesh, in the vicinity and inside the Lagrangian points. We demonstrate the accuracy of the proposed IB method, by comparing the numerical findings with those computed using a body-fitted mesh FE flow solver. Local refinement is an automated procedure implemented using FEniCS built-in functions (Plaza algorithm [32]).

Flow in an coronary artery bifurcation
In the first example, we simulate the blood flow in a coronary artery (CA) bifurcation (shown in Fig. 8). The immersed boundary (coronary artery) is embedded into a box with dimensions −0.0052 m ≤ x ≤ 0.01 m, −0.006 m ≤ y ≤ 0.006 m and 0 m ≤ z ≤ 0.0296 m. The inlet and two outlets of the CA geometry were extruded to ensure that the flow is fully developed.
We discretize the flow domain using tetrahedral elements. The box domain is discretized with a high quality tetrahedral mesh, generated using a uniform Cartesian grid (to generate the mesh we used the FEniCS build-in function BoxMesh). The tetrahedral mesh is locally refined in the vicinity and inside the immersed boundary (Figs. 8c-d). We use successively denser meshes to obtain a mesh independent numerical solution. The meshes used are listed in Table 4. The numerical results on the two finest meshes (for 627,040 and 3,641,633) are indistinguishable. Mesh generation is a straightforward and automated. The locally refined mesh with 3,641,633 elements is generated in less than 60 s (using a laptop computer with 2.7 GHz i7 quad-core processor and 16 GB internal memory), using the refine function in FEniCS. We define the mesh resolution close to the immersed object and in the interior based on the mean value of the surface elements area (see Section 3.2). The majority (85%) of the elements are located inside the immersed boundary (coronary artery).
We set the total time for the simulation to T = 0.8 s (one cardiac cycle), and the time step to dt = 5 × 10 −4 s. At the inlet ((red area located at the bottom surface in Fig.(8a)) we apply the pulsatile velocity waveform shown in Fig. 9 (we use a parabolic velocity profile since the Womersley number is small), and zero pressure boundary conditions at the two outlets (green area on the top surface in Fig.(8a)). At the remaining walls (yellow surface in Fig.(8a)), we apply noslip boundary conditions. We employ the Newtonian model for blood flow, with dynamic viscosity of µ = 0.00345 Pa · s and density of ρ = 1, 056 kg m 3 . To demonstrate the accuracy of our IB method, we compare the results obtained using this method with those for the body-fitted mesh. The body-fitted mesh consists of 154,681 nodes and 860,818 linear tetrahedral elements (P 1 /P 1 ). Table 5, lists the normalized root mean square error error norms at different time instances. The results obtained using our IB method and body fitted mesh reported in Table 5 are for practical purposes indistinguishable as in the bioengineering applications. The differences of under 5% would be overwhelmed by the biomechanical properties uncertainties. The decisive advantage   Table 5: L 2 error norm of the immersed boundary method numerical results against those computed using the body-fitted FE method at different time instances for the coronary artery bifurcation flow example.  of our method is the easiness of patient-specific computational grid generation with practically no increase in the 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 coronary artery bifurcation conducted here (Figure 8), this overhead was negligible less than 1% of the total computation time. Fig. 10 shows the streamlines in the coronary artery bifurcation at time t = 0.14 (diastole) and 0.25 sec (peak systole).

Blood flow in the ascending aorta
In the second 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. 11. The immersed boundary (ascending aorta) is embedded within a box (Fig. 11a) with dimensions −0.037 m ≤ x ≤ 0.035 m, −0.08 m ≤ y ≤ 0.04 m and 0 m ≤ z ≤ 0.12 m. The box domain is discretized with a high quality tetrahedral mesh, generated using a uniform Cartesian grid (to generate the mesh we used the FEniCS build-in function BoxMesh). The tetrahedral mesh is locally refined in the vicinity and inside the immersed boundary as shown in Fig.11c-d. We use successively denser meshes, indicated in Table 6, to obtain a mesh independent numerical solution. The numerical results on the two finest meshes (for 1,410,846 and 2,557,306) are indistinguishable. We refine the tetrahedral mesh close to the Lagrangian points and in the interior of the immersed boundary. As described in Section 2.5, the mesh resolution close to the immersed boundary (Lagrangian points) is based on on the mean area of the facets on the immersed object. Mesh generation and local refinement on the rectangular domain is a straightforward and efficient procedure. It took less than 30 s to generate the finest grid (information about this grid in the third row of Table 6) using a laptop computer with i7 quad-core processor with 16 GB internal memory.
We set the total time for the simulation to T = 0.8 s (one cardiac cycle), and the time step dt = 5 × 10 −4 s. At the inlet (red area in Fig. 11a) we apply the pressure waveform shown in Fig. 12, zero pressure boundary conditions at the two outlets (green area in Fig. 11a), and no-slip boundary conditions at the remaining walls (yellow surface in Fig. 11a). We employ the Newtonian model for blood flow, with dynamic viscosity of µ = 0.00345 Pa s and density of ρ = 1, 056 kg/m 3 . We use the diameter D and the maximum velocity U m as the characteristic length and velocity to calculate the Reynolds number.   To compare the flow velocity predicted using our IB method with the results from body fitted mesh, we interpolate our IB solution on a body fitted mesh consisting of 92,786 nodes and 504,628 linear tetrahedral elements. The comparison done at the time instances (every 0.1 s interval) indicated that the normalized root mean square error (NRMSE) of less than 5% for the three velocity components. We conducted also qualitative comparison of the predicted streamlines, pressure contours, and wall shear stress (WSS) contours at the time instances of t = 0.2 and t = 0.6 ( Fig.  13). No visually distinguishable differences between the results obtained using our IB method and approach using body fitted mesh can be found in Fig. 13. (every 0.1 sec) the IB numerical solution with the numerical results obtained using the body-fitted mesh. The normalized root mean square error (NRMSE) L NRMSE for the three velocity components is less than 5%. Fig. 13 shows the streamlines (left), pressure (middle) and wall shear stress (WSS) magnitude (right) contours at different time instances t = 0.2 and t = 0.6. The results were obtained using the mesh with 1,410,846 tetrahedral elements.

Blood flow in the ascending and descending aorta
As the final example, we consider blood flow in the ascending and descending aorta. This example differs from the previous one (blood flow in the ascending aorta) by the number of outlets considered (the brachiocephalic, left common carotid, and left subclavian arteries were removed).
This flow example is more challenging because the inlet of the immersed object is located in the interior of the box flow domain and not on the boundary (Fig. 14a). Therefore, we crop the rectangular domain such that the inlet of the immersed object (red circle in Fig. 14a) is on the boundary of the flow domain (yellow color in Fig. 14a) and the boundary conditions can be easily applied.   The immersed boundary (aorta) is embedded into a box domain (Fig. 14a) with dimensions −0.042 m ≤ x ≤ 0.0252 m, −0.0806 m ≤ y ≤ 0.025 m and −0.14 m ≤ z ≤ 0.0616 m. The box domain is discretized with a high quality tetrahedral mesh, generated using a uniform Cartesian grid (to generate the mesh we used the FEniCS build-in function BoxMesh). The tetrahedral mesh is cropped and elements are removed (the cropped domain is shown in Fig. 14b). Finally, the cropped tetrahedral mesh is locally refined in the vicinity and inside the immersed boundary (Fig. 14d).
We use successively denser meshes to obtain a mesh independent numerical solution. A mesh of 417,227 nodes and 2,409,460 linear tetrahedral elements (approx. 85% of the elements are located inside the immersed boundary) ensures a mesh independent numerical solution. We apply a uniform velocity profile U (plug flow) at the inlet (in red color Fig. 14a). We apply zero pressure boundary conditions at the outlet, and no-slip boundary conditions in the remaining boundaries. We employ the Newtonian model for blood flow, with dynamic viscosity of µ = 0.00345 Pa · s and density of ρ = 1, 056 kg/m 3 . The diameter D at the inlet, and the velocity at the inlet U, are used as the characteristic length and velocity to define the Reynolds number. For the inlet velocity value used in our simulations, the Reynolds number is Re = 310 for the present geometry. Fig. 15 shows the streamlines (left), pressure (middle) and wall shear stress (WSS) magnitude (right) contours at different time instances t = 0.2 and t = 0.6.

Conclusions
In this work, we present an immersed boundary method for internal flows, showing blood flow application. The proposed scheme combines the finite element (FE) method, which is used to numerically solve the Navier-Stokes equations, and the boundary condition enforced immersed boundary (BCE-IB) method to account for complex geometries.
The proposed method increases the computational efficiency of the IB method for internal flows, by reducing the number of unused mesh points through a sophisticated and efficient method for mesh refinement in the vicinity and inside the immersed boundary. This method facilitated generation of a computational mesh consisting of over three and half million elements (computation of flow in an aortic bifurcation see Fig. 8) in less than 30 seconds on an off-the-shelf laptop with quad-core i7 processor and 16 GB of internal memory. This very promising results make it possible to regard our adaptive mesh refinement method as compatible with the time constraints of clinical workflow (no firm definition here by the times of on order of under 1 minute have been reported in the literature [33]). Furthermore, we introduced an efficient method to reduce (crop) the flow domain, in flow cases where the immersed boundary inlet(s) and outlet(s) do not coincide with the boundaries of the flow domain (tetrahedral mesh).
The proposed method utilizes an efficient and accurate finite element solver based on incremental pressure correction scheme (IPCS). All simulations in this study were conducted using an off-the-shelf laptop with quad-core i7 processor and 16 GB of internal memory. The accuracy of the proposed scheme has been successfully verified through comparison with an analytical solution for the Poiseuille flow in a cylinder with 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 has been demonstrated through the solution of flow examples with complex geometry.