Theor. Comput. Fluid Dyn. (2014) 28:589–604 DOI 10.1007/s00162-014-0335-4 O R I G I NA L A RT I C L E Vijay Vedula · Stefania Fortini · Jung-Hee Seo · Giorgio Querzoli · Rajat Mittal Computational modeling and validation of intraventricular flow in a simple model of the left ventricle Received: 10 June 2014 / Accepted: 27 October 2014 / Published online: 9 November 2014 © Springer-Verlag Berlin Heidelberg 2014 Abstract Simulations of flow inside a laboratory model of the left ventricle are validated against experiments. The simulations employ an immersed boundary-based method for flow modeling, and the computational model of the expanding–contracting ventricle is constructed via image-segmentation. A quantitative comparison of the phase-averaged velocity and vorticity fields between the simulation and the experiment shows a reasonable agreement, given the inherent uncertainties in the modeling procedure. Simulations also exhibit a good agreement in terms of time-varying net circulation, as well as clinically important metrics such as flow-wave propagation velocity and its ratio with peak early-wave flow velocity. The detailed and critical assessment of this comparison is used to identify and discuss the key challenges that are faced in such a validation study. Keywords Hemodynamics · Left ventricle · Immersed boundary method · Validation · Color M-mode · Vortex propagation List of symbols V Tc t∗ Q ReLV W oLV FE EF Change in ventricular volume (ml) Cardiac cycle time (s) Non-dimensional time Flow-rate prescribed to the simulation Reynolds number Womersley number Vortex formation number (E-wave) Ejection fraction Communicated by Jeff D. Eldredge. V. Vedula · J.-H. Seo · R. Mittal (B) Department of Mechanical Engineering, Johns Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218, USA E-mail: mittal@jhu.edu S. Fortini Dipartimento Ingegneria Civile, Edile e Ambientale, Sapienza Università di Roma, Via Eudossiana 18, 00184 Rome, Italy G. Querzoli Dipartimento di Ingegneria Civile Ambientale e Architettura, Università Degli Studi di Cagliari, Via Marengo 3, 09123 Cagliari, Italy 590 DLV Um Am Dm ρ µ p u ω Γ CMM VE Vp ωp V. Vedula et al. Diameter of the base of the ventricle Peak area-averaged velocity through the mitral annulus Area of mitral orifice Diameter of mitral orifice Density of working fluid Dynamic viscosity of working fluid Pressure field Velocity field Vorticity, defined as the curl of the velocity field Circulation, defined as the area integral of the vorticity field Color M-mode Peak E-wave velocity measured at mitral annulus Mitral flow-wave propagation velocity Vortex propagation velocity 1 Introduction Continued advancements in imaging technologies augmented with ever-increasing computing power are enabling unprecedented analysis of cardiovascular flows [1–9]. Computational modeling of cardiac flows is also emerging as one of the viable noninvasive methods for analyzing cardiac function in health and disease [10–16]. The quantitative information available from high-fidelity numerical simulations can provide new ways of diagnosing cardiac-related diseases and help optimize surgical planning [16,17]. This provides motivation for developing a computational framework that would eventually facilitate rapid (and even real-time) patient-specific modeling of the total heart function. The earliest studies based on computational modeling of cardiac flows can be traced to Peskin et al. [18– 20] and McQueen et al. [21] which have employed the immersed boundary method. Several recent numerical studies have assumed canonical (such as semi-prolate-spheroid) left ventricles (LVs) in addition to other approximations on the expansion and contraction of the LV [4,5,10,11,13,16,22,23]. The use of simple ventricular models dates back to the seminal experiments of Bellhouse et al. [24,25], which demonstrated the utility of such models; such models have also been used in a number of recent experimental studies [12,23,26–30]. However, despite the fact that the computational modeling of intraventricular flows is now relatively well established, no study has, to our knowledge, focused on a comprehensive quantitative validation of such models. Saber et al. [1] and Schenkel et al. [2] have made qualitative comparison of their computed results with velocity fields reconstructed from two-dimensional MRI data sets for healthy hearts. Qualitative comparison of computed velocity and vorticity fields with experimental data has also been made by Doenst et al. [26] and Domenichini et al. [23]. The former has used a patient-specific model generated from cine MR to validate their model, while the latter performed validation using a simplified laboratory LV model. Recently, Krittian et al. [12] compared intraventricular flow patterns and ventricular wall deformation between in-vivo MRI data, laboratory experiment and a CFD-FSI-based simulation based on the Karlsruhe heart model [12]. However, they noted experimental limitations in obtaining realistic ejection fraction and Reynolds numbers; in addition, the comparison was mostly qualitative in nature [12]. With an eventual goal of developing a CFD-based method to analyze patient-specific heart function, our first step is to quantitatively validate the hemodynamic flow solver using a simple laboratory model of the left ventricle. Our numerical solver is based on a sharp-interface immersed boundary method [31–33], which has been validated for canonical flows [32] as well as complex flows such as flow past fish fins [34] and flapping insect wings [35]. The laboratory model of the left ventricle employed here has already been used earlier for the investigation of intraventricular flows [27–30]. We make quantitative comparison of the computed phaseaveraged velocity and vorticity profiles for this LV model and use this to assess the limitations and challenges associated with conducting such simulations. We then use our validated simulations to conduct an analysis of computed intraventricular flow; the analysis is based on velocity and vorticity distributions as well as the flow-wave propagation velocity (Vp ) obtained using “virtual” color M-mode (CMM) echocardiography. We relate Vp with respect to the propagation of the vortex ring (ωp ) which was earlier categorized in Vierendeels et al. [36]. Computational modeling and validation of intraventricular flow (a) (c) 591 (b) (d) Fig. 1 Experimental setup and configuration: a Schematic of the setup used in the experiment. b Simplified laboratory model of the left ventricle (LV) made of silicone rubber used in the experiment. c Time variation of the change in volume of the ventricle (V (ml)) with the indicated cardiac phases. d Top view of the reference planes along the longitudinal axis of the model LV 2 Methods 2.1 Experimental setup and parameters The experimental measurements have been performed in the pulse duplicator shown in Fig. 1a. This setup has been employed in previous studies [27–30], and therefore, we provide only a brief description here. The ventricle (labeled V in Fig. 1a) is flexible and made of silicone rubber that allows easy optical access. The ventricle is secured on a 56-mm-diameter circular plate (Fig. 1b) and connected to a constant head reservoir via two plexiglass tubes which serve as the mitral and the aortic conduits. Two one-way hydraulic valves are mounted inside these conduits in order to simulate the function of the native heart valves—the valve in the mitral conduit only allows flow into the ventricle [28,29] and the valve in the aortic conduit opens outward to allow only outflow. The entire ventricular assembly is housed inside a transparent plexiglass rectangular tank (labeled T in Fig. 1a). Changes in the volume of the ventricle are driven by moving the piston (labeled P in Fig. 1a), which is driven by a linear motor (labeled M in Fig. 1a). The linear motor is controlled by a personal computer by means of a speed-feedback servo-control, and this allows for the control of the time-rate-of-change of the ventricular volume. The procedure for 3D reconstruction of the ventricular flow used in the experiments is explained in detail in Fortini et al. [28]. The working fluid (distilled water) inside the ventricle is seeded with neutrally buoyant particles with a mean diameter of 30 µm. A high-speed camera (Mikrotron3620, F = 250 frames/s, 1,280 × 1,024 pixel resolution) is triggered by the motor to capture images of the seeded flow at chosen time-instants during the cycle. The acquired images are then analyzed to estimate the velocity fields over a regular grid by means of a feature tracking algorithm [29]. A series of 2D measurements over parallel planes are extracted and used to reconstruct the 3D flow. The temporal and spatial resolutions of this acquisition are 4 ms and 0.12 mm, 592 V. Vedula et al. respectively, and these are high enough to identify the vortical structures generated in the left ventricle and to follow their evolution during the cardiac cycle. Ventricular filling (diastole) is a biphasic phenomenon characterized by early filling (the “E-wave”) that corresponds to ventricular relaxation and a second peak (the “A-wave”) that is driven by the contraction of the left atrium. The interim phase, when the ventricle volume remains nearly constant, is referred to as “diastasis”. Filling is then followed by ejection (systole). The experimental setup attempts to reproduce these key features of left-ventricular diastole and systole. In Fig. 1c, the change in ventricular volume (V (ml)) is plotted as a function of normalized time, t ∗ = t/Tc where Tc is the period for one cardiac cycle. The key phases in the cardiac cycle relevant to the present investigation are also shown in Fig. 1c, and the sections where the velocity measurements were compared against the simulation data are indicated in Fig. 1d. The experiments are continued for 50 cardiac cycles, and the data are phase-averaged over these cycles. 2.2 Immersed boundary method (IBM) for flow simulation The flow solver employs a sharp-interface immersed boundary method based on the multi-dimensional ghostcell methodology (GCM) described in Mittal et al. [32,33], well suited for flows associated with complex moving and deforming bodies. Immersed boundaries are represented by surface meshes with triangular elements, and these boundaries are immersed in a Cartesian grid that covers a cuboidal domain. Further mathematical treatment of GCM and the implementation of the boundary conditions can be found in Mittal et al. [32,33]. As it is widely accepted that blood exhibits Newtonian behavior in large diameter vessels, the dynamics of blood flow in the left ventricle are modeled using the Navier–Stokes equations for unsteady, incompressible flows, ∇ · ū = 0 (1)   ∂ ū ρ + ū · ∇ ū = −∇ p + µ∇ 2 ū (2) ∂t where u is the velocity vector, p is the pressure, and ρ and µ are the density and the viscosity of the blood, respectively. The equations are discretized using a second-order, cell-centered, collocated grid arrangement of the primitive variables, u and p. A central-difference scheme is used to discretize the viscous terms, and a linear combination of the second-order central and upwind schemes is employed for the convection terms [37]. Both terms are treated implicitly using a second-order Crank–Nicolson method to eliminate the viscous stability constraint. A variant of the second-order fractional-step method proposed by Zang et al. [38] is used to solve the above governing Eqs. (1, 2) on a non-body-conforming Cartesian grid. In this version of the fractional-step method, the face-normal velocities are computed and separately updated from the cell-centered velocities. As shown by Zang et al. [38], this separate update of the face velocity results in discrete divergence of the velocity field to machine zero, even for a non-staggered primitive variable formulation. The solver has been tested extensively for fluid dynamics of complex three-dimensional stationary, moving and deforming bodies [4,5,16,32,34,35]. 2.3 Computational model of the left ventricle The present computational LV model is generated using biplanar images from the experiment along two orthogonal directions, referred to as A−A′ and B−B ′ in Fig. 1d. These images are enhanced by filtering followed by thresholding, and subsequently, the edges of the ventricle in each plane are detected from the zero-crossings of the second derivative of the image intensities [39] (see Fig. 2a). Once the edges were detected along A−A′ and B−B ′ planes, three-dimensional reconstruction is accomplished by fitting a circle at each axial plane. The center of the circle along each of these transverse planes is located at the intersection of A−A′ and B−B ′ planes while the radius is determined by an average of the distances from the center to the detected edges. Figure 2b shows the reconstructed LV at end-diastole (EDV) and end-systole (ESV). The segmentation approach described above is carried out for 1,500 frames during one cardiac cycle, resulting in a 4D (3D + time) reconstruction of the ventricular lumen. The surface of the LV model is represented by an unstructured grid with 124,340 triangular elements, and this time-evolving surface is “immersed” into the Cartesian grid as shown in Fig. 2c. The surface velocity is computed by differentiating the element coordinates in time and then prescribed to the solver as a velocity boundary condition on the lumen boundary. Computational modeling and validation of intraventricular flow 593 (a) (b) (c) Fig. 2 Details of the computational model employed in the current study. a Images of model LV acquired during experiment along A−A′ and B−B ′ planes are enhanced (Gaussian filter) and thresholded for ventricular edge detection. b Reconstructed LV from the detected edges of A−A′ and B−B ′ planes and approximating a circular cross-section axially. Nomenclature has been indicated for end-diastolic (EDV) and end-systolic (ESV) states. c Triangulated ventricle (124,340 elements) immersed into the background Cartesian grid The above procedure is subject to two main sources of uncertainty and error. The first source of error comes from the thresholding of the lumen boundary. As it is clear from Fig. 2a, the LV boundary appears as a pixelated region with a finite thickness due to imaging artifacts, and in principle, the position of the actual LV lumen could lie anywhere in the thin white pixelated band. Analysis of these images shows the thickness of this pixelated band is up to 2.0 mm which is about 8 % of the average LV radius. The second source of error is associated with the assumption of a circular cross-section at each axial location of the LV lumen, and we note a maximum deviation of the mean radius from the outer edges of either A−A′ or B−B ′ to be about 3 % of the mean LV radius. A major consequence of the above errors is that they introduce spurious temporal gradients in the endocardial velocity and consequently, in the computed flow-rate in and out of the LV. One way to alleviate this is to apply “smoothing” to the motion of the LV boundary; however, this introduces a number of ad hoc parameters and also does not ensure that the resulting flow-rate will match the flow-rate in and out of the LV as measured at the mitral inflow and aortic outlet of the LV. The matching of this velocity, especially the inflow at the mitral annulus, is crucial in order to be able to conduct a successful validation of the simulations. A decision was therefore made to directly match the segmented volume with the integral of a “target” flow-rate at the mitral annulus and adjust the LV motion to be consistent with this flow-rate. The “target” volume flow-rate at the mitral annulus was initially chosen to be the rate-of-change of volume prescribed to the piston motor of the experiment (Q in (motor), dotted line in Fig. 3a). However, simulations 594 V. Vedula et al. Fig. 3 Comparison between: a input flow-rate (Q in (motor)) and measured flow-rate (Q). (Q) is the uncertainty in the measured flow-rate prescribed to the simulation (Q) due to axisymmetric approximation. For clarity of comparison, the Q curve is shifted by 20 ml/s on the plot along the positive ordinate, b volume of LV  computed directly from the segmentation procedure (V0 ) and the integral of “target” or measured flow-rate at the mitral inlet ( Qdt) conducted using Q in (motor) as the target produced vortex propagation speeds that were noticeably lower than those in the experiment. The flow-rate was then estimated directly from the PIV measurements at the mitral annulus by using an average of the measurements along the A−A′ and C−C ′ planes assuming an axisymmetric flow; this flow-rate profile is also shown in Fig. 3a as a solid line (Q), and it is easily verified that the actual flow-rate profile at the mitral inlet was significantly different from Q in (motor) in both amplitude and phase. For systole, the flow-rate profile was estimated based on outflow measured at the A−A′ plane. An axisymmetric flow was then assumed for the outflow as well, and the contraction rate is multiplied by a factor (equal to the ratio of the total inflow to the integral of predicted systolic flow-rate) to match the net inflow or stroke volume. Also shown in Fig 3a is the error (Q) associated with the assumption of axisymmetry at the mitral inflow. This is computed as the difference between the flow-rates estimated from the profiles along the A−A′ and C−C ′ planes, and we note errors up to 11 % of the peak diastolic flow-rate. This error introduces inherent uncertainty into our computational model, and this needs to be kept in mind while making comparisons with the experimental data. Figure 3b shows the time variation of the LV volume, V0 , obtained directly  from the segmentation compared to the integral of the measured flow-rate at mitral inlet (denoted as Qdt on the plot). It is clear that there is a mismatch between the two, and in fact, analysis indicates this mismatch is up to 5.3 % of the mean LV volume. Due to the incompressibility constraint, the two volumes have to match each other at each time-step and we accomplish this by adjusting the LV boundary at each time-step. To elaborate, we use an iterative procedure where the segmented volume of the ventricle (V0 in Fig. 3b) is modified at  selected key frames to match this integral ( Qdt). In this iterative process, the radius of the circular crosssection at each axial location of the segmented LV is changed until the LV volume matches the target volume within a prescribed tolerance. We chose 125 frames over one cardiac cycle to apply this matching procedure, and this enables us to obtain a reasonably smooth variation in the resulting flow-rate curve. The iterative method converges rapidly for the current problem, and a pseudo-code for this procedure is provided in the “Appendix”. The presence of valves is modeled in the simulations by opening and closing the far ends of the mitral and the aortic annuli. During diastole, the mitral orifice is kept open by imposing a zero gradient on the velocity and a Dirichlet condition for pressure while the aorta is closed using a no-slip, no-penetration boundary condition. Conversely, during systole, the mitral annulus is closed while the aorta is opened by interchanging the aforementioned boundary conditions. For the sake of simplicity, the Dirichlet value for pressure at the inflow boundaries is set to 0; this will have no effect on the pressure gradients in the flow. It is to be noted that the phase-average in the simulations is computed over five cardiac cycles, while the experiment estimates the phase-average over 50 cycles. It is also noted that the phase-average in the simulations does not show any significant cycle-to-cycle variation beyond the first two cycles. The shorter duration of phaseaveraging in the simulations is due partly to the limitations of computational time and resources and also to the fact that there is greater inherent uncertainty and cycle-to-cycle variation in the experiment compared to the simulation which necessitates additional filtering. Computational modeling and validation of intraventricular flow 595 Table 1 Key parameters used for the simplified model left ventricle (LV) ReLV W oLV EF (%) E/A FE 8120 22.8 40.0 3.0 3.94 2.4 Non-dimensional parameters Table 1 indicates hemodynamically significant parameters for the current model. Following Fortini et al. [28], the pulsatile flow in the left ventricle is characterized in terms of two non-dimensional parameters: the Um Reynolds number (ReLV = ρ DLV ). In these expressions, ) and the Womersley number (W oLV = √TDLV µ c µ/ρ DLV is the diameter of the base of the left ventricle, Um is the peak, area-averaged velocity through the mitral annulus (=ratio of peak diastolic flow-rate to the annular area), and ρ and µ are the density and dynamic viscosity of the working fluid, respectively. The other key non-dimensional parameters for the ventricular flow are the ejection fraction [EF-defined as the ratio of stroke volume (EDV–ESV) to the end-diastolic volume (EDV)], the ratio of the peak velocities of early filling to atrial filling (E/A ratio) and the vortex formation d number [40] (FE = AmVD ) where Vd represents the flow volume discharged through the mitral orifice during m E-wave, and Am and Dm represent the area and diameter of the mitral orifice, respectively. The values of all these parameters for the current model are also reported in Table 1. It is noted that ReLV , W oLV and FE are within a reasonable physiological range for the current model and is therefore suitable for validation purposes. 2.5 Metrics for quantitative comparison of computational and experimental data We use the phase-averaged velocity and vorticity (ω, defined as the curl of the velocity) fields to compare computational and experimental data sets. Velocity profiles are compared along A−A′ plane at various phases of the cardiac cycle and at different cross-sections along this plane. The error in the velocity profiles is quantified using an index, E L1 = 100 Um (ū c f d − ū exp 1 ), which is defined as the L 1 norm of the difference in the velocity between the simulation and the experiment expressed as a percentage of Um . Another metric used for comparison between simulations and experiment is the time variation of the circulation inside LV cavity. The net non-dimensional circulation (Γ ) is defined as the area integral of the out-of-plane vorticity over the A−A′ plane and is equal to Γ = TSc S ωdS, where S is the integrated area along A−A′ plane which comprises only the ventricular cavity excluding sections of mitral and aortic channels. It is noted that the positive values of Γ are associated with counterclockwise circulation while its negative counterparts are associated with clockwise circulation. The color M-mode (CMM) of Doppler assessment of flow propagation velocity (Vp ) is routinely used in clinical cardiology to assess ventricular diastolic function [41]. As such, it is useful to assess the accuracy of the simulated flow via this modality. In clinical practice, the color M-mode is acquired by aligning the ultrasound probe in the direction from the ventricular apex to the mitral annulus, and the principle of Doppler shift is used to estimate the flow velocity directed toward the apex along the entire long axis of the ventricle [41,42]. In the present study, we follow a procedure that has been used in previously reported studies [5,16,36] to create a “virtual” CMM using simulation data. In this method, the longitudinal component of velocity along the mitral axis from base to apex is extracted directly from the simulation data. We note that this method of computing CMM based on computational modeling has not been validated, although a qualitative comparison with a simple Doppler ultrasound model has been provided in Seo et al. [5]. The velocity information obtained from this virtual CMM is then plotted with time as abscissa and the distance from the base to apex as ordinate and appropriately color-coded. Vp is then determined by the slope of the first aliasing velocity (transition from blue to red) where the aliasing limit is set to 50–70 % of the maximal velocity spread distally 4 cm into the LV cavity [42]. The ratio of Vp to peak E-wave velocity (VE ) was shown to be a potential marker of diastolic dysfunction [36,42,43], and we compare these metrics between the experiment and the simulation. In addition to this, we also probe the location of the core of the vortex ring (identified using λci -criterion [44,45]) superposed on the CMM. The slope of a linear fit to this curve represents the vortex propagation velocity, ωp which is also computed for comparison in the present analysis. 596 V. Vedula et al. Fig. 4 Results due to grid convergence (very coarse—64 × 64 × 128; coarse—96 × 96 × 192, medium—128 × 128 × 256, fine—256 × 256 × 512). a Comparison of velocity components along the centerline of the mitral orifice during the deceleration phase of the E-wave, t ∗ = 0.22. b Time variation of the total kinetic energy in the LV cavity (excluding mitral and aortic channels) for all the chosen grids 2.6 Grid refinement study In order to ensure adequate spatial accuracy, numerical simulations have been performed on four different computational grids: (a) very coarse (64 ×64 ×128), (b) coarse (96 ×96×192), (c) medium (128×128×256) and (d) fine (256 × 256 × 512). The size of the time-step (t) in the simulations is varied appropriately (0.0004 Tc − 0.0001 Tc ) for each chosen grid so that the maximum Courant–Friedrichs–Lewy (CFL) number is about 0.3. Figure 4a shows a comparison of the profiles of the three velocity components along the mitral axis during the deceleration phase of the E-wave (t ∗ = 0.22), and Fig. 4b shows a comparison of the total kinetic energy of the fluid in the LV cavity (excluding mitral and aortic channels) as a function of time for the chosen grids. A reasonable convergence was achieved on the fine grid for this highly transitional flow. However, due to limitations of computational time and resources, we have computed the flow on the fine grid only for one cycle. Hence, we use the fine grid only for understanding the flow and the vortical features while the medium grid is chosen for phase-averaging over multiple cycles and subsequent comparison with the experiment. It is to be noted that a simulation of one cardiac cycle on the medium grid takes about 100 h with 256 cores on a Cray XT5 system; this cluster has a total of 9,408 compute nodes, each with two six-core 2.6 GHz AMD Opteron processors. 3 Results 3.1 Vortical features We begin with an analysis of the three-dimensional (3D) vortical features of the intraventricular flow identified in the simulation data by the λci -criterion, which is defined as the imaginary part of the complex eigenvalues of the velocity gradient tensor [44,45]. The isosurfaces of this quantity shown in Fig. 5 for λci = 10 s−1 are color-coded by the longitudinal component of the velocity, and the corresponding cardiac phase is also indicated in each frame as a highlighted circle on the flow-rate curve. It is noted that during the first cycle, filling takes place into a static fluid chamber, i.e., the fluid inside the ventricular model is initially stationary without any remnants. It has been well established that under normal resting conditions, a vortex ring is formed inside the human left ventricle during filling [3,4,7,9,40]. In-line with this, we notice a strong vortex ring formed at the mitral orifice at the peak of E-wave in Fig. 5a. The formation number of this E-wave vortex ring is FE = 3.94, and we observe that this strong vortex ring pinches off from the base of LV and convects toward the apex. The trailing shear layers behind the vortex ring are nearly nonexistent, and this is consistent with the fact that FE is very close to the “optimal” value of 4.0 [40]. Early in the vortex ring formation and propagation (Fig. 5b), the ring interacts with the lateral wall of the ventricle causing it to tilt. Continued interaction of the ring with the lateral wall leads to disruption of the vortex ring (Fig. 5c), and a further tilting allows the left side of the vortex ring to fill up the central cavity (Fig. 5c) of the ventricle. Computational modeling and validation of intraventricular flow 597 Fig. 5 Time evolution of the three-dimensional vortical structures of the intraventricular flow visualized using λci -criterion [44,45] as the isosurfaces of imaginary part of complex eigenvalues of the velocity gradient tensor colored by the longitudinal component of velocity for value, λci = 10 s−1 at the indicated times Subsequently, the vortex impinges on the right wall and disintegrates completely leading to formation of smaller eddies that fill up the ventricular cavity (Fig. 5d). From this point onward, the flow exhibits traits of transition, with a cascade of eddies of various length scales as evident in Fig. 5d–g. However, due to the overall circulation being established in the clockwise direction, the fluid inertia pushes these dissipating eddies toward the aorta by the end of diastasis (Fig. 5e, f). Although the present model is very simplistic, flow features described here have also been identified in 4D patient-specific models [6,7,9]. Following this, the atrial filling due to A-wave produces a weak vortex near the mitral orifice which immediately interacts with the remnants of the E-wave vortex (Fig. 5g). The weak A-wave vortex does not seem to play any role in re-energizing the dissipating eddies. By the end of the filling, the stronger eddies which appeared during diastole seem to have been mostly dissipated. Systole ensues during which the fluid gets smoothly ejected through the outflow tract (Fig. 5h) that acts as a “sink.” The vortex structures that remain at the end of the diastole get stretched and aligned with this outwardly directed flow. 3.2 Validation with experimental data Phase-averaged velocity profiles along the A−A′ plane are compared between the simulation and the experiment at various phases of the cardiac cycle in Fig. 6. The vertical component of velocity is plotted along three horizontal cross-sections (H1, H2, H3), while the lateral component of velocity is plotted along three vertical cross-sections (V1, V2, V3) of A−A′ plane. The experimental data are shown as symbols, while the numerical data are plotted using solid lines. A visual comparison of the two sets of profiles indicates that the simulation reproduces all the key features of the experimental velocity profiles reasonably well. Some general trends observed are that the primary differences between the experimental and the simulated flow profiles are in the 598 V. Vedula et al. Fig. 6 Comparison of profiles of phase-averaged longitudinal (red) and lateral (blue) components of velocity between the computational data and that of the experiment (symbols) at various cross-sections of the LV along the A−A′ plane (see Fig. 1d). Reference vectors in frames A and D correspond to 10 cm/s along each direction (color figure online) Table 2 Quantification of error E L1 between simulation and experimental velocity profiles (Fig. 6) expressed as percentages of Um t∗ H2 H3 V2 V2 V3 0.172 11.5 2.0 0.22 7.4 8.4 0.28 3.7 6.9 0.36 1.7 5.4 0.56 5.4 7.6 0.84 6.3 1.9 Errors greater than 10 % have been italicized H1 0.7 1.4 6.7 7.4 4.0 1.5 1.6 4.0 9.7 9.4 5.6 6.5 6.2 10.5 15.0 13.7 5.4 8.2 3.8 6.7 6.2 8.6 3.9 5.5 vertical component of velocity in early- and mid-diastole and that these differences are mostly concentrated in the region near the lateral wall of the ventricle. It is also noted that the velocity profiles from the simulation are relatively smooth at all phases except for t ∗ = 0.56, where the simulation data are significantly “rougher” than the experimental ones. Quantitative comparisons of the velocity profiles are performed in Table 2 that shows EL1 error between the simulated and the experimental velocity profiles for each of the designated sections on the A−A′ plane and at the corresponding cardiac phase. We notice that except for two profiles (H 1 at t ∗ = 0.172 and V 2 at t ∗ = 0.22, 0.28, 0.36), the difference between these data sets is <10 % elsewhere. The highest E L1 error is 15 %, and this occurs for V 2 at t ∗ = 0.36. In Figs. 7 and 8, we compare contours of phase-averaged out-of-plane vorticity on the A−A′ and B−B ′ planes at the indicated non-dimensional times, respectively. The plots in the top row correspond to the simulation, while the bottom row shows the corresponding plots for the experiment. A visual comparison indicates a reasonable agreement in the overall vorticity pattern between the simulation and the experiment at all these phases of the cardiac cycle on both these planes. In Fig. 7, the comparison of the clockwise (blue contours) Computational modeling and validation of intraventricular flow 599 Fig. 7 Comparison of the computed phase-averaged out-of-plane vorticity (top row) with the corresponding phase-averaged data from the experiment (bottom frame) along the A−A′ plane (see Fig. 1d) at various non-dimensional times Fig. 8 Comparison of the computed phase-averaged out-of-plane vorticity (top row) with the corresponding phase-averaged data from the experiment (bottom frame) along the B−B ′ plane (see Fig. 1d) at various non-dimensional times. The dark circle indicates the location of the core of the vortex used for computing ωp rotating vorticity is particularly good; the shape as well as the location of the clockwise vorticity matches well at these phases of diastole. On the other hand, the comparison between the counterclockwise (red contours) rotating vorticity may be considered good only up to t ∗ = 0.22. At t ∗ = 0.28, it is noted that the shape of this 600 V. Vedula et al. Fig. 9 Comparison of non-dimensional circulation, Γ along A−A′ plane between the simulation and the experiment over the entire cardiac cycle (color figure online) Fig. 10 Comparison of “virtual” CMM between a simulation and b experiment for the present LV model. Vp is measured by the slope of the thick yellow colored line, while ωp is measured by the slope of the linear fit to the locus of the position of vortex ring indicated by symbols. Comparison between the numerical values is reported in Table 3 vortex is not predicted well and at t ∗ = 0.36, the simulation does not accurately capture the location and size of this vortex. Beyond this phase in the cardiac cycle, the magnitude of the peak vorticity drops to exceedingly low values, thereby making a comparison of this quantity very difficult. The comparison of both clockwise and counterclockwise vorticity agrees reasonably well along B−B ′ plane in terms of magnitude, phase and the location of the vortex as shown in Fig. 8. However, a careful observation reveals a slight tilt of the vortex ring in the experiment (t ∗ = 0.22 and 0.28), while the simulation predicts a nearly symmetric vortex structure. In Fig. 9, we compare the time variation of Γ computed over the A−A′ plane between the simulation and the experiment. The experimental data indicates a low-level oscillation in this quantity during early diastole (t ∗ < 0.15) before it settles into a slowly varying negative value during the rest of the diastolic phase indicative of a net clockwise circulation during this period. This is consistent with Fig. 7, which shows the dominance of the clockwise portion of the vortex ring in diastole. The peak of the E-wave is associated with an increased negative circulation, whereas the peak of the A-wave corresponds to an increase toward positive circulation. The beginning of the systole shows a rapid reversal in the direction of circulation followed by decay at endsystole. The simulations predict all of these trends reasonably well. The largest quantitative mismatch between the simulation and experiment for this quantity is in late diastole, where the simulations tend to over-predict the magnitude of circulation. This is consistent with our previous observation that the counterclockwise vortex is weaker in the simulation than in the experiment. The quantitative match between the simulation and experiment during systole is quite good. The final comparison between the simulation and the experiment is performed using “virtual” CMM and the associated metrics, Vp and ωp . Figure 10a, b shows color-maps of axial velocity as a function of height and time for the computational and experimental data sets, respectively. A comparison of the two plots shows Computational modeling and validation of intraventricular flow 601 Table 3 Comparison of metrics computed using “virtual” CMM shown in Fig. 10 Parameter Simulation Experiment % Diff. VE (cm/s) Vp (cm/s) VE /Vp ωp (cm/s) 16.8 5.6 3.0 3.6 15.5 5.1 3.0 3.9 8.7 8.3 0.3 8.6 that while the initial E-wave propagation is similar for the two, the wave in the simulation propagates deeper (∼9.7 cm) into the LV than for the experiment (∼9 cm). This is also evident from Fig. 6a–c from the velocity profiles for the sections, H 1−H 3. The CMM for the simulation also shows higher fluctuations and reversed flow during diastasis and the A-wave (t ∗ = 0.45−0.7) when compared to that of the experiment. The locus of the position of the vortex core (indicated in Fig. 8 as a dark circle) is also superposed on the CMM plots in Fig. 10, and it indicates that the vortex ring lags behind the flow for this LV model in both simulation and experiment. Both Vp and ωp can be estimated from the slopes of the respective lines in the Fig. 10, and these are reported in Table 3. Also reported in the table is the peak E-wave velocity (VE ) along the probe direction at the mitral orifice level. It is noted that the difference in VE between the simulation and experiment is about 8.7 % of the experimental value. This is comfortably within the uncertainty of the flow-rate as shown in Fig. 3 and discussed in Sect. 2.3. It is interesting to note that the values of Vp and ωp are also different from the experiment by about 8–9 %. Not surprisingly, the value of VE /Vp , a quantity that is considered to have diagnostic significance [36,42,43], is matched very well between the experiment and the simulation (Table 3). 4 Discussion A detailed quantitative as well as qualitative comparison of simulation results with those of the experiment has been presented. The comparison between the simulation and the experimental data is generally reasonable, with some computed quantities (and metrics) agreeing better with the experiment than others. As with any computational model, there exist sources of uncertainty and an attempt has been made to identify and even quantify these uncertainties. This further provides an opportunity to assess the degree to which modeling uncertainties and limitations affect the flow prediction. First, it is generally noted that the differences between the simulated and experimental velocity profiles (Fig. 6) are mostly below 10 % (Table 2), with only a few selected profiles exceeding this error value. It should be noted that the uncertainty in the estimation of the diastolic flow-rate (as described in Sect. 2.3) is about 11 % of the peak diastolic flow-rate. Thus, the level of discrepancy between the two velocity profiles is very much in-line with this inherent uncertainty in the modeling procedure. It is also noted that the majority of the larger differences in the two velocity profiles (Fig. 6) are localized to the region near the right (or lateral) wall of the ventricle. This mismatch is clearly correlated with the mismatch of the portion of the vortex ring closer to the lateral wall as evident from Fig. 7. It can be seen here that the right vortex decays faster in the simulation due to its interaction with the wall while it persists much longer in the experiment. This vortex–wall interaction is highly sensitive to the precise location and motion of the lateral wall relative to the inflow jet. Our assumption of a locally circular cross-section (Sect. 2.3) as well as our method to match the ventricular volume to the measured flow-rate is expected to introduce discrepancies in the exact location of the ventricular walls and would modify the vortex–boundary layer interaction in this region. Furthermore, it is noticed from Fig. 6a, b that there is a horizontal shift in the mitral jet between the experiment and simulation, with the jet in the simulation being directed closer to the right wall. Given that the jet in the simulation is centered along the axis of the mitral annulus (Fig. 6a), this implies that the difference is due to a slight jet tilt in the experiment from the centerline of the mitral annulus. A slight tilt in the mitral flow is also evident in the vorticity structure shown in Fig. 8 along B−B ′ plane. The cause for this tilt in the experiment is not clear but likely has to do with a slight asymmetry in the flow exiting the mitral inflow pipe that is not captured in the simulation. The leftward shift of the jet in the experiment provides a larger area for the growth of the near-wall portion of the vortex ring and leads to a stronger counter-clockwise vortex in that region (Fig. 7b–d). 602 V. Vedula et al. During the late stages of diastole (Fig. 6d, e), it is noted that the peaks in the computational velocity profiles are larger than those in the experiment; this has two possible explanations. First, as seen from the CMM diagrams, the jet penetrates deeper into the LV in the simulations than in the experiment. This forces the simulated jet into the narrower apical region of the LV leading to a stronger deflection of the flow toward the left wall in Fig. 6. Second, while the experimental results have been phase-averaged over 50 cycles, the simulations are averaged over only five cycles. Given the transitional nature of the jet, it is expected that averaging over many cycles will smoothen out the sharp variation in the velocity profiles. We have limited ability to test this hypothesis given the significant computational expense associated with these high-resolution simulations. It is well known that early ventricular filling is a suction phenomenon during which the velocity increases from the base to apex and this can be captured clinically using color-coded Doppler M-mode echocardiography [36,42,43]. Such a “virtual” CMM is presented in Fig. 10 for both computational and experimental data. Immediately, one can identify a strong E-wave flow propagating from the basal region into the ventricular cavity and the velocity with which it propagates is called the “flow-wave” propagation velocity (Vp ) [36,42,43,46]. This is not equal to the velocity of the blood particles which, in fact, tend to move at much higher velocities within the vortex ring. Instead, it represents the delay in ventricular filling partly due to varying degree of expansion of the ventricle along its length together with the roll-up of the shear layer in the form of a vortex ring. From Table 3, it is clear that the propagation speed of the vortex ring (ωp ) is slower than that of Vp which is not in-line with the conclusion of Vierendeels et al. [36] that the vortex ring propagates at the same speed as that of the flow. However, their model assumes axisymmetric filling into a long and slender prolate-spheroidal shaped LV model, whereas the present LV model has a low aspect ratio with a large cavity at the center. Moreover, the filling takes place asymmetrically into this cavity from the mitral annulus, thus allowing one side of the vortex to fill up the cavity while subjecting the other to wall interaction thereby creating this velocity difference. 5 Conclusions A comprehensive quantitative validation of intraventricular flow between a highly resolved simulation and an experiment has been conducted. The major sources of modeling uncertainty are identified and quantified. Among these, the uncertainty in the estimation of the instantaneous flow-rate as well as the non-axisymmetric nature of the mitral flow in the experiment was found to be the most significant. The phase-averaged velocity profiles and vorticity contours from the simulation were found to match reasonably well with the experiment given the uncertainty inherent in the modeling procedure. The trends in the time variation of net circulation in ventricular flow were also predicted well by the simulations. Flow-wave propagation velocity (Vp ) and the ratio, VE /Vp computed from “virtual” color M-mode (CMM) as well as the vortex propagation velocity, showed a very good agreement between simulation and experiment. The study demonstrates the challenge of modeling the flow even in a relatively simple experimental model of the left-ventricle and underscores the need for continued quantitative validation to clearly understand the limitations of such simulations. Acknowledgments This research is supported by the National Science Foundation through Grants IOS-1124804 and IIS1344772. Computational resources were provided through TeraGrid Grant TG-CTS100002. The experimental work was supported by the Italian MIUR, through Grant PRIN-2009J7BL32. Appendix LV volume matching A pseudo-code for matching the volume of LV computed after segmentation with that of the integral or “target” volume is presented here. Computational modeling and validation of intraventricular flow 603 References 1. Saber, N.R., Wood, N.B., Gosman, A.D., Merrifield, R.D., Yang, G.Z., Charrier, C.L., Gatehouse, P.D., Firmin, D.N.: Progress towards patient-specific computational flow modeling of the left heart via combination of magnetic resonance imaging with computational fluid dynamics. Ann. Biomed. Eng. 31, 42–52 (2003) 2. Schenkel, T., Malve, M., Reik, M., Markl, M., Jung, B., Oertel, H.: MRI-based CFD analysis of flow in a human left ventricle: methodology and application to a healthy heart. Ann. Biomed. Eng. 37(3), 503–515 (2009) 3. Töger, J., Kanski, M., Carlsson, M., Kovács, S.J., Söderlind, G., Arheden, H., Heiberg, E.: Vortex ring formation in the left ventricle of the heart: analysis by 4D flow MRI and Lagrangian coherent structures. Ann. Biomed. Eng. 40, 2652–2662 (2012) 4. Seo, J.H., Mittal, R.: Effect of diastolic flow patterns on the function of the left ventricle. Phys. Fluids 25(11), 110801:1– 110801:21 (2013) 5. Seo, J.H., Vedula, V., Abraham, T., Mittal, R.: Multiphysics computational models for cardiac flow and virtual cardiography. Int. J. Numer. Methods Biomed. Eng. 29(8), 850–869 (2013) 6. Mihalef, V., Ionasec, R.I., Sharma, P., Georgescu, B., Voigt, I., Suehling, M., Comaniciu, D.: Patient-specific modeling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac CT images. J. R. Soc. Int. Focus 1, 286– 296 (2011) 7. Le, T.B., Sotiropoulos, F.: On the three-dimensional vortical structure of early diastolic flow in a patient-specific left ventricle. Eur. J. Biomech. B/Fluids 35, 20–24 (2012) 8. Le, T.B., Sotiropoulos, F., Coffey, D., Keefe, D.: Vortex formation and instability in the left ventricle. Phys. Fluids 24(9), 091110:1–091110:2 (2012) 9. Chnafa, C., Mendez, S., Nicoud, F.: Image-based large-eddy simulation in a realistic left heart. Comput. Fluids 94, 173– 187 (2014) 10. Cheng, Y., Oertel, H., Schenkel, T.: Fluid–structure coupled CFD simulation of the left ventricular flow during diastole phase. Ann. Biomed. Eng. 33, 567–576 (2005) 11. Domenichini, F., Pedrizzetti, G., Baccani, B.: Three-dimensional filling flow into a model left ventricle. J. Fluid Mech. 539, 179–198 (2005) 604 V. Vedula et al. 12. Krittian, S., Schenkel, T., Janoske, U., Oertel, H.: Partitioned fluid–solid coupling for cardiovascular blood flow: validation study of pressure-driven fluid–domain deformation. Ann. Biomed. Eng. 38(8), 2676–2689 (2010) 13. Pedrizzetti, G., Domenichini, F.: Nature optimizes swirling flow in the human left ventricle. Phys. Rev. Lett. 95, 108101:1– 108101:4 (2005) 14. Taylor, C.A., Draney, M.T.: Experimental and computational methods in cardiovascular fluid mechanics. Ann. Rev. Fluid Mech. 36, 197–231 (2004) 15. Wong, K.K.L., Kelso, R.M., Worthley, S.G., Sanders, P., Mazumdar, J., Abbott, D.: Cardiac flow analysis applied to phase contrast magnetic resonance imaging of the heart. Ann. Biomed. Eng. 37(8), 1495–1515 (2009) 16. Zheng, X., Seo, J.H., Vedula, V., Abraham, T., Mittal, R.: Computational modeling and analysis of intracardiac flows in simple models of left ventricle. Eur. J. Biomech. B/Fluids 35, 31–39 (2012) 17. Yang, W., Feinstein, J.A., Marsden, A.L.: Constrained optimization of an idealized Y-shaped baffle for the Fontan surgery at rest and exercise. Comput. Methods Appl. Mech. Eng. 199(33–36), 2135–2149 (2010) 18. Peskin, C.S.: Flow patterns around heart valves: a numerical method. J. Comput. Phys. 10, 252–271 (1972) 19. Peskin, C.S.: Numerical analysis of blood flow in the heart. J. Comput. Phys. 25, 220–252 (1977) 20. Peskin, C.S., McQueen, D.M.: A three-dimensional computational model of blood flow in the heart: I. Immersed elastic fibers in a viscous incompressible fluid. J. Comput. Phys. 81, 372–405 (1989) 21. McQueen, D.M., Peskin, C.S.: A three-dimensional computational model of blood flow in the heart: II. Contractile fibers. J. Comput. Phys. 82, 289–297 (1989) 22. Baccani, B., Domenichini, F., Pedrizzetti, G., Tonti, G.: Fluid dynamics of the left ventricular filling in dilated cardiomyopathy. J. Biomech. 35, 665–671 (2002) 23. Domenichini, F., Querzoli, G., Cenedese, A., Pedrizzetti, G.: Combined experimental and numerical analysis of the flow structure into the left ventricle. J. Biomech. 40, 1988–1994 (2007) 24. Bellhouse, B.J.: Fluid mechanics of a mitral valve and left ventricle. Cardiovas. Res. 6, 199–210 (1972) 25. Bellhouse, B.J., Bellhouse, F.H.: Fluid mechanics of the mitral valve. Nature 224, 615 (1969) 26. Doenst, T., Spiegel, K., Reik, M., Markl, M., Hennig, J., Nitzsche, S., Beyersdorf, F., Oertel, H.: Fluid-dynamic modeling of the human left ventricle: methodology and application to surgical ventricular reconstruction. Ann. Thorac. Surg. 87, 1187– 1195 (2009) 27. Espa, S., Badas, M.G., Fortini, S., Querzoli, G., Cenedase, A.: A Lagrangian investigation of the flow inside the left ventricle. Eur. J. Biomech. B/Fluids 35, 9–19 (2012) 28. Fortini, S., Querzoli, G., Espa, S., Cenedese, A.: Three-dimensional structure of the flow inside the left ventricle of the human heart. Exp. Fluids 54(11), 1609:1–1609:9 (2013) 29. Querzoli, G., Fortini, S., Cenedese, A.: Effect of the prosthetic mitral valve on vortex dynamics and turbulence on the left ventricular flow. Phys. Fluids 22, 041901:1–041901:10 (2010) 30. Vukćević, M., Fortini, S., Querzoli, G., Espa, S., Pedrizzetti, G.: Experimental study of an asymmetric heart valve prototype. Eur. J. Biomech. B/Fluids 35, 54–60 (2012) 31. Mittal, R., Iaccarino, G.: Immersed boundary methods. Ann. Rev. Fluid Mech. 37, 239–261 (2005) 32. Mittal, R., Dong, H., Bozkurttas, M., Najjar, F.M., Vargas, A., von Loebbecke, A.: A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries. J. Comput. Phys. 227, 4825–4852 (2008) 33. Seo, J.H., Mittal, R.: A sharp-interface immersed boundary method with improved mass conservation and reduced spurious pressure oscillations. J. Comput. Phys. 230, 7347–7363 (2011) 34. Dong, H., Bozkurttas, M., Mittal, R., Madden, P., Lauder, G.V.: Computational modeling and analysis of the hydrodynamics of a highly deformable fish pectoral fin. J. Fluid Mech. 645, 345–373 (2010) 35. Zheng, L., Hedrick, T.L., Mittal, R.: A multi-fidelity modeling approach for evaluation and optimization of wing stroke aerodynamics in flapping flight. J. Fluid Mech. 721, 118–154 (2013) 36. Vierendeels, J.A., Dick, E., Verndock, P.R.: Hydrodynamics of color M-mode Doppler flow wave propagation velocity V (p): a computer study. J. Am. Soc. Echocardiogr. 15, 219–224 (2002) 37. Ghias, R., Mittal, R., Dong, H.: A sharp interface immersed boundary method for compressible viscous flows. J. Comput. Phys. 225, 528–553 (2007) 38. Zang, Y., Street, R.L., Koseff, J.R.: A non-staggered grid, fractional step method for time-dependent incompressible Navier– Stokes equations in curvilinear coordinates. J. Comput. Phys. 114, 18–33 (1994) 39. Rangayyan, R.M.: Detection of regions of interest. In: Rangayyan, R.M., Neuman, M.R. (eds.) Biomedical Image Analysis, pp. 363 CRC Press–Taylor and Francis, Boca Raton (2004) 40. Gharib, M., Rambod, E., Kheradvar, A., Sahn, D.J., Dabiri, J.O.: Optimal vortex formation as an index of cardiac health. Proc. Natl. Acad. Sci. 103(16), 6305–6308 (2006) 41. Nagueh, S.F., Appleton, C.P., Gillebert, T.C., Marino, P.N., Oh, J.K., Smiseth, O.A., Waggoner, A.D., Flachskampf, F.A., Pellikka, P.A., Evangelista, A.: Recommendations for the evaluation of left ventricular diastolic function by echocardiography. J. Am. Soc. Echocardiogr. 22(2), 107–133 (2009) 42. Garcia, M.J., Ares, M.A., Asher, C., Rodriguez, L., Vandervoort, P., Thomas, J.D.: An index of early left ventricular filling that combined with pulsed Doppler peak E velocity may estimate capillary wedge pressure. J. Am. Coll. Cardiol. 29(2), 448– 454 (1997) 43. De Boeck, B.W.L., Oh, J.K., Vandervoort, P.M., Vierendeels, J.A., van der Aa, R.P.L.M., Cramer, M.J.M.: Colour M-mode velocity propagation: a glance at intra-ventricular pressure gradients and early diastolic ventricular performance. Eur. J. Heart Fail. 7, 19–28 (2005) 44. Chong, M.S., Perry, A.E., Cantwell, B.J.: A general classification of three dimensional flow fields. Phys. Fluids 2(5), 765– 777 (1990) 45. Mittal, R.: Planar symmetry in the unsteady wake of a sphere. AIAA J. 37(3), 388–390 (1999) 46. Greenberg, N.L., Vandervoort, P.M., Firstenberg, M.S., Garcia, M.J., Thomas, J.D.: Estimation of diastolic intraventricular pressure gradients by Doppler M-mode echocardiography. Am. J. Physiol. Heart Circ. Physiol. 280, H2507–H2515 (2001)