1 Introduction

Just in the United States more than 100,000 heart valve surgeries are performed annually on patients with heart valve disease to replace native valves with prosthetic heart valves [1]. Approximately half of the replaced valves are bileaflet mechanical heart valves (BMHV) (Fig. 1), which are particularly attractive to relatively young patients due to their durability and long life span as they are usually made from pyrolytic-carbon. However, all current BMHV designs are prone to a number of complications including among others increased risk of thromboembolism and possible hemorrhage [64]. Due to the propensity of BMHV recipients to thromboembolic complications, patients are required to take anti-coagulant medication which if poorly managed can lead to life-threatening hemorrhage. These complications are believed to be strongly associated with the complex, non-physiologic blood flow patterns induced by BMHVs [21, 63, 64] and for that ongoing research efforts focus on understanding BMHV hemodynamics at physiologic conditions and quantifying the link between the hemodynamic environment and the potential for thromboembolic complications [21].

Fig. 1
figure 1

A typical clinical quality bi-leaflet mechanical heart valve (St. Jude Regent)

Full size image

Numerous in vivo and in vitro experimental studies have been carried out to better understand the flow patterns and the mechanical environment induced by BMHV [10, 17, 20, 44, 45, 55, 65]. Many of these experiments use state-of-the-art particle image velocitometry (PIV) to measure the instantaneous and phase-averaged flow patterns downstream of the valve leaflets [6, 10, 38, 39]. It is important to recognize, however, that PIV experiments, regardless of their sophistication and resolution, can only provide 2D cross-sections through the very complex flow field induced by the valve leaflets as they interact with the incoming pulsatile flow. High-resolution 3D measurements are very challenging if not impossible to obtain with PIV especially in in vivo studies.

Numerical simulations provide the only viable tool for quantifying the BMHV-induced hemodynamics at the level of detail and resolution required for establishing the link between the mechanical environment experienced by blood cells and the potential for thromboembolic complications. The numerical simulation of such flows, however, poses a major challenge to even the most advanced computational techniques available today. BMHV flows at physiologic conditions take place in complex geometries with compliant walls, involve geometric features at disparate spatial scales (from the scale of the aorta diameter to the scale of the valve hinges and leakage jet), are dominated by fluid–structure interaction (FSI) and the pulsatile flow effects, and undergo transition to turbulence. It is important to note that most of BMHV simulations have been carried out in the aortic position rather than the mitral position. An important hemodynamic reason for this choice is that the shear stresses are higher in the aortic position [64] and, thus, the potential of BMHV induced blood cell damage is greater in the aortic position. Moreover, the geometry/movement of the aorta is less complicated relative to the left ventricle downstream of the mitral valve so simulations in the aortic position are less challenging from the computational standpoint. Pulmonic and tricuspid valve flows have been found to be similar to aortic and mitral valve flows, respectively, but with lower overall velocity magnitudes [64]. Therefore, this review mainly focuses on the BMHV flows in the aortic position.

In this review paper we: in Sect. 2 present an overview of the computational methods that have been proposed in the literature for simulating BMHV flows and summarize major recent algorithmic advances that have permitted the first, high-resolution FSI simulations of BMHV flows at physiologic conditions; in Sect. 3 discuss recent simulation results and the insights gained through them about the BMHV flow physics; and in Sect. 4 discuss the next frontiers toward developing patient-specific computational tools for studying BMHV hemodynamics in vivo.

2 Governing equations and numerical methods

2.1 Fluid equations

The governing equations for the blood flow through mechanical heart valves are the 3D, unsteady incompressible continuity and Navier–Stokes equations, which in compact tensor notation read as follows:

$$ {\frac{\partial u_{i}} {\partial x_{i}}} =0 $$
(1)
$$ {\frac{du_{i}} {dt}} =-{\frac{\partial p} {\partial x_{i}}}+{\frac{1} {Re}}\;{\frac{ \partial ^{2}u_{i}} {\partial x_{j}\partial x_{j}}} $$
(2)

where u i are the Cartesian velocity components, p the pressure divided by the density ρ, and Re the Reynolds number of the flow based on a characteristic length and velocity scale. d/dt is the material derivative defined as:

$$ {\frac{\rm d} {{\rm d}t}}(\cdot )={\frac{\partial } {\partial t}}(\cdot )+u_{j}{\frac{\partial } {\partial x_{j}}}(\cdot ) $$
(3)

In Eq. 2 blood has been assumed to be Newtonian. For the purpose of this paper, which deals exclusively with flows in large arteries (e.g. aorta), this assumption is valid. It is important to keep in mind, however, that there are small regions of the flow domain (valve hinges and leakage jet during closure) within which non-Newtonian effects could become important and should be taken into account in the governing equations even for a BMHV implanted in the aortic position. Modeling of these fine, albeit potentially important from a hemodynamic standpoint, flow features is beyond the scope of this review and will not be discussed herein.

Equation 2 need to be solved in a domain defined by the aortic lumen and the left ventricle within which the BMHV leaflets and valve mounting mechanism are immersed. In the in vivo setting, the aortic and left ventricle walls are compliant and deform dynamically within the cardiac cycle. The valve leaflets are also free to pivot around their hinges and open and close periodically in response to the cardiac flow pulse. Therefore, Eq. 2 need to be solved in a domain enclosed by a moving boundary that contains moving immersed boundaries undergoing arbitrarily large deformations and boundary conditions need to be satisfied at the interfaces between the blood and the aortic wall and the blood and the valve leaflets. In what follows, we present a brief review of numerical methods for handling such complex FSI problem. It is important to note that in this paper we will simplify somewhat the problem by assuming that the aortic wall is rigid and only the BMHV leaflets are free to move. In addition, the motion of the left ventricle, which creates the physiologic pulse that drives the blood flow through the aorta, has been replaced by a physiologic inflow waveform with a plug flow profile prescribed at the inlet of the computational domain several diameters upstream of the BMHV (see Fig. 2 for the geometry of the simplified problem and Fig. 3 for a typical inflow waveform). The methods that will be reviewed below, however, are, at least in principle, suitable for handling the full FSI problem involving the compliant aortic/left ventricle walls and the moving valve leaflets. Such methods can be broadly classified as: (1) moving grid methods; and (2) fixed grid methods.

Fig. 2
figure 2

The 3D BMHV geometry (23 mm St. Jude Regent) including the housing and the leaflets in a straight aorta (left) and the definition of leaflet angle (right). Taken from [4]

Full size image
Fig. 3
figure 3

BMHV simulation [4]. Physiologic inflow waveform (dashed line) and comparison of the calculated leaflet kinematics (solid line) with experimental observations [10] (circles). Adopted from [4]

Full size image

Moving grid methods are methods in which the computational grid is fitted to and moves/deforms with the moving boundary. The movement of the grid is taken into account by using the arbitrary Lagrangian Eulerian (ALE) formulation of the governing equations [15]. This method has been previously applied to simulate the flow through mechanical heart valves [7, 8]. However, a significant restriction of the ALE approach stems from the fact that the mesh must conform to the moving boundary at all times and as such it needs to be constantly displaced and deformed following the motion of the boundary. Updating the mesh at every time step could be, however, quite challenging and expensive especially for complicated 3D problems. Cheng et al. [8], for instance, had to use interpolation between two previously generated meshes to obtain the intermediate mesh for a given leaflet angle and then applied an elliptic solver to the entire mesh to smooth it. The difficulties with ALE methods are further exacerbated in problems involving large structural displacements, as is the case with the BMHV leaflets. In such cases, obtaining smooth and well-conditioned computational meshes at every time step is far from trivial if not impossible and frequent remeshing may be the only option. Because of these inherent difficulties, the ALE approach is not the best choice for simulating BMHV flows, which are geometrically complex and involve large structural displacements.

Fixed grid methods are becoming increasingly popular in recent years due to their enormous versatility in simulating FSI problems involving large structural discplacements [4, 10, 13, 24, 54, 58]. In such methods the entire fluid computational domain is discretized with a single, fixed, non-boundary conforming grid system (most commonly a Cartesian mesh is used as the fixed background mesh) while the structural domain is discretized with a separate grid, which can move freely inside the fluid domain. The effect of a moving immersed body on the fluid is accounted for by adding, either explicitly or implicitly, body forces to the governing equations of fluid motion so that the presence of a no-slip boundary at the location of the solid/fluid interface can be felt by the surrounding flow. Since the grid used to discretize the fluid domain does not have to conform to the moving immersed boundaries, such methods are inherently applicable to moving boundary problems involving arbitrarily large structural displacements such as mechanical heart valves.

The earliest work to apply a fixed grid method to simulate heart valve flows is Peskin’s pioneering work with the immersed boundary (IB) method [48]. In this method the presence of the immersed deformable solid boundary on the surrounding fluid grid nodes is accounted for by adding a body force in the Navier–Stokes equations. The body force is distributed on all nodes of the fixed background grid via a discrete delta function that has the effect to smear, or diffuse, the solid boundary over several fluid grid nodes in the vicinity of the boundary. Because of this inherent property of the method, Peskin’s method, which has also been applied to simulate the flow in a complete heart model [49], is known as a diffused interface method. The original IB method is only first order accurate [48] but a variant of the method that is formally second-order accurate and combines adaptive mesh refinement to increase resolution in the vicinity of immersed boundaries has also been proposed [28].

The fictitious domain method [27] is another related fixed grid, diffused interface method that has also been applied to heart valve simulations [11–14, 58–60]. In this method the immersed solid is, as in Peskin’s IB method, free to move within the fluid mesh but the two domains are coupled together at the solid/fluid interface through a Lagrange multiplier (or local body force) [60]. The fictitious domain method had been applied to simulate flow in a 2D model of the native valve [12] as well as in a 3D trileaflet heart valve at relatively low, non-physiological Reynolds number (peak systole Re = 900) [11, 13, 14]. A major issue with the fictitious domain method is that it cannot yield accurate results for the viscous shear stresses on the solid boundary [60]. To remedy this major limitation, a combination of the fictitious domain method with adaptive mesh refinement has also been proposed [59, 60]. This enhanced fictitious domain method has been applied in 3D FSI simulations assuming, as in previous studies, that the geometric symmetries of the valve will also be respected by the induced flow field [58]. Due to the high computational costs of this approach, however, it has yet to be used to carry out a full 3D simulation of heart valve flows at physiologic conditions.

As mentioned above, both the IB and fictitious domain methods are diffused interface techniques causing the smearing of the solid/fluid interface. Because of this property, such methods typically require increased grid resolution in the vicinity of the boundary for accurate results and their application to high Reynolds number flows can be impractical from the computational standpoint. To remedy this situation, a class of sharp-interface immersed boundary methods has recently been developed and are attracting increasing attention in the literature—see for example [9, 18, 31, 35, 40, 56, 62] among many others and some recent publications from our group [4, 24–26]. In these methods the immersed boundary is treated as a sharp interface and its effect on the fluid can accounted for in a variety of ways. For example, in the cut-cell methods [57] the shape of grid cells in the vicinity of the boundary is modified to produce a locally boundary-fitted mesh. In immersed interface methods [35, 36, 62] the jump conditions, caused by the discrete delta function in the classical IB method, are incorporated into the finite difference scheme to avoid the approximation of the discrete delta function by a smooth function and eliminate the smearing of the interface. In the hybrid Cartesian/Immersed-Boundary (HCIB) method [25], developed by our group, boundary conditions are reconstructed in the vicinity of the immersed boundary via an appropriate interpolation scheme along the local normal to the boundary. Ge and Sotiropoulos [24] develop a novel paradigm integrating the HCIB method with body-fitted curvilinear grids. Their method, which was dubbed the Curvilinear Immersed Boundary (CURVIB) method, is ideally suited for simulating heart valve flows in anatomic geometries. The empty aorta is discretized with a boundary-conforming curvilinear mesh and the valve leaflets are treated as sharp, immersed boundaries within the background curvilinear mesh. A major issue for the efficient implementation of the CURVIB method, and for that matter for other sharp-interface methods, is the efficiency of the algorithm for classifying at every time step the location of the nodes of the background grid relative to the moving immersed boundary. This issue was successfully addressed by Borazjani et al. [4] who incorporated a new and very efficient algorithm in the CURVIB method to classify the fluid domain grid nodes into fluid, solid, and immersed boundary. For a detailed discussion of various sharp-interface methodologies the reader is referred to [25] and the recent review paper by Mittal and Iaccarino [41].

2.2 Leaflet equations

The motion of the leaflets of mechanical heart valves is governed by the angular momentum equation around the hinge axis, which after non-dimensionalization reads as follows:

$$ {\frac{\partial ^{2}\theta} {\partial t^{2}}}+\zeta {\frac{\partial \theta} {\partial t}}={\frac{1} {I_{\rm red}}}C_{\Im } $$
(4)

In the above equation, θ = θ(t) is the leaflet angle defined as in Fig. 2, which can vary between θmin and θmax. I red is the leaflet reduced moment of inertia defined as:

$$ I_{\rm red}={\frac{\rho_{\rm s}I_{xx}} {\rho _{\rm f}D^{5}}} $$
(5)

where, D is diameter of the aorta, ρf the density of the fluid, ρs the density of the leaflet material, and I xx  = ∫r 2dV (with \(r=\sqrt{(y-y_{\rm h})^{2}+(z-z_{\rm h})^{2}}; y_{\rm h}\) and z h position of the hinge axis; and dV as the infinitesimal volume element) is the product of inertia around the leaflet’s hinge axis. The damping coefficient is defined as

$$ \zeta ={\frac{cD} {I_{xx}U}} $$
(6)

where c is the damping factor due to hinge friction, and U the peak bulk flow velocity at inlet. In most studies the hinge friction has been neglected due to its small value relative to the flow forces and lack of experimental value in most of BMHV studies [4, 5, 16, 43, 61]. The moment coefficient exerted by the flow on the leaflet is defined as

$$ C_{\Im }={\frac{\Im} {\rho _{\rm f}U^{2}D^{3}}} $$
(7)

where, \(\Im\) is the moment around the hinge axis. Note that the gravity has been neglected since we have assumed that the hinge axis is in the direction of gravity. Therefore, gravity has no effect on the moment around the hinge axis.

Equations 4 comprise a system of second-order ordinary differential equations. Numerically these equations are typically solved by first transforming them into a system of first-order ordinary differential equations as follows:

$$ {\frac{\partial \theta} {\partial t}} = \omega $$
(8)
$$ {\frac{\partial \omega} {\partial t}}={\frac{{C}_{\Im}} {I_{\rm red}}}- \zeta \omega $$
(9)

2.3 Boundary conditions and coupling for FSI problems

The fluid and structure dynamics are coupled together at the fluid/structure interface Γ by the following boundary conditions:

$$ {\mathbf{u}}={\mathbf{U}}={\mathbf{r}}\times {\mathbf{n}}{\frac{\partial \theta} {\partial t}} \quad {\rm at} \quad \Gamma $$
(10)

where, \({\mathbf{u}}\) and \({\mathbf{U}}\) is the velocity of fluid and solid at the interface, respectively. \({\mathbf{r}}\) the position vector from the hinge and \({\mathbf{n}}\) is the normal vector at the fluid/solid interface Γ. Equation 10 couples the Eulerian velocity field of the fluid with the Lagrangian description of the motion of the solid surface. Other than the kinematic boundary condition Eq. 10 the fluid and structure domains are also coupled together at the interface by the dynamic boundary condition, i.e. the force exerted on the solid is equal but opposite the direction the force exerted on the fluid, which shows itself on the right hand side of Eq. 4. A practical approach to simulate the FSI, comprising of Eqs. 2, 4, and 10, is the partitioned approach [19]. In this approach the FSI problem is partitioned into two separate fluid and structure domains. Each domain is treated computationally as an isolated entity and is separately advanced in time. The interaction effects are accounted for through boundary conditions at the interface [4, 19].

The partitioned approach can be implemented either in a loose or strong coupling fashion [4]. The domains are loosely coupled (LC-FSI) if the boundary conditions at the interface are obtained from the domain solutions at the previous time level (explicit in time) i.e. the flow-imparted moment \(C_{\Im}\) in the right hand side of equations (4) is evaluated from the previous time level flow field solution. They are strongly coupled (SC-FSI) if the interfacial boundary conditions are obtained from the domain solutions at the current time level (implicit in time). This is achieved by performing several sub-iterations per physical time step until the FSI equations (Eqs. 2, 4 with the kinematic boundary condition equation 10) have converged within a desired tolerance at the n + 1 time level. Assuming that the solutions for both the fluid \({\mathbf{q}}^{n}=( p^{n}\quad {\mathbf{u}}^{n}) ^{T}\) and structure θn domains are known at time levels n and n−1, the strong-coupling algorithm determines the solution at the next time step n + 1 as follows [4]:

  1. (1)

    Set \(\widetilde{\theta}^{1}\equiv \theta^{n}\) and \(\widetilde{\mathbf{q}}^{1}\equiv {\mathbf{q}}^{n}.\) Starting from ℓ = 1, loop over ℓ (steps a to e) until convergence is achieved:

    1. (a)

      Using the known position of the structure \(\widetilde{\theta}^{\ell }\) to prescribe boundary conditions for the fluid domain, solve the Navier–Stokes Eq. 2 to obtain \(\widetilde{{\mathbf{u}}}^{\ell }\) and \(\tilde{p}^{\ell }.\)

    2. (b)

      Calculate the force on the structure domain as \(\widetilde{C_{\Im}}^{\ell }=C_{\Im}(\widetilde{\theta}^{\ell },\widetilde{{\mathbf{q}}}^{\ell })\)

    3. (c)

      Solve the structure domain Eq. 9 with the trapezoidal rule as

      $$ {\frac{\widetilde{\omega}^{\ell+1}-\omega^{n}} {\Updelta t}} ={\frac{\widetilde{ C_{\Im}}^{\ell}+{C_{\Im}}^{n}} {2}} $$
      (11)
      $$ {\frac{\widetilde{\theta}^{\ell +1}-\theta^{n}} {\Updelta t}}= {\frac{\widetilde{ \omega}^{\ell +1}+{\omega}^{n}} {2}} $$
      (12)
    4. (d)

      Check for convergence of the structure solution:

      $$ ||\widetilde{\omega}^{\ell +1}-\widetilde{\omega}^{\ell }|| < \varepsilon $$
      (13)

      where ɛ is a preset convergence ratio set equal to ɛ = 10−6 and || . || is the infinity norm in all our simulations.

    5. (e)

      If not converged, increment ℓ by one and return to step a to continue the iteration loop. If converged, then end the loop and go to 2.

  2. (2)

    Set \(\theta^{n+1}=\widetilde{\theta}^{\ell +1}, {\mathbf{u}}^{n+1}=\widetilde{{\mathbf{u}}}^{\ell +1}, p^{n+1}=\widetilde{p}^{\ell +1}\) and \(C_{\Im}^{n+1}=\widetilde{C_{\Im}}^{\ell +1}\)

The loose coupling is similar to the above strong coupling algorithm but only one sub-iteration (ℓ = 1) is performed. For a detailed discreption of the SC- and LC-FSI the reader is referred to Borazjani et al. [4].

2.4 Stability of the FSI coupling for MHV simulations

The SC-FSI and LC-FSI couplings have different stability and computational cost. LC-FSI is computationally attractive since the fluid and structure domains are solved only once per time step. LC-FSI, however, due to the explicit nature of time advancement is generally less stable than the SC-FSI. For a given FSI problem the stability of the coupling method depends on the level of interaction between the fluid and the structure solutions in the given problem. If the solution of one domain is very sensitive to small changes in the solution of the other domain then there is high interaction. As shown by Borazjani et al. [4] usually there is a high interaction when the added mass/inertia of the system is similar to the actual mass/inertia of the system. This is normally the case for mechanical heart valves, which are characterized by their low inertia of the leaflets, e.g. I red = 0.001 for a St. Jude Regent 23 mm BMHV [4]. In such cases, not only the LC-FSI but also the SC-FSI is not stable. To remedy this problem Borazjani et al. [4] used under-relaxation to stabilize the SC-FSI iterations. The convergence of the SC-FSI sub-iterations were found to be greatly dependent on the under-relaxation coefficient. The under-relaxation coefficient was dynamically evaluated using the Aitken acceleration technique [42], which greatly reduced the number of SC-FSI sub-iterations needed for convergence. The SC-FSI typically converged within 4–5 sub-iterations using the Aitken acceleration technique [4].

Borazjani et al. [4] used a theoretical analysis for a simple FSI problem to explain the findings of their numerical experiments regarding the stability of different FSI coupling methods for their BMHV simulations. They showed that the ratio of the added mass to the mass of the structure as well as the sign of the local time rate of change of the force or moment imparted on the structure by the fluid determine the stability and convergence of the FSI coupling. This explains the striking finding from their MHV simulations that the LC-FSI coupling is, as one would expect, unstable during the valve opening phase but the same coupling is stable and yields very similar solution to that obtained by the SC-FSI algorithm during the valve closing phase. This is particularly surprising when one takes into account the relative complexity of the instantaneous flow during the opening and closing phases—the flow being well organized and laminar during the opening phase and very complex, turbulent-like during valve closing. Their stability analysis clarified this seemingly paradoxical finding by showing that during the closing phase the sign of the rate of change of the flow-induced moment on the valve leaflets is such that flow effects tend to alleviate the adverse stability effects of the added mass term, which arise due to the very low inertia of the valve leaflets [4]. Furthermore, with their analysis the stabilizing role of under-relaxation is also clarified and an upper bound of the required for stability under-relaxation coefficient was derived [4].

3 Recent simulations and insights into the flow physics

The early simulations of the flow in BMHV did not considered moving leaflets and were carried out for fixed leaflets and steady inflow. Shim and Chang used a finite-element code to simulate the 3D flowfield in a tilting disk valve in the half open position [51] and a Björk–Shiley valve [52]. Kris et al. [33] used a finite-volume, artificial compressibility method with overset grids to solve the 3D Reynolds-averaged Navier–Stokes (RANS) equations closed with an algebraic turbulence model to simulate the flow through a Björk–Shiley tilting disk valve fixed at the 30° angle and Re ranging from 2,000 to 6,000. King et al. [32] used the commercial code finite-element FIDAP, which is based on on the Galerkin method of weighted residuals was used to model the flow through one quarter of the valve with fixed leaflet during the first half of the systole with peak Re = 1,500. Ge et al. [22, 23] carried out direct numerical simulations at Re = 750 and 6,000 with the leaflets in the fully opened position on fine girds with about 1.5 million grid nodes. These computations showed that the flow is highly 3D even with the fixed leaflets and questioned the validity of the computationally expedient assumption of flow symmetry.

Most of the recent work focused on moving leaflets and pulsatile flow. A considerable amount of this work has been carried out using 2D models of the actual BMHV problem. van Loon et al. [59, 60] and de Hart et al. [12] have performed full FSI simulations deformable slender structures in 2D, the "leaflets" of a native valve, with a fictitious domain finite element method. Stijnin et al. [53] have used fictitious domain method for 2D simulation of mechanical heart valves at peak systole Re = 750. Bluestein et al. [2] have performed 2D unsteady RANS simulations for fixed leaflet position to study platelet activation. Cheng et al. [7] and Krishnan et al. [34] carried out 2D FSI simulations for valve closure and Pedrizzetti and Domenichini [46, 47] during valve opening. Rosenfeld et al. [50] have performed 2D simulations of a tilting disk valve in the mitral position. A 3D FSI simulation has been carried out by Cheng et al. [8] using only one quarter of the valve (quadrant symmetry assumption) and a grid with about 200,000 nodes to discretize the flow domain. van Loon et al. [58] have also carried out 3D FSI simulations of a tissue valve with symmetry assumption. In spite of considerable progress made and many new insights into the problem gained by these studies, full 3D FSI simulations of BMHV in anatomically realistic aorta geometries, at physiologic conditions, and at resolution sufficiently fine to resolve hemodynamically relevant scales of motion were out of reach of modern computational methodologies up until very recently. Note for instance that de Hart et al. [11, 13, 14] have performed 3D simulations of tissue valves for a peak systole Reynolds number Re = 900, a value that is only a fraction of the actual physiologic range (Re = 5,000–6,000). The first attempts to simulate the flow at physiologic, pulsatile conditions were reported only recently. Tai et al. [54] used the immersed object method, which adds a body force to the momentum equations such that the desired velocity is obtained on the object boundary, with overlapping grids to carry out FSI simulations of bileaflet MHV on an unstructured mesh (86,000 nodes, 450,000 elements in four zones) with an artificial compressibility solver enhanced with multigrid. They preformed the simulations under physiologic conditions but the coarseness of the computational mesh did not allow them to obtain insights into the underlying physics of the flow. Nobili et al. [43] carried out FSI simulations of a BMHV using the commercial code Fluent on a relatively course mesh (1.2 million tetrahedral and 900,000 hexahedral). Even though the computational mesh used in this study was finer than previous FSI simulations, it was still far too coarse to capture the rich dynamics of the flow throughout the cardiac cycle as revealed by recent laboratory experiments [10]. Consequently, the simulations of Nobili et al. yielded a simulated flow environment downstream of the valve leaflets that was significantly simpler and less rich dynamically than observed in experiments.

The CURVIB method of Ge and Sotiropoulos [24] was the first method to be successfully applied to perform a highly resolved (10 million grid nodes to discretize the empty aorta) direct numerical simulations of physiologic pulsatile flow in a BMHV mounted in a straight aorta (Fig. 2) and to validate the numerical simulations with detailed laboratory experiments—see [24] for the details of the method and [10] for validation and discussion of BMHV flow physics. The simulations were carried out in a domain that extended four diameters upstream of the valve, where the physiologic wave form was prescribed, and ten diameters downstream of the leaflets where outflow conditions were imposed. The CURVIB method was shown to be able to reduce the discrete divergence of the velocity field to machine zero at all instants in the cardiac cycle. It should be noted, however, that both in [10] and [24] the simulations were carried out by prescribing the kinematics of the valve leaflets from experimental measurements and as such the FSI aspects of the BMHV problem were not considered. In [4] the CURVIB approach was extended to develop a coupled FSI formulation that is applicable to problems involving multiple, moving, rigid bodies of complex geometry undergoing arbitrarily large structural displacements. Borazjani et al. [4] reported the first full FSI high-resolution, direct numerical simulation of a BMHV under physiologic conditions (peak systole Re = 6,000), which could capture not only the leaflet kinematics (see Fig. 3) but also all the flow physics in excellent agreement with the experimental results (see Fig. 4).

Fig. 4
figure 4

BMHV simulation [4] compared with the PIV measurements of Dasi et al. [10]. Instantaneous out-of-plane vorticity contours on the mid-plane of the valve (a) from simulation [4] and (b) from experimental measurements [10]. The contour levels are identical. c Instantaneous vortical structures visualized by iso-surfaces of q-criteria. The dots on the inflow waveform shown at the bottom of each column indicate the time instant during the cycle for that column. Taken from [4]

Full size image

The CURVIB-FSI method [4] along with the experiments of Dasi et al. [10] provided the first comprehensive insights into the instantaneous hemodynamic environment induced by the BMHV leaflets at physiologic conditions and at hemodynamically relevant scales, at least for the in vitro configuration with a straight, axisymmetric, rigid-wall aorta. The following major conclusions regarding the physics of the flow emerged collectively from the high resolution experiments of Dasi et al. [10] and FSI simulations of Borazjani et al. [4].

  1. 1.

    During the acceleration phase the flow is dominated by large-scale, coherent instabilities and organized unsteadiness associated with the roll-up of the valve housing shear layer into the aortic sinus and the formation of two shear layers from the valve leaflets. For approximately the first half of the accelerating flow phase the flow exhibits very little variability from one cycle to another and the ensemble-averaged vorticity fields are nearly identical to instantaneous realizations.

  2. 2.

    The onset of significant cycle-to-cycle variations in the vorticity field is triggered by the breakdown of the leaflet shear-layers and the emergence of von Karman like vortex shedding, which occurs at approximately the middle of the accelerating phase.

  3. 3.

    At peak systole the shear layers are still reasonably well defined but their coherency is rapidly diminished as the sinus and leaflet shear layers undergo an explosive transition to a small scale turbulent state downstream of the valve. The deceleration and closing phase flow fields show little evidence of coherent flow patterns with the flow almost entirely dominated by multiple small scale eddies and complex vortical interactions.

  4. 4.

    Unlike previously believed, the chaotic, turbulent-like state that emerges in the decelerating phase never laminarizes again. Instead the chaotic state decays slowly and its remnants are washed out at the start of the accelerating flow phase shortly after the valve opens again.

  5. 5.

    The measured leaflet kinematics revealed significant variation from cycle-to-cycle, especially during the valve closing phase, and clearly showed that the leaflet motion never reaches a periodic state. Furthermore, the measured leaflet kinematics show that the valve leaflets neither open nor close symmetrically during the cardiac cycle, the asymmetry being more pronounced during the closing phase [10]. The simulations captured the observed asymmetry (Fig. 5) as well as other subtle kinematical features, such as the leaflet rebound (see the insets in Fig. 5), in spite of the fact that in the simulations the leaflets are geometrically symmetric and exposed to a symmetric incoming flow. The asymmetry in the experiments can be due to different sources, such as imperfection in the geometry of each leaflet, apparatus-specific noise and asymmetric disturbances at the inflow, etc. However, the simulations showed, for the first time, that the observed asymmetries in the computed leaflet kinematics are related to asymmetries in the flow due to natural flow instabilities [4].

Fig. 5
figure 5

Comparison of the upper and lower leaflet kinematics during the opening (top) and closing (bottom) phases. In each figure the inset shows the asymmetric rebound of the leaflets. During the opening phase the calculations are carried out with SC-FSI while both LC and SC-FSI algorithms are stable during closing. Taken from [4]

Full size image

Ge et al. [21] analyzed the results of the aforementioned experiments and simulations to investigate the mechanical environment experienced by blood elements under physiologic conditions. Until recently the Reynolds shear stress was widely considered as as the key metric that needed to be minimized to optimize different MHV designs [64, 21], mainly because of the different studies pointing at the damage caused by high shear stress—see Leverett et al. [37] for a collection of shear stress threshhold data that cause red blood cell damage. The Reynolds shear stress has been frequently used as equivalent to the viscous shear stress and Leverett et al. [37] lumped both stresses under the common name "shear stress". Ge et al. [21] showed that the so-called Reynolds shear stresses neither directly contribute to the mechanical load on blood cells nor is a proper measurement of the mechanical load experienced by blood cells and the instantaneous viscous stress provides the proper measurement of the mechanical forces in BMHV flows. They also showed that the overall levels of the viscous stresses are apparently too low (<15 N/m2 during the cardiac cycle) to induce damage to red blood cells but could potentially damage platelets (≈10 N/m2 [30, 37]). Their analysis, however, was restricted to the flow downstream of the valve leaflets and thus did not address other areas within the BMHV where potentially hemodynamically hazardous levels of viscous stresses could still occur (such as in the hinge gaps and leakage jets). They also suggested the concept of local maximum shear [29, 55] to establish a coordinate-independent metric for viscous stress tensor and underscored the significance of three-dimensionality of the flow for accurate quantification of viscous stresses.

Finally, Borazjani [3] and Borazjani et al. [5] applied the CURVIB-FSI method to carry out the first high-resolution 3D FSI simulations of a BMHV implanted in an anatomically realistic aorta obtained from MRI and compared the results with the straight aorta case (see Fig. 6). This study showed the significant effect of the aorta geometry on leaflet kinematics and the mechanical environment experienced by blood cells and underscored the need for patient-specific FSI simulations [3, 5].

Fig. 6
figure 6

Simulations of a BMHV implanted in an anatomic aorta [3, 5]. Left instantaneous out-of-plane vorticity contours on the midplane of the valve. Right 3D instantaneous vortical structures visualized by iso-surfaces of q-criterion

Full size image

4 Future outlook

The review of recent work presented in this paper underscores the major progress made in the last few years in our ability to simulate numerically BMHV flows at physiologic conditions and at resolution sufficiently high to start probing the links between valve fluid mechanics and thromboembolic complications. A major computational challenge that has yet to be tackled in this regard is the development of computational models that can elucidate the hemodynamics in microscopic regions of BMHV designs, such as the valve hinges and the leakage jet during closure, which could induce hemodynamic stresses large enough to damage blood cells. Such models should be inherently multi-scale, due to the large disparity in the macro- (aorta diameter ∼cm) and micro-scales (typical size of the gaps in the hinge region ∼102 μm), and also account for the two-phase, non-Newtonian nature of blood.

Another limitation of all existing computational models is that they have thus far treated the aorta as a rigid-wall vessel. This assumption is obviously incorrect but has been adopted so far for computational expedience since the main emphasis of previous work was on simulating and understanding the hemodynamics induced by the moving leaflets alone. The compliance of the aortic wall, however, could play an important role in the BMHV hemodyanmics and needs to be taken into account by developing complete FSI models that resolve both the valve motion and the deformation of the aorta in a coupled manner. Finally, the computational tools need to be coupled with state-of-the-art medical imaging modalities to develop a patient-specific computational framework that will allow surgeons to optimize the implantation of mechanical valves in a virtual surgery environment.

Even though these challenges are significant or even daunting, the progress we have made so far coupled with the rapidly increasing power of modern massively parallel computational platforms and advances in medical imaging allow us to be more than optimistic. The computational advances needed to meet these challenges are well within the reach of our present-day capabilities and will define the future research agenda in the area of computational hemodynamics for mechanical heart valves.