3D Fluid-Structure Interaction Simulation of Aortic Valves Using a Unified Continuum ALE FEM Model

Due to advances in medical imaging, computational fluid dynamics algorithms and high performance computing, computer simulation is developing into an important tool for understanding the relationship between cardiovascular diseases and intraventricular blood flow. The field of cardiac flow simulation is challenging and highly interdisciplinary. We apply a computational framework for automated solutions of partial differential equations using Finite Element Methods where any mathematical description directly can be translated to code. This allows us to develop a cardiac model where specific properties of the heart such as fluid-structure interaction of the aortic valve can be added in a modular way without extensive efforts. In previous work, we simulated the blood flow in the left ventricle of the heart. In this paper, we extend this model by placing prototypes of both a native and a mechanical aortic valve in the outflow region of the left ventricle. Numerical simulation of the blood flow in the vicinity of the valve offers the possibility to improve the treatment of aortic valve diseases as aortic stenosis (narrowing of the valve opening) or regurgitation (leaking) and to optimize the design of prosthetic heart valves in a controlled and specific way. The fluid-structure interaction and contact problem are formulated in a unified continuum model using the conservation laws for mass and momentum and a phase function. The discretization is based on an Arbitrary Lagrangian-Eulerian space-time finite element method with streamline diffusion stabilization, and it is implemented in the open source software Unicorn which shows near optimal scaling up to thousands of cores. Computational results are presented to demonstrate the capability of our framework.


INTRODUCTION
The World Health Organization (WHO, 2014) has identified cardiovascular disease as the major cause for death in the world. Therefore, developing new ways to support early diagnosis of cardiac dysfunction is of vital importance. In vivo and in vitro studies offer valuable information on the relationship between the blood flow (hemodynamics) and cardiac disease, and advances in computational fluid dynamics (CFD) and high performance computing (HPC) enable the usage of computer simulation as an important tool to further enhance our understanding of this relationship.
The field of cardiac modeling is extensive, and highly interdisciplinary. It is therefore important to be clear on what the research is aiming for. Our goal is to develop a framework for simulating the intraventricular blood flow, where specific properties such as fluid-structure interaction (FSI) of the aortic valve can be implemented in a modular way without extensive efforts. In Spühler et al. (2015) we focus on the aspect of fluid mechanics, and present a computational model of the blood flow in the left ventricle (LV) of the heart. The movement of the wall is based on ultrasound measurements and an Arbitrary Lagrangian-Eulerian (ALE) space-time finite element method is used to simulate the blood flow by solving the incompressible Navier-Stokes equations. The opening and closing of the mitral and aortic valves are modeled by time-dependent velocity and pressure boundary conditions. In this paper, we present an extension of this work by embedding different geometrical models of aortic valves in the LV and the aorta. Prototypes of a biological valve and bileaflet mechanical heart valve (BMHV) are modeled. While surgical treatments of valvular diseases are firmly established, many decisive factors for the performance of the implant are not fully understood yet. Numerical simulations provide an important insight to the interaction between the blood flow and the leaflets which can be applied to optimize the design of BMVHs or improve technologies as transcatheter aortic valve replacement (Wu et al., 2016). The fluid-structure interaction problem is described by a unified continuum model, using the conservation laws for mass and momentum and a phase function, which is a novel approach for simulating valve motions. The Navier-Stokes equations are solved by an ALE space-time finite element method with streamline diffusion stabilization implemented in Unicorn (Hoffman et al., 2012), which is part of the open source software framework FEniCS-HPC (Jansson, 2013). This paper is structured as follows. In section 2 we describe the different components and functions of an anatomical aortic valve. section 3 explains the mathematical equations and the numerical method. In section 4, we specify the mechanical and biological aortic valve model we use in our simulations. The numerical results are presented in section 5 and we conclude our paper in section 6 by summarizing our findings and discuss possible steps of future work.

MODELING THE AORTIC VALVE
The left ventricle possesses a mitral and an aortic valve, each of them consisting of two and three leaflets respectively. The valves ensure unidirectional flow and prevent the blood to flow back. The opening and closing of the valves are mainly controlled by the pressure gradient between the ventricle and the adjacent chamber. One edge of the biological leaflet is completely attached to the inner wall of the heart. The free edge of the mitral valve is connected to the papillary muscles by the chordae tendineae. The aortic leaflets do not have such fibrous tissue connections and they open and close passively due to the flow.
The nomenclature of the different components of the aortic root can vary remarkably as revealed by Sievers et al. (2012). We apply the definitions proposed in Sievers et al. (2012), as indicated in Figure 1. The aortic root is situated between the left ventricle and the ascending aorta, and is bordered by the annulus and the sinotubular junction. The three bulges just above the annulus are referred as sinus of Valsalva. The aortic valve contains three leaflets which are attached to the aorta wall. The point of contact where two leaflets meet at the root wall is called commissure and the surface of contact at the free edge is known as coaptation.

MATHEMATICAL MODEL AND NUMERICAL METHOD
In order to put our approach in context, we review different models for simulating the FSI of the blood flow around the aortic valve. Usually, the structure model is formulated in the Lagrangian coordinate system whereas fluid flow is described in the Eulerian coordinate system. At the common interface of the two models, the following kinematic and dynamic constraints have to be satisfied by the velocity u and the stress τ : u f = u s (kinematic constraint, continuity of the velocity), (1) τ s · n = τ f · n (dynamic constraint, continuity of the normal stresses). (2) The subscript indicates whether the variable is defined in the solid (s) or in the fluid (f ) part respectively, while n is a unit vector normal to the interface. We denote vectors and matrices with bold letters. FSI simulations can roughly be categorized as moving or fixed mesh methods and partitioned or monolithic approach as presented in Borazjani et al. (2008).

Discretization of the Coupled Problem
For fixed mesh methods, the fluid and structure domains are discretized in a non-boundary conforming matter. Since the structure is spatially disconnected from the fixed background mesh, it is crucial to efficiently trace and move the interface between the solid and the fluid domain. The interface can be discretized with a set of markers and tracked by a Lagrangian method (front tracking) or represented by contours or level sets of a scalar function (front capturing). Fixed mesh methods were pioneered by Peskin and McQueen (Peskin, 1972;McQueen and Peskin, 2000) introducing the concept of immersed boundary methods, where body forces are imposed on the fluid domain to account for the interaction between the fluid and the structure. Large structural deformations are manageable, but the solution at the interface can be diffuse. This disadvantage can be lessened by e.g., increasing the mesh resolution in the vicinity of the immersed boundary as done by Griffith (2012), or by treating the boundary as a sharp interface as in e.g., Borazjani et al. (2008); Udaykumar et al. (1999); Mittal and Iaccarino (2005); Gilmanov and Sotiropoulos (2005), and Xia et al. (2009). Fictitious domain methods is a another class of fixed mesh methods where Lagrange multipliers account for the kinematics constraints between the fluid and solid domain, see e.g., Glowinski et al. (1999);van Loon et al. (2005); De Hart et al. (2003), and Astorino et al. (2009).
In moving mesh methods, the computational mesh conforms to the deformation of the solid domain, and is typically represented by an Arbitrary Lagrangian Eulerian formulation. The strength of the moving grid methods can be found in its accuracy and clearly defined coupling condition, as the mesh is aligned with the fluid-structure interface. A good smoothing algorithm or local remeshing is needed to keep the quality of the computational mesh.
Reviewing the literature of aortic valve simulations, we came across the following work which apply an ALE approach: Bolger et al. (2007) and Penrose and Staples (2002) simulate the flow past a geometrically reduced mechanical valve prosthesis taking advantage of its symmetrical form; in Dumont et al. (2007) two commercially available bileaflet mechanical heart valves are compared regarding hemodynamics and thrombogenic performance; Guivier-Curien et al. (2009) employs particle image velocity measurements to quantitatively and qualitatively compare experiments and numerical simulations; the FSI model of Choi and Kim (2009) provides detailed flow information and leaflet behavior of a BMHV; Morsi et al. (2007) analyzes the fluid dynamics of a trileaflet heart valve but only for the initial opening phase.

Coupling Strategies
Depending on whether the structure and fluid problems are solved simultaneously or separately, the FSI solver can be classified as monolithic or partitioned. The FSI approach is called monolithic if the fluid and solid problems are solved as one single system where no matching of the data is required at the interface.
In a partitioned approach, there are two different solvers simulating the fluid and the solid part respectively. If the coupling between the solvers is explicit in time then the coupling is loose. The loose coupling has low computational cost, but the simulation may become unstable. To overcome these instability issues, the partitioned problem can be formulated implicit in time, introducing an iteration loop at each time step until a dynamic equilibrium between the fluid and solid is achieved.
Data exchange between the fluid and solid part in this implicit algorithm is called a strong coupling.

Unified Continuum Model
We now specify our ansatz, which corresponds to a monolithic, moving mesh method. An elaborate description can be found in  at full length. Here, we only describe the main features.
Where the size of the vessel is much larger than the size of a red blood cell, the blood flow can be modeled as an incompressible Newtonian fluid (Quarteroni et al., 2014). The governing equations are the Navier-Stokes equations. The dynamic viscosity is chosen as µ = 0.0027Pa · s and the blood density ρ = 1, 060kg/m 3 (Di Martino et al., 2001). In small domains, as the region around the revolute joints of a mechanical heart valve, non-Newtonian effects might have to be incorporated in the model, but these flow features are not targeted in this work.
With the aim of establishing a framework that allows for general formulation and implementation of different models, while applying adaptive error control for realistic 3D applications, a so-called unified continuum model for FSI was developed. The model is described by the conservation laws of mass and momentum for an incompressible continuum, where a stress and phase variable define the properties of the continuum.
Let t ⊂ R 3 be a time-dependent domain with t ∈ I : = [0,t]. Our goal is to determine u(x, t) : t → R 3 , where t encompasses both the solid and the fluid domain and u defines the fluid velocity in the fluid part and the deformation velocity in the structure part: Here τ is the stress tensor and m identifies the mesh velocity in the ALE formulation. In the solid, we choose m to be the material velocity of the structure. In the remaining part of the mesh, m is determined by the mesh smoothing algorithm applied to uphold the quality of the mesh. The constitutive laws are defined via the stress term, where the phase function θ is set to zero in the solid domain and to one in the fluid domain: The kinematic constraint u f = u s is satisfied implicitly by the continuity of the velocity field u for the unified continuum. The dynamic constraint is weakly enforced by applying integration by parts on the stress term and setting it to zero in the weak formulation.
This approach allows us to use the same discretization method, stabilization technique and mesh deformation algorithm as for a pure fluid problem.

Time and Space Discretization
Let 0 : = t 0 < t 1 < · · · < t N : =t be a sequence of discrete time steps, with associated time intervals I n : = (t n−1 , t n ] of length k n : = t n − t n−1 .
We introduce the space-time slab S n : = t n × I n , and let T n = {K} denote the spatial discretization of t n . U n is the discrete velocity, P n is the discrete pressure, and h n specifies the maximal diameter of the cells K ∈ T n .
We choose the finite element function space of piecewise linear functions W n ⊂ H 1 ( t n ), where We identify the discrete solution for velocity and pressure aŝ U = (U, P), the discrete stress for both the fluid and the solid as T , the discrete mesh velocity as M, and the test function aŝ v = (v, q). In time, we choose U to be piecewise linear, and P, v and q to be piecewise constant. Based on these definitions and assuming homogeneous Dirichlet boundary condition for the velocity, the spatially and temporally discretized variational formulation of Equation (3) reads as follows: for each space-time slab S n , find (U n , P n ) : = (U(t n ), P(t n )) with U n ∈ W n 0 and P n ∈ W n , such that: To stabilize the convection dominated problem (3), we use a simplified Galerkin/least-square method, where we drop the time derivative and the diffusion term, and we define SD δ as FIGURE 2 | To detect collision we calculate the distance d ij between the leaflets L j and L i for i, j = 1, 2, 3 (A). As soon as the minimal distance is below a certain threshold, the valve opening is closed. A 2D-surface (blue) is included to model a proper closure (B) and the geometric model of the valve opening is closed by marking the cells directly attached to the 2D-surface as solid (C).

Smoothing Algorithms
Due to the fluid-structure interaction of the aortic valve and the pumping blood flow from the left ventricle, it is crucial for an ALE-method to have a suitable method to adjust an existing mesh. There are different ways to enhance and optimize the quality of the mesh, which may involve e.g., swapping faces and edges, or changing the number of vertices.
Meshing algorithms, which involve change of topology or the number of mesh cells, are not suitable for time-dependent, Frontiers in Physiology | www.frontiersin.org parallel computing. Therefore, it is preferable to use a mesh adaptivity method which omits the necessity or at least minimizes the frequency of remeshing.
To keep a good mesh quality, while limiting the computational cost, our solver combines a linear and a nonlinear mesh smoothing algorithm. The linear smoother accounts for the rough overall re-distribution of the vertices, while the nonlinear smoother optimizes locally the mesh based on the quality of the cells.

Linear Smoother
The linear smoother solves a linear elastic equation in the fluid domain for the mesh velocity, which corresponds to a Poisson equation with Dirichlet boundary conditions given by the structure velocity on the fluid-structure interface, where the vertices are diffusively relocated over the domain. Although it is a simple and fast method, there is no guarantee that improvement is achieved since the equation does not take into consideration the quality of the cells in the mesh.

Nonlinear Smoother
To locally enhance distorted cells, we describe the deformation of the mesh using a nonlinear elasticity problem, and weight the stiffness of the model by a quality measure Q(K) of each cell K in the mesh T n : where d specifies the dimension of the spatial domain and ||.|| F the Frobenius norm. F denotes the deformation gradient between K ∈ T n and a scaled equilateral reference cell.
By weighting the equation by Q(K) and advancing the partial differential equation toward its equilibrium, the mesh is improved toward its goal of optimal shape. A more detailed description is elaborated in .
To limit the computational cost, the nonlinear smoother is stopped after a certain number of "pseudo" time stepsk before a stationary solution is obtained. Depending on the quality of the mesh T n , the total number of pseudo time steps can be adapted to achieve a desired quality.

Modeling of Contact
In order to simulate the closing of a heart valve, an algorithm needs to be implemented to both detect collision and to simulate contact. Our approach is derived from the idea to describe the fluid-structure interaction as a unified continuum. We model contact implicitly by switching fluid cells to solid cells as soon as contact is detected. Collision is detected by solving an Eikonal equation for the distance between two solid surfaces.
In order to detect contact between two leaflets of a native valve, we calculate the minimal distance d min : = min ij,i =j {d ij } between the leaflets L j and L i for i, j = 1, 2, 3, as illustrated in Figure 2A. To model a proper closure of the leaflets, we include a 2D-surface in our volume mesh, which covers the entire valve opening, and as soon d min is below a certain threshold, all cells directly attached under the 2D-surface are marked as solid, as shown in Figures 2B,C. Since the closing moment of a healthy valve is very short, we argue that it is acceptable to cover the whole opening at once. The contact is released at the beginning of the subsequent contraction phase (systole) of the left ventricle. Frontiers in Physiology | www.frontiersin.org

Computational Tools
Nowadays high performance computing is an essential part of computational science. The Heart-FSI solver is implemented in the HPC branch (Jansson, 2013) of the open source FEM library DOLFIN (Logg and Wells, 2010) and the adaptive flow solver Unicorn (Hoffman et al., 2012). Both libraries have successfully been used to efficiently solve large scale industrial problems as described in e.g.,  and Vilela de Abreu et al. (2016).
The simulations were performed on Beskow, a Cray XC40 system, where each node has two CPUs (Intel E5-2698v3) with 16 cores. All volume meshes are created in ANSA (2014), a computer-aided engineering tool for pre-processing.

VALVE MODELS
In the subsequent paragraphs, we describe how we model native and bileaflet mechanical heart valves (BMHV) embedded in the left ventricle and ascending aorta. For each case, we detail the geometry and specify the material as well as the initial and boundary conditions.

Geometry
The geometry of the aortic root has been studied, where geometrical parameters are optimized to resemble the function of a trileaflet valve (Swanson and Clark, 1974). Our model is based on such an optimized geometry proposed by Thubrikar (1990).
We generate a computer-aided design (CAD) model of an idealized native aortic root based on a small set of parameters which can be personalized to a particular patient. The aortic root generator is a set of Python scripts for Rhinoceros 5 (Rhinoceros, 2016) that outputs an aortic root in a fully open valve configuration, as presented in Figure 3A. The model parameters are illustrated in Figures 3B-D. We assume that the aortic root has a threefold symmetry around the z-axis and label the rotational angle by β as depicted in Figure 3B. The plane P A at the annulus and the plane at the sinotubular junction P S are assumed to be parallel. The inner radius at the annulus R A , the inner radius at the sinotubular junction R S and the length of the leaflet in the open position h l , are used to define a truncated cone as shown in Figure 3C. To find the leaflet attachment and the leaflet surface the cone is cut by the plane P c , which is defined by three points P 1 c , P 2 c , and P 3 c , see Figures 3C,D. These points are determined by the height of the commissure h c and the opening angle β. We attach a cylinder with radius R S to the aortic root to model the beginning of the ascending aorta. Geometrical parameters for the sinus of Valsalva are not considered as modifiable yet. The thickness is acquired by copying, scaling and translating surfaces. The parameter values used in the simulations are listed in Table 1.

Material
The leaflets are made of a very thin, flexible and inextensible material. The fibers in an aortic leaflet are aligned in the circumferential direction (Swanson and Clark, 1974), and the mechanical properties vary in different parts of the aortic valve (Kasyanov et al., 1985). In the framework of this work, it is sufficient to assume the solid material to be homogeneous. As material model we choose an incompressible, neo-Hookean material. At this point of development, the material parameters are set to µ s = 3.3 · 10 3 MPa and ρ = 1, 000kg/m 3 . Although these parameters do not conform with realistic values yet, typical characteristics of the flow and valve dynamics can be captured.

Initial and Boundary Conditions
Even though in the initial geometry the valve is in a fully open position, the leaflets are pushed into a starting configuration to facilitate the movement of the leaflets. In order to remove excessive leaflet material resulting from the deformation, we prescribe a constant, initial stress in radial direction such that the material behaves like a contracting balloon which was stretched. The starting position for our simulations with initial radial stress 4 Pa is shown in Figure 4A. The stress is reset for the FSI simulation.
We only consider the two major phases, systole and diastole, and one heart cycle lasts for 1.124 s. The inflow profile is flat and the magnitude is adopted from the left ventricle flow simulations presented in Spühler et al. (2015). At the end of systole, the direction of the inflow is inverted to create a backflow which is physiologically consistent and helps the valve to close. The timedependent inflow magnitude is plotted in Figure 4B. Diastole starts when the valve is closed and the inflow is set to zero. A homogeneous Dirichlet boundary condition for the pressure is set at the outlet.

Geometry
Pathological conditions caused by valvular dysfunction in the form of a narrowing of the valve opening (stenosis) or insufficient closing of the leaflets, reduce the efficiency of the heart. To restore the hemodynamics function, the native heart valve may need to be repaired or even replaced by an artificial implant. Since the first clinical implantation of an artificial valve by Dr. Charles A. Hufnagel in 1952, many different mechanical and bioprosthetic valves have been developed. Due to their wear resistance, the bileaflet mechanical heart valves (BMHV) are most widely favored as aortic valve replacement.
FIGURE 6 | All 2-D images are visualized using these 2-D cuts in Paraview (Ahrens et al., 2005). The plane for the BMHV (A) is defined by its origin in (0.278, −1.65, 1.05) and normal (0, 1, 0), and the plane for the native valve (B) by its origin in (−0.179, 0.0, 1.05) and normal (0, 1, 0). Figure 5A, a typical BMHV is made of a circular housing and two semicircular discs, which are mounted in the housing through a hinge mechanism. Both leaflets are rotating passively in response to the fluid dynamics resulting from the periodic contraction and expansion of the left ventricle.

As can be seen in
Since feasibility, but not clinical validation is the focus of this paper, a detailed geometric model of a mechanical valve is secondary at this stage of investigation. Therefore, the BMHV models are reduced to the leaflets only, embedded in an idealized aorta as depicted in Figure 5B. The geometry is simplified in such a way that no contact between the leaflets occurs.
To use numerical simulations in order to study the flow through a mechanical prosthetic heart valve began in the early 1970s. Since then, many simulations of the flow dynamics around a BMHV have been conducted with the aim to elucidate and eliminate complications as thromboembolism. Simulating flow dynamics in the vicinity of a heart valve is a challenging task. The flow is pulsative and undergoes transition to turbulence. Patient-specific framework and the computational models should account for the multi-scale nature of the flow and deformability of the wall.

Material
We apply the same material model as for the native valve and set the material parameters to µ s = 6.5 · 10 5 MPa and ρ = 1, 000kg/m 3 .

Initial and Boundary Conditions
Contrary to the native valve, we simulate the fluid-structure interaction of the leaflets and the hemodynamics of the left ventricle (LV) conjointly. A detailed description of the boundary conditions for the numerical simulation of the blood flow in the LV can be found in Spühler et al. (2015). We define a rotational axis by fixing two edge points of each leaflet. The hinges on which   Frontiers in Physiology | www.frontiersin.org the leaflets are placed limit the rotational angle so that the BMHV is properly opened and closed. To mimic this mechanism, we set a threshold for the opening and closing angles respectively. As soon as a leaflet exceeds this angular barrier during systole or diastole, its position is locked. The leaflets are released from the fully open position if the mean pressure above the valve exceeds the mean pressure under the valve, and is disengaged from the closed position as soon as a new heart cycle starts. The maximal angular opening is set to 45 • .

RESULTS
In this section we present the numerical results for the native and bileaflet mechanical valves. The 2-D cuts for the native valve and the BMHV are specified in Figure 6. Since we do not model the coronary arteries, which originate from the sinus of Valsalva, and the flow within the aortic root is almost quiescent during diastole, the results are based on the first heart cycle. The results for the native and mechanical valve are presented from meshes with 248 ′ 980 and 783 ′ 823 ′ vertices respectively. At the beginning of the simulation and during diastole, the time step size k n is set such that the Courant-Friedrichs-Lewy (CFL) number is 0.5. During systole, we have to reduce k n such that the mesh smoothing algorithms can maintain the mesh quality. No remeshing is required but in the worst case the CFL number had to be reduced to 0.01 to bypass a sensitive phase of large and fast deformation of the leaflets. To advance the solver one time step, the momentum and continuity equation, the Eikonal equation for contact detection, and the linear and non-linear elasticity equations for mesh smoothing have to be solved. When distributing ∼ 2, 000 vertices per core, each sub-problem is solved in less than 0.5 s but its total time is about 5 s. The latter can slightly vary depending on the quality of the mesh since the cost of the non-linear elastic smoother is higher when the mesh quality is low.

Simulation Results of the Native Valve
First, we examine the opening and closing movement of the aortic valve, which can be divided into four phases (Bellhouse and Talbot, 1969;Labrosse et al., 2010). A rapid valve opening time (RVOT) is followed by a period when the valve stays widely opened (quasi-steady phase). The valve first closes steadily and then rapidly due to reversed flow (RVCT) in the very end of systole. All these stages can be observed in our simulations by measuring the geometric orifice area (GOA), which is calculated by determining the area of the surface used for closing the valve. The time-dependent GOA is depicted in Figure 7 and matches well the dynamics captured in Labrosse et al. (2010) During RVOT, the fluid is accelerated over the whole domain flowing toward the outlet. As observed in De Hart et al. (2003), even the blood residing in the sinus cavity is washed out as shown in Figures 8A-C. During the subsequent period, as the valve reaches and stays in the fully opened position, the flow is dominated by a strong, central jet. The flow starts to decelerate at about t = 0.2 s when the valve is still completely opened, and at about t = 0.25 s the flow in the sinus cavity does not flow toward the outflow anymore, see Figure 9A. A small vortex starts to form at the tip of the backside of the leaflet, as depicted in Figure 10. Computing Lagrangian coherent structures, (Shadden et al., 2010) can distinguish two flow domains in this phase of deceleration. They observe a boundary between the strong  outflowing jet and the regions with recirculating flow. These features can also be observed in our simulations as visualized in Figure 11.
During the closing phase, two different vortices can be observed, as shown in Figures 9B,C, 12. One vortex is located just above the leaflet and the other one within the sinus cavity. Although they are rotating in counter directions, both of them drive the valve to close. The vortex within the sinus cavity merges to a streamline flow and only the vortex at the tip of the leaflet is left. A fully reversed flow in the ascending aorta is modeled by altering the inflow condition and a complete closure of the valve is achieved.
High stress has been connected to leaflet damage and failure. To analyze the stress distribution in the leaflets of our model, the von Mises stress τ v in logarithmic scale is computed for the same time instances as the velocity and pressure fields in Figures 8, 9, We also visualize the stress distribution at the moment when the valve has just been closed at t = 0.442. As can be observed in Figure 13, regions with high stress can be localized to the attachment lines, commissure and leaflet belly. However, due to the low mesh resolution, this is only a qualitative analysis of the stress distribution.
No elaborated studies to analyze mesh sensitivity have been conducted yet. So far, we have only investigated to what extent the point of contact is affected by mesh refinement. For this purpose, the mesh is uniformly refined in the vicinity of the aortic root and we observe that the point of contact does slightly differ as listed in Table 2.

Simulation Results of the BMHV
To examine the valvular kinematics, we calculate the opening and closing angle as well as the rotational velocity of the leaflets. The rotational angle of the right and left leaflet is defined as depicted in Figure 14A, and the results are presented in Figure 14B. We observe that both leaflets are slightly open at first and accelerate and decelerate linearly while opening. The right leaflet precedes the left leaflet in the opening phase. This kinematic variation is of course strongly influenced by the geometry of the aorta. They then stay in their fully opened position until they close very rapidly, mainly due to backflow.
The geometry of a BMHV generates three jets, namely one central jet flowing through the gap between the leaflets and two side jets. During the end of the rapid opening phase, vortex rings are shed from the tip of the leaflets due to the difference in the velocity magnitude of the central jet and the two side jets. The vortex rings travel downstream a short distance before they vanish. Snapshots of the velocity field using line integral convolution (LIC) are visualized in Figure 15A. Figure 15B provides a closer view of the recirculation areas We use the open source code Saaz to calculate λ 2 for our simulations (King et al., 2011). The threshold λ 2 is manually adjusted until we can differentiate coherent vortex structures as shown in Figure 15C. The velocity vectors are added to indicate the rotational direction. The vortex observed at t = 0.1 at the right leaflet merges after a very short time into a recirculating flow with opposite direction (t = 0.11) and separates from the leaflet (t = 0.115). Meanwhile, a clockwise vortex is developed at the outer part of the left leaflet (t = 0.12), which eventually entails a neighboring, counter-clockwise rotating flow (t = 0.124). The former is swept off downstream, while the latter stays attached to the leaflet. When the valve has reached the fully opened position, no further vortices are developed.

CONCLUSION
The aim of our research is to develop an open source modular framework for modeling and simulating the blood flow in the heart. In the present work we place prototypes of a native and mechanical aortic valve between the left ventricle and the aorta.
We model both the fluid-structure interaction of the valve and the contact problem in the framework of a unified continuum. This approach to simulate the valvular dynamics is unique and has the advantageous properties that the whole problem can be described by a set of partial differential equations for which the same numerical methods are applicable. Furthermore, no instability issues due to the fluid-structure coupling is encountered. All algorithms are implemented in the FEniCS-HPC software framework optimized for parallel computing.
We generated a CAD model of an idealized native aortic root based on a small set of parameters, where we leave for future work to adapt the geometry to patient-specific data, and to connect the native aortic root to the left ventricle. The bileaflet mechanical heart valve is reduced to the leaflets only, which are embedded in a simplified geometric model of the aorta. In contrast to the native valve, we simulate the fluid-structure interaction of the leaflets and the hemodynamics of the left ventricle conjointly. The next step is a more realistic geometric model of the BMHV.
The weak point of our approach is the degradation of the mesh quality under large mesh deformations. All the simulations were conducted without remeshing, but we usually had to reduce the time step such that the mesh smoothing algorithms could comply with the deformation. The small time step size increased the computational time. Remeshing the volume mesh is an alternative, but not ideal for parallel computing. Thus, this limitation has to be addressed.