An interface‐tracking unified continuum model for fluid‐structure interaction with topology change and full‐friction contact with application to aortic valves

An interface tracking finite element methodology is presented for 3D turbulent flow fluid‐structure interaction, including full‐friction contact and topology changes, with specific focus on heart valve simulations. The methodology is based on a unified continuum fluid‐structure interaction model, which is a monolithic approach, where the fundamental conservation laws are formulated for the combined fluid‐structure continuum. Contact is modeled by local phase changes in the unified continuum, and computational results show the promise of the approach. The core algorithms are all based on the solution of partial differential equations with standard finite element methods, and hence any general purpose finite element library which can leverage state of the art hardware platforms can be used for the implementation of the methodology.


Heart valve simulation
In the wake of medical imaging, high-performance computing, and big data, in silico medicine has emerged as a powerful tool for personalized medicine. With the introduction of fractional flow reserve (FFR) by computed tomography, the US Food and Drug Administration in 2014 granted the first approval of cardiovascular simulations for routine clinical use, to assess the functional severity of coronary artery stenosis. Clinical trials have demonstrated good agreement with invasive measurements, and a significant reduction in the number of invasive procedures performed on patients. 1 Heart valve disease constitutes a major challenge for society, for which computational methods have shown promise to improve clinical diagnosis and treatment. Patient-specific simulations have the potential to be used both in risk assessment for valve repair or replacement, for personalized clinical interventions, and for the design of prosthetic valves. 2 Quantities of interest in a heart valve simulation may be the deformation of the opening and closing valve tissue and the associated internal and surface stresses, and the blood flow dynamics including shear stress, pressure, localized turbulence, and flow stagnation. Needless to say, a detailed analysis of such quantities requires a fully coupled fluid-structure interaction (FSI) model of the valve with the ability to simulate both laminar and turbulent flow.

Purpose and scope
To model a patient-specific heart valve, there are multiple challenges related to data acquisition, uncertainty quantification, geometry description, and boundary conditions. But here we focus exclusively on the FSI aspects of the simulation problem, which include the discretization of deforming structure models with contact and topology changes, coupling the structure models to fluid models such that fundamental conservation laws are respected, and the simulation of flow separation and turbulence. For simulations to be practical, the numerical methods must be robust and accurate, and the computational algorithms must be efficient. Specifically, the algorithms need to be designed for modern hardware platforms and software environments. The purpose of this article is to present an FSI methodology which addresses the points above.

Arbitrary Lagrangian-Eulerian methods
Traditionally, the equations of structure (solid) mechanics are based on a Lagrangian (material) description of the material, where each coordinate represents a material point which is tracked in time. By contrast, the equations of fluid mechanics typically use an Eulerian (spatial) description, where each coordinate corresponds to a spatial point, where the flow is recorded over time. In an Arbitrary Lagrangian-Eulerian (ALE) method, the conservation laws of an FSI model are expressed with respect to an arbitrary reference coordinate system, connected to the discretization of the equations. Often a Lagrangian description is retained for the structure mechanics, so that the ALE reference coordinate system in the structure mechanics model is equal to the material coordinate system. ALE methods for fluid dynamics were first developed in the context of finite difference methods, [3][4][5][6][7] while finite element methods were used to model structure mechanics. Inspired by this work, ALE methods were then developed also for finite element methods, [8][9][10] opening for monolithic methods based on finite element discretization of the complete FSI system in one single computational mesh. The main challenge of ALE methods is the deformation of the fluid dynamics mesh as a result of structure displacements, ranging from small deformations to rotation, contact, and topology changes. Mesh smoothing techniques are designed to handle moderate mesh distortion without changing the mesh connectivity, while large deformations of the mesh require methods based on local or global remeshing. 11,12 For a detailed description of ALE methods see one of the several reviews available. 13

Immersed boundary methods
A different approach to the FSI problem is taken for immersed boundary (IB) methods, a class of methods that originate from the work by Peskin, 14 inspired by Viecelli. 15 Here the structure is immersed in the fluid, and as the structure is convected by the fluid it exerts forces on the fluid as a response to its own deformation. The fluid mesh is not deformed since an Eulerian coordinate system is used to formulate the fluid mechanics equations. Hence, mesh smoothing or remeshing is not necessary, even under large deformations. The IB method was developed in the context of finite difference methods which are not well suited for mesh deformation, and thus the IB method offered a way to circumvent this problem. On the other hand, the indirect representation of the structure leads to a loss of accuracy compared with ALE methods. For more information, see one of the multiple reviews on the IB method. 16 Lagrangian coordinates for the structure and Eulerian coordinates for the fluid. FD methods were first developed for thin-walled structures, 18 but were later extended to elastic solid bodies. 19,20 The fluid-structure coupling constraints are enforced by a Lagrange multiplier, which can be interpreted as the surface force exerted on the fluid and structure over the fluid-structure interface. In the discretization of the resulting saddle-point problem, the finite element approximation space of the Lagrange multiplier must satisfy an inf-sup condition. 21 The main strengths and weaknesses of the FD method are similar as for the IB method. [22][23][24]

Overlapping domain methods
The ALE method is an example of an interface tracking technique, where the fluid-structure interface is explicitly represented by a matching computational mesh. By contrast, the IB and FD methods are interface capturing techniques, in which the effect of an implicit fluid-structure interface is distributed over a nonmatching computational mesh.
Overlapping domain (OD) methods aim to combine the strengths of the two techniques, by combining a background mesh for the fluid with an embedded mesh that contains the structure and extends into the fluid with an explicit representation of the fluid-structure interface. Hence, OD methods provide an improvement in accuracy compared with IB and FD methods with respect to the approximation of the geometry of the structure and its impact on the fluid dynamics.
The challenge for OD methods is instead the fluid-fluid coupling of the background mesh and the embedded mesh, which can be accomplished, for example, by a domain decomposition technique, Nitsche's method for nonmatching interfaces, or through blending of the approximation spaces by a partition of unity method. 25-29

Contact models
Since our target is heart valve simulation, the ability to model contact is of key importance. The contact problem can be divided into collision detection and contact modeling, where the former refers to an algorithm that activates the latter. Collision detection is often based on either testing for intersecting structures, or computing the shortest distance between structures. 30 Once collision is detected and the contact surface is determined, the contact model prescribes how the colliding structures should interact with each other, typically in the form of inequality constraints on the contact surface stresses, enforced, for example, by penalty methods or Lagrangian multipliers. A frictionless contact model only takes normal stresses into consideration, whereas a frictional contact model also includes tangential stresses. Clearly, the methods used to discretize the colliding structures strongly impact the performance and accuracy of both the collision detection algorithms and the contact models. 31 Specifically, the Lagrange multiplier saddle-point problem needs to be addressed. 32 In the context of FSI, specific challenges arise. For interface tracking ALE methods, the extreme mesh deformation associated with contact must be addressed, for example, by enforcing a minimal distance between colliding structures. 33 On the other hand, interface capturing methods suffer from the implicit representation of the surface geometry of the colliding structures. 34,35 More recently the idea of a contact domain, or a contact medium, was suggested. [36][37][38] Here the space between two colliding structures is filled by a third medium, from which collision can be detected and contact models enforced.

Heart valve simulation
The four valves in the heart open and close in response to pressure changes in systole and diastole, the contraction and relaxation of the heart muscle, to allow blood to flow between the chambers of the heart without backflow. The aortic valve and the mitral valve, both connected to the left ventricle, most often require clinical intervention due to incomplete valve opening (stenosis) or closing (regurgitation). Blood flow through the aortic valve is characterized by complex flow and turbulence, [39][40][41] with thin leaflets that undergo large deformations and contact, hence a challenging FSI problem. State of the art FSI methods for 3D heart valve simulation include methods based on ALE in the form of integration of a space-time variational multi scale, slip interface, topology change and isogeometric discretization 42 , and variants of the IB and FD methods. 43,44 Some direct comparisons of the different methodologies have been performed for model problems, which indicate that interface tracking methods in general show a higher accuracy than interface capturing methods, 45,46 and that monolithic methods are more accurate than weakly coupled partitioned methods. 47 Flexibility is the main counterargument for interface capturing methods, to circumvent the mesh deformation problem, and flexibility is also the reason a partitioned approach is popular, since separate methods and software can be chosen for the fluid and structure models.

Unified continuum FSI
A monolithic interface tracking method for heart valve simulation is presented, in which the idea of a contact medium is extended to FSI to allow for contact and topology changes. The method builds on the unified continuum FSI (UC-FSI) methodology, 48,49 where the fluid and structure subdomains are discretized by a finite element method as one single unified continuum. Preliminary work has shown promise in the area of biomechanics, 50,51 and in this article a detailed description is given of the methodology with a specific focus on heart valve simulation.
The main idea to model contact is to let the fluid subdomain function as a contact medium for colliding structure subdomains. In this contact medium, Eikonal equations are solved to generate global distance fields from which collision can be detected. Once two structures are in contact based on a threshold condition with respect to the distance fields, the contact model is activated so that the continuum phase of the contact medium is changed from a fluid phase to a contact phase in a certain contact region. The phase change is realized by locally modifying the constitutive law for the mesh elements in this contact region. In fact, if the contact phase is defined by a structure constitutive law, this corresponds to a topology change, where the two colliding structures are treated as one single structure, composed of the two colliding structures and the activated contact medium, possibly with different constitutive properties.
The stabilized finite element discretization of the fluid dynamics model can be interpreted as an implicit large eddy simulation, or with sufficient mesh resolution a direct numerical simulation. Hence, both laminar and turbulent flow can be simulated by this model. 52,53

Implementation
All the core algorithms of the methodology build on the solution of standard partial differential equations with the finite element method. Therefore, the methodology can be implemented in any general purpose finite element software. The computational results we report in this article are generated from an implementation in the open source software FEniCS-HPC, which shows near optimal parallel scaling up to tens of thousands of compute cores. 54,55 Specifically, the partial differential equations to solve are the conservation laws for the unified fluid-structure continuum, elasticity equations for mesh smoothing, and Eikonal equations to determine distance fields between the structures of the model.

Outline
In Section 1 the motivation and background for the article are given, and in Section 2 the unified continuum fluid-structure interaction model is described, with specific focus on the space-time finite element formulation, the mesh smoothing algorithm, and the contact model. Computational results are presented in Section 3, both for a model problem and an aortic valve. The article concludes with Section 4, where the results are discussed and put into context.

Unified continuum FSI model
Let Ω analogous to a multiphase flow problem. The unified fluid-structure continuum is here assumed to be F I G U R E 1 The unified continuum fluid-structure interaction model, with the map −1 t ∶ Ω 0 → Ω t from the Lagrangian coordinate system to an ALE reference coordinate system corresponding to the accumulated mesh deformations at time t isothermal, so that no energy equation is necessary, and incompressible with constant densities over each phase, so that mass conservation is reduced to the constraint of a divergence free velocity field. These approximations are common in the field of heart valve simulations.
At initial time t = 0, the geometry of the unified continuum Ω 0 is assumed to be a polytope domain which is discretized by a conforming simplicial finite element mesh  h,0 = {K}. If Ω 0 cannot be represented by a simplicial mesh, a geometry approximation error is introduced at time t = 0. Each element K ∈  h,0 is either a fluid or a solid element, determined by a phase function defined over the unified mesh being the collections of fluid and solid elements, respectively. As the unified continuum domain Ω t deforms over time, the vertex coordinates of the corresponding mesh  h,t change, but the mesh connectivity stays constant. Hence, there exists a bijective map between the evolving fluid and structure subdomains Ω for any K i ∈  h,t f and K j ∈  h,t s . Over time the mesh  h,t is deformed according to a mesh velocity m(x, t), which for the structure submesh  h,t s is determined by an approximation of the material velocity, and for the fluid submesh  h,t f is given by a mesh smoothing algorithm. Using standard notation, 13 the mesh velocity defines a map −1 t from the Lagrangian coordinate system corresponding to Ω 0 to an ALE reference coordinate system at time t corresponding to the accumulated mesh deformations (see Figure 1), The UC-FSI model can then be expressed in ALE form as where (x, t), u(x, t), p(x, t), (x, t), f(x, t) are the density, velocity, pressure, Cauchy stress tensor, and external force. Constitutive laws for the fluid and structure phases are defined through the density and the Cauchy stress tensor , which is decomposed into a deviatoric stress tensor and a mechanical pressure p = −1∕3tr( ), with tr( ) being the trace of the Cauchy stress tensor, and I the identity tensor. The mechanical pressure acts as a Lagrange multiplier for the constraint (4) in both the fluid and structure subdomains, but the deviatoric stress tensor is different in the two subdomains.
If a phase function is defined by the deviatoric stress tensor can be expressed as with f and s being the fluid and structure deviatoric stress tensors, respectively, and similarly for the density Note that the phase function (x, t) is a discontinuous piecewise constant function over the elements of the mesh  h,t .

UC-FSI constitutive models
The UC-FSI model described above combined with suitable constitutive models can be used to simulate a wide range of fluid-structure interaction problems. In this article, constitutive models are chosen with the purpose to simulate fluid-structure interaction in an aortic valve. Blood flow is known to exhibit non-Newtonian properties, but at the scale of ventricular blood flow these effects are typically small. Therefore, in the context of heart valve simulation often a Newtonian model is used, where f is the dynamic viscosity of the blood, and is the strain rate tensor. Hence, with this model the constitutive parameters for the blood are f and f . The aortic valve consists of complex anisotropic material, but here we make the simplifying assumption that the aortic root and the leaflets can be approximated by a neo-Hookean solid, for which we use the following rate formulation, where s is the shear modulus. Therefore, the constitutive model of the aortic valve is determined by the density s and the shear modulus s . In summary, the UC-FSI conservation laws (3)-(4) take the form together with appropriate boundary conditions and initial conditions, where the fluid and structure constitutive laws are selected based on the phase function (x, t) defined by (6). The coupling conditions of the fluid-structure interaction problem take the form of constraints over the fluid-structure interface Γ fs t , enforcing continuity in the velocity field (the kinematic constraint) and the normal stresses (the dynamic constraint) .

Space-time formulation of UC-FSI
In a semidiscretization approach, where the UC-FSI model is discretized with respect to its spatial variable using a finite element method, the result is a system of ordinary differential equations to be solved by a suitable time stepping domain Ω t and its associated mesh for each t ∈ I n = (t n−1 , t n ) (left); and the corresponding space-time mesh over the slab where (x,t) = (x n−1 , t) are the coordinates of the tensor product space-time slabŜ n = Ω n−1 × I n (right) method. Note however that since the domain Ω t deforms over the time subinterval I n , the finite element integrals must be recomputed for each t ∈ I n , even the integrals that correspond to the linear terms in (12)-(13) (see Figure 2). Alternatively, a space-time finite element method 56-58 is chosen to discretize the space-time domain with Q f and Q s being the naturally defined fluid and structure space-time subdomains.
For each subinterval I n = (t n−1 , t n ), denote by x n the spatial coordinates of the domain Ω n = Ω t n , and by x n−1 the spatial coordinates of the domain Ω n−1 = Ω t n−1 (see Figure 2). Then define a coordinate map F n ∶Ŝ n → S n , from the tensor product reference space-time slabŜ n = Ω n−1 × I n with coordinates (x,t) = (x n−1 , t), to a tilted space-time slab S n with coordinates where the corresponding Jacobian matrix of the spatial map, denoted by J n , is defined by its components with ij being the Kronecker delta which represents an identity matrix. Hence, the spatial integration measure over S n can be expressed in the local reference coordinates as with | det(J n )| being the absolute value of the determinant of J n , and any function v(x, t) can be written in reference coordinates asv with the gradient∇v = J −1 n ∇v. Remark 1. The map F n can be used to determine the ALE reference coordinate system for any t ∈ I n , with the map −1 t constructed from the recursion (x n , t n ) = F n (x n−1 , t n ).

Finite element approximation
For each t ∈ I, define a finite element space of piecewise polynomial functions of order k over Ω t , and define the corresponding space-time finite element space over each tensor product space-time slabŜ n , The space-time finite element space over the tilted space-time slab S n is expressed by the map F −1 n ∶ S n →Ŝ n , Global finite element spaces can then be constructed over the space-time domain Q, either discontinuous or continuous between the space-time slabs S n , where the number of degrees of freedom on each S n is lower for a continuous approximation, since By contrast, a Neumann boundary condition for the stress, is implemented through the variational form.
With (⋅, ⋅) X the standard L 2 (X) inner product, a finite element method to approximate the UC-FSI system (12)-(13) can then be expressed as: find U ∈ [V (1,1) h,g D (Q)] d and P ∈ W (1,0) h (Q), such that for all v ∈ [W (1,0) h,0 (Q)] d and q ∈ W (1,0) h (Q). The velocity approximation is continuous piecewise linear in space and time, and the pressure approximation is continuous linear in space but discontinuous piecewise constant in time. A time stepping approach is used to compute U n and P n for each subinterval I n , where U n (x) = U(x, t n ), and P n (x) = P(x, t) for all t ∈ I n . For each subinterval I n , the approximation of the discontinuous piecewise constant stress tensor T ∈ [W (0,0) h (Q)] d×d is computed from (11) by a trapezoidal time stepping method, whereas the mesh velocity M ∈ [V (1,0) h (Q)] d is set equal to the material velocity M n = U n in Q s , and determined by a mesh smoothing algorithm in Q f . For the nonlinear terms in (27) a midpoint quadrature rule is used for the time integrals, which introduces the mean velocity U n = 0.5(U n + U n−1 ) in the convection velocity over each I n . The resulting nonlinear algebraic system is solved by a quasi-Newton approach, using Krylov subspace methods to solve the linear systems of equations. In our monolithic approach it is crucial to keep a good mesh quality. Therefore, a small time step is chosen. Applying a quasi-Newton approach contributes to reduce the simulation time since the Jacobian does not need to be assembled at each time step. The Krylov method used is the generalized minimal residual method (GMRES) with an algebraic multigrid (AMG) preconditioner.
The kinematic constraint of continuity in the velocity field is satisfied by construction in the UC-FSI model, whereas the dynamic constraint of continuity in normal stresses is satisfied weakly through the variational form (27), due to the jump terms that appear for all internal element facets in the mesh as a consequence of partial integration in the stress term. 59 Remark 2. For the computations in this article, the Jacobian J n in (27) is approximated by the identity matrix, which simplifies the formulation and makes the computations of the integrals cheaper, but introduces a quadrature error.

Stabilization
The role of the stabilization term SD (U, M, P; v, q) is both to stabilize the system with respect to its saddle-point nature to allow for equal order interpolation of the velocity and pressure approximations, and to avoid spurious oscillations if the system is convection dominated. Here a reduced form of a Galerkin least squares (GLS) stabilization is used, where = ( 1 , 2 ) are stabilization parameters, Note that in the structure subdomain Q s the convection velocity U − M is small, so that which is the part of the stabilization that provides the necessary stability with respect to the saddle-point nature of the problem. In the fluid subdomain Q f , standard stability estimates for the method can be proven. 60 Methods for heart valve simulation must be capable to model both laminar and turbulent flow, since blood flow in the heart is characterized by high-velocity jets through the valves and vortex breakdown in the chambers. Further, to capture the deforming chambers and the pulsating flow, computational resolution of the dynamics is necessary, which disqualifies traditional turbulence modeling based on computing statistical averages.
The stabilization of the fluid flow equations in the UC-FSI model (27) can be interpreted as a residual based turbulence model, 52 which does not rely on any flow dependent parameters and, therefore, does not need to be reparameterized for different flow regimes. This methodology is well tested in various fields of science, 53 including for the simulation of ventricular blood flow and heart valves. 61-63

Mesh smoothing
For each space-time slab S n the mesh velocity M n is constructed from the material velocity U n in the structure subdomain S s n , and by mesh smoothing in the fluid subdomain S f n . The quality of each element K in the mesh can be quantified by defining a transformation map and its associated Jacobian J K with respect to a reference simplexK, here chosen as the simplex with one vertex at origo and the other vertices at hKê i for each unit Cartesian basis vectorê i , with hK being the diameter of K at time t = 0 (see Figure 3).
The Frobenius norm || ⋅ || F and the determinant det(⋅) are used to define the following scale invariant quality measure, 64 F I G U R E 3 Map from a reference triangleK to an element K and its associated Jacobian J K (left), and a triangle mesh with each element K colored based on its quality measure Q(K), with darker colors for higher values (right) with Q(K) ≥ 1 for noncollapsed elements, and optimal quality Q(K) = 1 for orthogonal Jacobian matrices J K , corresponding to transformation maps in the form of translations and rigid rotations. An associated quality function q(x) is defined by Here two mesh smoothing methods are described, one relatively cheap method which corresponds to a diffusive redistribution of the mesh velocity by solving a Poisson equation, and one more expensive method based on an elastic analogy, where the fluid mesh is modeled as an elastic solid. The diffusive mesh smoothing method takes the form: find , with the Dirichlet boundary condition In the elastic analogy, the mesh is modeled as an elastic solid under relaxation from internal stresses due to strains in each element K with respect to the reference simplexK. Hence, the level of local strain in an element K can be directly connected to the quality measure Q(K), with zero strain if Q(K) = 1, corresponding to translations and rigid rotations with respect toK.
The elasticity equations are formulated by introducing a pseudotime variable s > 0, with an associated pseudotime derivative D s = d∕ds. A semidiscretization approach is employed, where a finite element method is used to discretize the spatial variable, which leads to the following system of ordinary differential equations: for each s > 0, find M n ∈ [V (1) with the K subscript indicating component K of the piecewise constant functions. The initial condition for the mesh velocity is taken from the previous space-time slab, M n (s = 0) = M n−1 , and the deformation tensor is initialized by the Jacobian J K acting as a deformation gradient, that is, Note the inverse of the quality function q −1 (x) which acts as a weight for the strain I − B −1 in (36), with the effect that for the elements of poorest quality the mesh velocity is kept small, to equidistribute the quality of the elements over the mesh.
The system (36)-(37) is solved for a number of pseudotime steps by a trapezoidal method until the maximal Q(K) for all K is below a given tolerance, after which the resulting mesh velocity M n is used to deform the mesh Ω In order to keep a good mesh quality, while limiting the computational cost, we combine the "Poisson" and the "elastic analogy" mesh smoothing algorithms. First, the simpler smoother accounts for the rough overall redistribution of the vertices, while the elastic smoother locally optimizes the mesh based on the quality of the cells.
Remark 4. The reference simplexK is here chosen to be equal to a scaled reference finite element, which enables the reuse of standard finite element mapping functions to construct the Jacobian J K . The choice of the scaling factor hK means that each element K seeks to keep its size over the time interval I. Remark 5. Equation (37) is derived from the relations B K = J K J T K and D s J −1 K = −J −1 K ∇M n K , with B K analogous to a left Cauchy-Green deformation tensor based on the deformation gradient J K , and with ∇M n K being a material velocity gradient.

Contact model
To model contact between structures, the first step is to know when to activate the contact model, that is, collision detection. In a UC-FSI model, the fact that the structures are embedded in a monolithic unified continuum mesh is used to compute global distance fields, by which collision can be detected. Specifically, for each t ∈ I the distance d i (x, t) from a structure surface Γ i in the fluid subdomain Ω f t satisfies an Eikonal equation which can be expressed as with the boundary condition d i | Γ i = 0. Consequently, for each t ∈ I n a discrete distance field D n i (x) can be computed by the finite element method: is an artificial viscosity which is used to stabilize the finite element method, and the resulting nonlinear algebraic system is solved by a quasi-Newton method.
The distance between two structure surfaces Γ i and Γ j can then be computed from their respective distance fields as so that when d ij (x) < d TOL for some x ∈ Ω f n−1 , the contact model is activated over a contact region Ω c n−1 . The choice of the threshold d TOL > 0 and the contact region Ω c n−1 ⊂ Ω f n−1 determines how local the contact model is, which influences the accuracy, efficiency, and stability of the method. For the examples in this article, the choices of thresholds and contact regions are specified together with the rest of the problem descriptions. Once collision is detected, the phase function n (x) = (x, t n−1 ) is modified for all x ∈ Ω c n−1 , to reflect a phase change in the contact region (see Figure 4). The contact phase could be the same as in the structure part of the domain Ω s n−1 , but also a different new contact medium to model various forms of friction contact. In this article we set the contact phase equal to the phase of the structure part of the domain, corresponding to a full-friction contact, where no sliding contact is possible. The development of the UC-FSI contact model is inspired by previous work in the field, 36-38 but the context of 3D turbulent flow fluid-structure interaction is novel with several promising results in the biomechanics area. 50,51,63

Implementation
The core algorithms of the UC-FSI methodology all rely on the solution of partial differential equations by standard finite element methods, which means that in principle it can be implemented in any generic finite element software. Although, for the problem class targeted in this article, 3D turbulent flow fluid-structure interaction, the bare size of the computational problem demands that the hardware platform is sufficiently powerful and that the software implementation is capable enough to leverage this computational power.
In this article, the reported simulation results are generated by a software implementation in FEniCS-HPC, 54,55 based on a hybrid MPI/OpenMP parallel programming model and the PETSc linear algebra library. 65 The hardware platform was Beskow, a Cray XC40 system, where each node has two cpus (Intel E5-2698v3) with 16 cores, and the simulations are performed based on unicorn-0.

Mesh sensitivity study of the UC-FSI contact model
In previous work the UC-FSI model without contact has been validated for benchmark problems in both 2D 48 and 3D, 49 which is here complemented by a study of UC-FSI with contact for a 3D model problem under mesh refinement. The 3D model problem takes the form of a bouncing sphere, inspired by a related 2D study. 66  No-slip (zero velocity) boundary conditions are used at the bottom and side walls, and at the top wall a zero pressure boundary condition is used, modeling a tank without a lid. In Figure 5 the problem is described in more detail, with the constitutive parameters of the sphere and the fluid in Table 1. A constant time step size k n = 0.7 × 10 −3 is chosen, which is a conservative but robust choice in order to ensure the quality of the mesh during the simulation. Note that this is just a model problem constructed to investigate the properties of the UC-FSI contact model, for which the constitutive parameters do not translate to any specific physical experiment. Therefore, we do not use explicit units in the text or in Table 1.
The contact model is determined by the collision detection threshold set to d TOL = 0.06, and the contact region Ω c n−1 ⊂ Ω f n−1 is chosen as all elements for which the collision detection threshold is satisfied and the velocity is negative in the vertical direction. Based on preliminary test simulations, this choice of d TOL was the optimal value for the chosen material property of the contact phase, in order to ensure bouncing with good mesh quality. In the activated contact region the continuum phase is changed from the fluid phase to the contact phase, which is here identical to the structure phase. When a positive vertical velocity is detected in an element K ∈ Ω c n−1 it is removed from the contact region, to avoid a sticking effect, and its continuum phase is changed back to a fluid phase so that K ∈ Ω f n−1 . In Figure 6 the effect of the contact model is illustrated for the bouncing sphere, where the contact region Ω c n−1 appears and disappears based on the collision detection threshold d TOL . The stresses in the sphere are visualized by the von Mises stress from which we observe a high local stress at impact, which generates compression waves that travel back and forth in the vertical direction. A qualitatively similar phenomenon was reported for the von Mises stress in the case of a 2D elastic disk. 66 After impact and the subsequent bounce off the ground, the sphere recovers its original shape due to its elastic properties (see Figure 7). Hence, the UC-FSI contact model behaves qualitatively as expected, both reproducing the elastic recovery of the shape of the sphere after the bounce, and with respect to the propagation of internal stresses in the sphere.
To investigate the sensitivity of the UC-FSI contact model with respect to mesh resolution, a mesh sensitivity study is performed, where the vertical coordinate of the center of mass of the sphere is tracked over four bounces. The finite element mesh is locally refined four times in a region below the sphere. Apart from the discretization error, the mesh resolution also impacts the contact region which in turn influences the contact model. In Figure 8 the center of mass of the sphere is plotted for the five different meshes, where the results for the three finest meshes show only a small sensitivity to the mesh resolution. A slight decrease in the maximal height can be observed over time, which is expected since the sphere experiences a drag force from the fluid in which it is immersed. This decrease in height is much more pronounced for the coarsest mesh, which likely is due to a stronger drag force resulting from artificial numerical diffusion, which disappears with finer resolution.

Simulation of an aortic valve
We here focus on the UC-FSI methodology and its suitability for heart valve simulations, whereas a detailed analysis of the simulation results in a biomedical context is given in previous work. 63,67 The geometry of the aortic valve was F I G U R E 7 The deformation of the sphere in 2D cut-planes, plotted for the same time instants as presented in Figure 6, showing the deformation both before (left) and after (right) the impact Based on a small set of parameters, 63 an idealized native aortic root is created for which the height is 250 mm, the inner base radius at the annulus is 20 mm, the inner radius at the root 22 mm, and the thickness of the leaflet is 1 mm. We use the incompressible Navier-Stokes equations to model the blood flow, where we set the dynamic viscosity to f = 0.0027 Pa⋅s and the blood density to f = 1060 kg/m 3 . 69 A critical part of our approach is to maintain the quality of the deforming mesh under deformation without remeshing which requires small time steps, here a constant time step of k n = 10 −3 seconds is chosen. A homogeneous Dirichlet boundary condition is used for the material velocity of the aortic root, and a homogeneous Dirichlet condition for the blood pressure is used at the outflow boundary. The two main phases of the heart cycle are systole, in which the heart muscle contracts and blood flows from the left ventricle through the open aortic valve, and diastole, where the heart muscle relaxes and the aortic valve closes. In systole the inflow blood velocity boundary condition is constructed from the outflow in a left ventricle simulation, 63 which at the end of systole is inverted to create a backflow that assists in closing the valve, consistent with heart physiology. At initial time t = 0 seconds, the starting positions of the three leaflets L i (i = 1, 2, 3) are shown in Figure 9, which also correspond to the stress free configurations of the leaflets. The surface that connects the leaflet edges is marked in Figure 9, here referred to as the valve opening.
The minimal distance between the leaflets is computed as where d ij (x) = D n i (x) + D n j (x), and D n i (x) is the discrete distance field constructed by solving an Eikonal equation with a zero Dirichlet boundary condition on the ventricle side of leaflet L i (see Figure 10). Here the contact region Ω c n−1 is determined a priori as all elements that are attached to the valve opening from the ventricle side (see Figure 9). When the minimal distance is shorter than a threshold d min < d TOL , the contact model is activated so that the phase of all elements in the contact region is changed from the fluid phase to a contact phase, here identical to the structure phase with the same properties as the leaflet material. Similar to the previous application case of the bouncing sphere, the threshold value is heuristic and here set to d TOL = 3 mm. At the onset of systole the contact model is deactivated. Mesh sensitivity of the contact model activation time is studied for three different meshes, uniformly refined in the vicinity of the aortic root. It is observed that the difference is less than 2% (see Table 2). In Figure 11, the finite element mesh is depicted in open and closed positions for "Mesh 3" (248 ′ 980 vertices).
In order to understand the dynamic behavior of the leaflets, we measured the geometric orifice area (GOA) over time. The time-dependent GOA is depicted in Figure 13, reproduced from our previous work. 63 In Figure 12 the opening and closing of the aortic valve is visualized, which compare well to the expected phases described in the literature; 70,71 a rapid opening of the valve, ca. 10% of systole, which then stays wide open for ca. 50% of systole, after which the valve closes steadily under ca. 30% of systole and finally achieves full closure under ca. 10% of systole with the help of backflow and a vortex formed at the tip of the leaflet (see Figure 13).

DISCUSSION
In this article, an interface tracking finite element methodology is presented for fluid-structure interaction problems characterized by 3D turbulent flow, contact, and topology changes. Specifically, heart valve simulations fall into this category. The methodology is based on a unified continuum fluid-structure interaction model, where contact is modeled by local phase changes. Previous computational studies have validated the unified continuum methodology for fluid-structure interaction, including turbulent flow, and from the computational results presented in this article it is clear that the methodology can be extended to also model contact and topology changes. The core algorithms are all based on the solution of partial differential equations with standard finite element methods. Hence, any general purpose finite element library which can leverage state of the art hardware platforms can be used for the implementation of the methodology.
The results in this article show that although an interface tracking mesh strategy is used, contact and topology changes can be modeled. In fact, in this framework these are to some degree equivalent. In case the mesh deformation is too extreme, there is always the option to use local or global remeshing, followed by a projection of the fields to the new mesh.
Note that the choice of constitutive models and finite element approximation spaces is not limited by the methodology, but is guided by the target application, in this article heart valve simulation. For thin structures, the choice of a piecewise linear finite element approximation of the incompressible solid mechanics needs to be used with care, due to potential shear-locking effects. In this article our focus was the FSI coupling and the contact model with topology changes, but in future work we will investigate the sensitivity of the leaflet model with respect to mesh resolution and shear-locking effects. Further, a more realistic material model of the leaflets will be studied in future work, and we will also investigate how to improve the modeling of the commissure, to avoid artificial stiffening of the leaflets, where they attach to the aorta.
The effect of mesh resolution on the fluid dynamics is also of great interest, and something that we will study in future work. In the present article the laminar boundary layers on the leaflets may be locally underresolved, and increased mesh resolution of the shear layers forming at separation from the leaflets may reveal more detailed turbulent structures. To optimize the computational resources, we plan to investigate local mesh resolution in the context of adaptive strategies.