On the Impact of Fluid Structure Interaction in Blood Flow Simulations: Stenotic Coronary Artery Benchmark

We study the impact of using fluid-structure interactions (FSI) to simulate blood flow in a large stenosed artery. We compare typical flow configurations using Navier-Stokes in a rigid geometry setting to a fully coupled FSI model. The relevance of vascular elasticity is investigated with respect to several questions of clinical importance. Namely, we study the effect of using FSI on the wall shear stress distribution, on the Fractional Flow Reserve and on the damping effect of a stenosis on the pressure amplitude during the pulsatory cycle. The coupled problem is described in a monolithic variational formulation based on Arbitrary Lagrangian Eulerian (ALE) coordinates. For comparison, we perform pure Navier-Stokes simulations on a prestressed geometry to give a good matching of both configurations. A series of numerical simulations that cover important hemodynamical factors are presented and discussed.


I Introduction
The application of computational fluid dynamics (CFD) to blood flow is a rapidly growing field of biomedical and mathematical research. Currently the development of numerical methods in hemodynamics is evolving from a purely academic tool to aiding in clinical decision making [36,9]. In particular the investigation of blood flow in stenosed arteries can help to shape medical treatment. For example virtual/computed Fractional Flow Reserve (cFFR), can evaluate the physiological significance of sclerotic plaque [25,8]. Moreover, the correct reconstruction of wall shear stress (WSS) is of crucial importance for the cell signaling and in consequence for the stenosis development [23] or for the assessment of rupture of brain aneurysms [12]. These two factors are studied in detail.
Hemodynamical simulations face numerous challenges, that are connected with measurements, medical image segmentation, mathematical modelling and development of numerical methods [29]. In this work we confine to an idealized geometry and a Newtonian fluid and focus on the comparison between a compliant wall vs. a rigid vessel wall. As a model example we have chosen a curved channel that resembles a large human artery. Stenotic and non-stenotic configurations are considered. The channel is preloaded by first considering a steady inflow to reach a certain physiological pressure and diameter before starting the pulsatile hear cycle. We investigate the aforementioned important clinical hemodynamical factors cFFR and WSS.
In the physiology of vessel macrocirculation the compliance plays a crucial role. However for individual arteries the problem is still open and there is no gold standard that can be adapted into clinical practice. For a discussion on assumptions in hemodynamic modeling of large arteries we refer to [35].
Numerous numerical studies were performed to compare rigid with compliant vessel simulations. The results report reduction of WSS for compliant vessels compared to rigid walls. In the case of a carotid bifurcation the authors of [27] could not observe significant changes in the flow patterns. In contrast relatively large vessels can undergo substantial radial wall motion, for example in the case of coronary arteries. The motion consists of bulk deformation and wall compliance that results in notable changes of flow characteristics [21,38].
After the brief introduction in Section I we present the monolithic formulation of the FSI problem, see Section II. The setting of the simulations is the topic of Section III. The numerical method is described in Section IV. Section V is dedicated to the presentation of investigated hemodynamical factors and the discussion of numerical results obtained on relevant test cases. Finally in Section VI we present a conclusion of this work.

II Model description
We consider a 3-dimensional domain Ω ⊂ R 3 , that represents a part of a vessel. The domain is partitioned in the reference configuration Ω = F ∪ I ∪ S, where F is the fluid domain, S the solid domain and I = ∂F ∩ ∂S is the fluid structure interface.
The velocity field v and the deformation field u are split into fluid v f := v| F , u f := u| F and solid v s := v| S , u s := u| S counterparts respectively. The pressure variable p f only exists on the fluid domain.
The boundary of the fluid domain Γ f := ∂F \ I is split into the inflow boundary Γ in f and the outflow boundary Γ out f . Similarly the solid boundary Γ s = ∂S \ I is split into inflow Γ in s and outflow Γ out s boundaries. Boundary conditions are described in Section III.

II.1 Fluid material model
Although, blood exhibits many unique properties, i.e. in certain regimes blood shows a non-Newtonian behaviour, we confine this work to considering an incompressible Newtonian fluid with the viscosity of µ f = 0.033 g · cm −1 s −1 and the density f = 1 g · cm −3 . The assumption of a Newtonian model for blood rheology is widely accepted for large and medium vessels, see e.g. [29]. The flow is governed by the Navier-Stokes equations

II.2 Solid material model
Arterial walls consist of heterogeneous layers with significant difference in physical properties. The principle layer construction of arterial wall consist of intima (inner layer), media (middle layer), and adventitia (outer layer). For detail account we refer to [16] and [20]. We briefly describe how the elastic constitutional law used in this work is derived.
Since arteries hardly change their volume within the physiological range of deformation [10], they can be regarded as incompressible or nearly incompressible materials. This motivates the application of a multiplicative decomposition of the deformation tensor (F) into its volumetric part J 1 3 and the deviatoric partF: Its associated modified deviatoric Cauchy Green tensorC then has the struc-tureC Thereby the free-energy function Ψ can be split in a volumetric and deviatoric part as described in [20] or [18]: Artery and vein walls consist of elastin and colagen fibres. The measurements of stress-strain curve exhibit stiffening effects at higher pressures due to the collagen fibres, c.f. [16]. Whereas under low loading of the artery the properties of elastin dominates. This motivates to model the artery as pseudoelastic material. Following [13] and [20] we employ an exponential deviatoric energy functional For the volumetric part of the energy functional the simple convex energy function is used as stated in [5,4]. By assuming hyperelastic stress response we obtain the second PiolaKirchhoff stress tensor Σ s : with the material parameters µ s = 44.2 kPa and γ = 20 as well as κ s = 4998 kPa. Similar values have been used in [3]. Moreover the solid density equalsˆ s = 1.2 g · cm −3 . The equation for the conservation of momentum in the elastic vessel wall is then given byˆ where we formulated the hyperbolic problem as a first order system in time by introducing the solid velocity v s and the deformation field u s as separate variables. The solid problem is formulated in the Lagrangian coordinates on the reference domain S, which is not moving with time. Hence, the densityˆ s = 1.2 g · cm −3 is the reference density which is unaffected by any compression or extension.

II.3 Fluid-structure interactions
One aim of this work is to study the impact of coupled fluid-structure interactions on typical blood flow configurations found in medical applications.
We therefore couple the Navier-Stokes equations (1) with the elastic material law (7) via coupling conditions on the common interface I. The coupled dynamics of a fluid-structure interaction problem leads to a free boundary problem with moving domains and a moving interface.
In the following, we briefly sketch the derivation of the coupled fluidstructure interaction problem. For details we refer to the literature [32,Chapter 5]. To overcome the mismatch of a Navier-Stokes equations (1) defined in Eulerian coordinates, e.g. on the evolving fluid domain F(t), and the solid problem (7) derived on the fixed reference domain S we use the well established concept of Arbitrary Lagrangien Eulerian (ALE) coordinates. See [14] or [32,Chapter 5] for a detailed derivation. We denote by F the fixed fluid reference domain and by T f (t) : F → F(t) the ALE map. Then, the velocity and the pressure can be mapped onto the fixed domain by definingv f : f . This allows to transform the Navier-Stokes problem in its variational formulation onto the fluid reference domain where φ and ξ are test functions. We denote by F f :=∇T f the gradient of the deformation variable and by J f := det(F f ) its determinant. Most characteristic feature of the ALE formulation is the appearance of the domain v that takes care of the implicit motion of the fluid domain. Furthermore, the Cauchy stress tensor is mapped to the reference domain which gives rise to the Piola transform J fσf F −T f . The reference stress is given bŷ It remains to describe the construction of the ALE map. Typically the ALE map is defined by means of an artificial fluid domain deformation u f via whereû f is an extension of the solid deformation u s from S to the fluid reference domain F. The most simple choice for the extension operator is to use a harmonic extension by implicitly solving the vector Laplacian For a discussion of this extension operator we refer to the literature [32, Sections 3.5.1 and 5.3.5]. Hereby, T f can be considered as a natural extension of the Lagrange-Euler map T s (x, t) := x + u s (x, t) such that we will skip the subscripts f and s when denoting the deformation T (x, t) := x + u(x, t), its gradient F = ∇T and its determinant J = det(F). In the following we skip all hats referring to the use of ALE coordinates.
As the fluid reference domain F and the Lagrangian solid domain S do not move, they always share the well defined common interface I. Here, we require the continuity of velocities, which is denoted by the kinematic coupling continuity of normal stresses, denoted as dynamic coupling condition and finally, the geometric coupling condition which says that the evolving domains F(t) and S(t) may not overlap and may not separate at the interface.
Since the velocity field and the deformation field are continuous across the interface, we formulate the coupled fluid-structure interaction problem using global solution fields v ∈ H 1 (Ω) 3 and u ∈ H 1 (Ω) 3 . Hereby, the kinematic coupling condition and the extension condition u f = u s in (9) are strongly realized as parts of the function spaces. The dynamic coupling condition is realized by testing the variational formulations of the Navier-Stokes equations an the solid problem by one common continuous test function φ ∈ H 1 (Ω) 3 . The resulting variational system of equations is given by For detailed account of monolithic formulations for fluid structure interactions we refer to [32].

III Simulation setup
We perform pulsatile blood-flow simulations by numerically solving system (10). The simulation setup that includes geometry and material parameters has been inspired by the Benchmark paper [3]. Note however that in some of the cases the assumption of rigid vessel walls is employed such that u = 0 and F = I on Ω. Then system (10) reduces to the standard incompressible Navier-Stokes equations. In what follows we distinguish between FSI-and NS-cases, respectively.

III.1 Geometry
The geometry of the computational domain reflects an idealized coronary artery. We show a sketch of the geometry in Figure 1. It consists of three parts, two straight and one curved tube. The straight parts are aligned with the x− and y− axes for inflow and outflow, respectively. The centerline of the curved section is a part of a circle with center in (1, 0, 0) and radius R = 1. It can be indicated as a parametrization ψ C : The fluid domain F is a cylinder around the centerline ψ C (s) with radius r F = 0.15 cm. Furthermore, the solid domain (i.e. the vessel wall) is a cylindrical outer layer of width 0.06 cm, i.e. the outer radius is given by r S = 0.21 cm. These dimensions correspond the dimensions of realistic arteries. We consider two stenotic and one non-stenotic configuration. The stenosis is modelled by a reduction of the inner radius r F on the second part of the curved boundary. This can be described by the parametrization Moreover we proceed twofold. First, we consider a symmetric stenosis such that the fluid domain is a cylinder of radius r F (s) around the centerline ψ C (s). Second, we modify the centerline in such a way that the stenosis is only on the inner side of the curve. This is achieved by shifting the centerline to match the radius variation, i.e. where χω(s) is the direction of the shift In both stenotic configurations, the resulting reduction of the radius equals 0.04 cm and the lumen area shrinks by 46% from 0.15 2 π cm 2 ≈ 0.07 cm 2 to 0.11 2 π cm 2 ≈ 0.038 cm 2 . Computational geometries are summarized in Table 1. Figure 1 shows Geometry 3.

III.2 Boundary conditions
The blood flow is enforced by a Dirichlet inflow condition on the inflow boundary Γ in f given by We prescribe a parabolic inflow profile where the maximum value v max (t) is presented in Figure 2. The temporal inflow profile v max (t) = 1 28.3 is an approximation, by means of a Fourier series of order 8, of a typical coronary velocity profile provided by [1]. The inflow profile consists of three stages: -Ramp phase (0 s ≤ t ≤ 0.1 s) with increasing inflow rate.
For the outflow of the fluid on Γ out f in the FSI case we apply absorbing boundary conditions: where p * = E 1−ν 2 . The condition is equipped with the reference pressure p e which is chosen such that a vascular pressure of 80 mmHg at t = 0.3 s is achieved. The absorbing condition explicitly disallows pressure waves to reenter the domain, for details see [26].
The structure is fixed at the inflow boundary and is allowed to move in y-z direction at the outflow boundary.
In order to compare both rigid and elastic computations, we run the FSI simulation first and then prescribe the resulting reference pressure profile P ref (t) at the outflow boundary in the rigid wall case such that This is the usual do-nothing outflow condition including a pressure offset, see [19].
The first two phases of (11) are designed to reach the intravascular pressure of 80 mmHg. This procedure resembles prestressing of the vessel wall. We first start the FSI simulation on the undeformed stress free reference domain and prestress the geometry by the ramp and steady phases. The resulting geometry at t = 0.3s is then extracted to compute the FSI as well as the Navier-Stokes simulation on this deformed geometry. This geometry corresponds to a reconstructed artery geometry from MRI or CT images, as during the measurements the blood flows with the pressure of at least 80 mmHg through the blood vessel. As the stress free reference domain is known here a priori, prestressing procedures as discussed in [17] do not need to be applied.

IV Numerical approximation, solution and implementation
The realization of a numerical framework for monolithic fluid-structure interactions is very challenging and a detailed description is not possible in one manuscript. We therefore give a brief description and refer to the relevant and detailed literature. See [32] for a comprehensive overview.
We employ a very strict form of the ALE formulation, where all equations are solved on the reference domain, no mesh update is used. This prevents the necessity of projections between moving meshes and also it is the only straightforward approach for obtaining higher order discretization in time [34]. In principle this approach allows for direct Galerkin discretizations of the monolithic variational formulation (10) and the choice of finite element spaces should be based on the following considerations -The finite element mesh should resolve the fluid-structure interface I in reference framework. -In the fluid domain, the velocity-pressure pair V h × Q h should fulfils the inf-sup condition to cope with the incompressibility constraint. Or, if a non-stable finite element pair is used, stabilization terms must be added. Our realization is based on equal order tri-quadratic elements for pressure and velocity, enriched with the local projection stabilization method [6]. -To reach a balanced approximation of the velocity-deformation pairing the same function space is used for the global deformation field.
To summon up, the discrete solution In time we use a variant of the Crank-Nicolson time discretization scheme that gives better stability properties by an implicit shifting. We refer to [34] for details. By restricting (10) to the fully discrete setting, a system on nonlinear algebraic equations arises in each time step. As nonlinear solver we employ a Newton scheme with an analytical evaluation of the Jacobian, see [30] or [32,Section 5.2.2] for details on the derivation.
The resulting linear systems of equations are very large and extremely ill-conditioned with condition numbers that are by far larger than those of fluid and solid equation on their own, see [31,2] for numerical studies. The approximation of these systems is still a great challenge, in particular if it comes to 3d applications. Only few fast solvers for the nonlinear setting are available [31,2,22]. Our approach is based on a multigrid solution that appears to be superior in 3d. In [15] we present the solution approach, which is based on a partitioning of the Jacobain based on two simple strategies: 1. Within the Navier-Stokes equations we neglect those parts of the Jacobian that come from the derivative with respect to the domain extension u f . The resulting nonlinear solver is an approximated Newton method that however still solves the original problem since the residual is exact. In [32,Section 5.2.3] and [15] we have found that such an approximated Newton solver is even more efficient, despite slightly increased iteration counts. This is due to the very costly evaluation of the full Jacobian that is not required in our approach. ) .
The combination of these two modifications allows for a natural splitting of the Jacobian within the multigrid smoother. Further, it allows to apply very simple iterations of Vanka-type, as smoother in the geometric multigrid preconditioner, that are easy to parallelize. In [15] we have demonstrated the efficiency of the approach in different 3d configurations. The implementation is based on the finite element toolkit Gascoigne3D [7].

V Simulations
We finally use the presented finite element framework for a computational analysis of several hemodynamic parameters that are relevant to answer clinical questions. We focus on the dependence of these parameters on the elasticity of the vessel walls. The additional effort of fully coupled fluid-structure interactions over pure Navier-Stokes simulations is immense and should be justified by corresponding effects. The behavior of three specific parameters is investigated. The wall shear stress (WSS) plays an eminent role in several applications. It is used as indicator to model the growth of atherosclerotic plaques [39,33] but also when it comes to deriving measures to evaluate the risk of plaque rupture [24] or the rupture of aneurysms [12]. Both the minimum wall shear stress and the distribution of the wall shear stress on the vessel walls are important measures. In Section V.1 analyze the effect of the different complexities, Navier-Stokes and fluid-structure interactions, on the WSS distribution in all three different geometries.
Furthermore we analyze the computational fractional flow reserve (cFFR). The fractional flow reserve (FFR) is a technique that measures pressure differences across a stenosis. In medical practice, the FFR is determined by inserting a catheter in the artery and measuring flow parameters at maximum blood flow. During the procedure, the tip of the catheter, where the sensor is located, is retrieved, such that measurements are available along the affected section of the vessel. Healthy vessels should give a pressure ratio close to 1. If the ratio drops below 0.8, i.e. a 20% drop in pressure, the stenosis is considered to be severe [37]. Aim of the computational fractional flow reserve is to replace the risky intervention by computer simulations based on medical imaging.
Finally, we discuss the amplitude of the pressure oscillation during one heart cycle. In clinical observation is usually observed that the amplitude significantly drops after a stenotic region in a blood vessel. By comparing Navier-Stokes simulations with coupled fluid-structure interactions, we will show that this effect can only be described by considering the fully coupled model including elastic vessel walls.

V.1 Wall shear stress (WSS)
The tangential component of the surface force at the vessel wall is denoted as wall shear stress (WSS). By means of the Cauchy's theorem we have where n is a unit outward normal vector to the vessel wall and τ is corresponding tangential vector. Note that WSS is a vector and its often confused with its magnitude |WSS| 2 which is a scalar quantity denoted by wss.
The plots of the wall shear stress are presented on the boundary of the fluid domain, see Figure 3, Figure 5 and Figure 4. The surrounding solid is removed from these plots such that only the fluid domain is given.
We use the same scale ranging from 0 Pa to 20 Pa in all figures. The distribution of wss is always presented for 3 different points in time. First, when the pulsatile flow reaches its minimum inflow pressure, then, at maximum pressure and finally for an intermediate pressure value. Each specific situation is shown from two different angles, such that the inner and outer surfaces are well visible.
Comparing the coupled fluid-structure interaction model (left) with the pure Navier-Stokes flow (right) always shows much higher wss values in the Navier-Stokes case. Regions of very high wss are only found in the outer surface in the case of Navier-Stokes. In the case of the FSI model, the values are smoothly spread. In general, the vessels are widened and show a larger diameter of the lumen in the case of FSI. Prestressing of the domains for the Navier-Stokes simulations is based on the deformation resulting from the increasing inflow of the ramp phase (11) to reach pressure of 80mmHg. Figures 5 and 4 show a shift of the wall shear stress distribution. High values are found all around the stenosed parts, on the inner and the outer surfaces of the curved areas. Again, the Navier-Stokes model is not able to yield a proper distribution of wss around the cylinder surface and it only concentrates in the outer surface of the curved section. To sum up, it appears to be important to consider elastic fluid-structure interactions, if the spatial distribution of the wall shear stress is of interest. The Navier-Stokes model is nearly unaware of the stenosis and always concentrates the WSS on the outer surface of the curved region. If rupture locations [12] or localized growth processes are of interest [39,33], the use of FSI is essential.

V.2 Computational fractional flow reserve (cFFR)
The computational fractional flow reserve (cFFR) is the ratio of the maximum blood pressure after the stenosis (distal) and before the stenosis (proximal). We denote by p d the distal pressure and by p a the proximal pressure such that the computational fractional flow reserve is defined by In our configuration, compare with Figure 1, we evaluate the distal and proximal values in p a = p(0, 0, 0) and p d = p(1.7, 1, 0).
Since the cFFR varies slightly with time the maximum value over the cardiac cycle is computed. Further computational aspects of cFFR are covered by the recent benchmarking paper [11]. In healthy vessels we expect cFFR ≈ 1 and in the medical practice, a stenosis with cFFR > 0.8 is considered to be functionally non significant [37]. No medical intervention, e.g. by placing a stent would be required.
The results of cFFR computations are presented in Table 2. The first striking observation is that the Navier-Stokes model is not able to reflect the presence of the stenosis at all. A pressure drop of about 5% is observed for all configurations and is only due to the curvature of the domain. The FSI model is well able to yield cFFR ≈ 1 in the case of the healthy configuration and shows a loss of about 5% in both stenotic geometries. Based on these simulations, the stenosis would not be considered to be severe and in the need of an intervention. These results clearly show that a pure Navier-Stokes simulation is not able to serve as computational basis for replacing the medical FFR procedure by simulations.

V.3 Pressure amplitude
The third quantity of interest is the dynamic pressure amplitude before and after the stenosis, i.e. a temporally resolved analysis of the proximal and distal pressures p a and p d evaluated in the coordinates as mentioned in (14). We study the progress of the pressure to investigate possible damping effects of a stenosis. The change of pressure profile after strongly stenotic regions, is reported as an assessment tool for cardiovascular risk [28].
We show the evolution of the pressures p a (t) and p d (t) over the third heart beat in Figure 6. On the left, we show results for the fully coupled FSI models, on the right we give the corresponding pressure lines for the Navier-Stokes case. Comparable to the cFFR study given in the previous section, the most striking observation is the invariance of the Navier-Stokes solution to the kind of vessel geometry. For all cases, healthy vessel, centered stenosis and non symmetric stenosis, the Navier-Stokes results are nearly identical and always indicate a loss in pressure amplitude of about 7.5%. In contrast, the FSI solution is able to preserve the oscillation in the case of the healthy vessel and shows a drop of about 5% for the two stenotic cases. A closer look reveals that the distal pressure lines are very similar in all cases. The slight oscillations in the distal pressure in case of the FSI simulation are not numerical instabilities. Instead they stem from the oscillatory behavior of the elastic vessel wall.
Finally, Table 3 collects the amplitudes before and after the stenosis. This third study also shows enormous qualitative and quantitative discrepancies between the simulation results depending on the model under consideration.

VI Conclusions
We have studied a prototypical geometry that represents a large and curved artery. Three different configurations are considered: a healthy vessel with a diameter (the lumen) of 0.15 cm and two cases where a stenosis is given in the curved area of the vessel. First, the stenosis is centered around the vessel centerline, finally, a shifted variant where the stenosis is concentrated on the inner surface. The blood flow is driven by a time-dependent flow rate with values that are clinically relevant. For all configurations we study the effect of the elasticity in the vessel walls, i.e. we compare a pure Navier-Stokes simulation with a fully coupled fluid-structure interaction system.
Three different indicators that are also relevant in clinical decision making are investigated: the distribution of the wall shear stress that is made responsible for stenosis growth and risk of rupture, the computational fractional flow reserve that is used to estimate the severity of a stenosis and the amplitude of the pressure oscillation that also measures the severity of a stenosis. In all cases we observe that the simple Navier-Stokes model is not able to depict the effect of the plaque. In particular the pressure lines are nearly identical for all three geometrical configurations. The FSI model is however well able to replicate clinical observations. For instance, the energy, in terms of pressure oscillations, is fully preserved throughout the curved region, if elasticity of the vessel walls is taken into account.
Although the study covers only an idealised geometry and boundary conditions, the proposed regime corresponds to physiological blood flow. Therefore, we stress that the compliance of the vessel has significant impact on clinical hemodynamical factors. In medical practise one has to be aware of the difference in the results from CFD with FSI and NS models and the limitations of the later.