 Research
 Open Access
 Published:
Computational fluid dynamics simulations of blood flow regularized by 3D phase contrast MRI
BioMedical Engineering OnLine volume 14, Article number: 110 (2015)
Abstract
Background
Phase contrast magnetic resonance imaging (PCMRI) is used clinically for quantitative assessment of cardiovascular flow and function, as it is capable of providing directlymeasured 3D velocity maps. Alternatively, vascular flow can be estimated from modelbased computation fluid dynamics (CFD) calculations. CFD provides arbitrarily high resolution, but its accuracy hinges on model assumptions, while velocity fields measured with PCMRI generally do not satisfy the equations of fluid dynamics, provide limited resolution, and suffer from partial volume effects. The purpose of this study is to develop a proofofconcept numerical procedure for constructing a simulated flow field that is influenced by both direct PCMRI measurements and a fluid physics model, thereby taking advantage of both the accuracy of PCMRI and the high spatial resolution of CFD. The use of the proposed approach in regularizing 3D flow fields is evaluated.
Methods
The proposed algorithm incorporates both a Newtonian fluid physics model and a linear PCMRI signal model. The model equations are solved numerically using a modified CFD algorithm. The numerical solution corresponds to the optimal solution of a generalized Tikhonov regularization, which provides a flow field that satisfies the flow physics equations, while being close enough to the measured PCMRI velocity profile. The feasibility of the proposed approach is demonstrated on data from the carotid bifurcation of one healthy volunteer, and also from a pulsatile carotid flow phantom.
Results
The proposed solver produces flow fields that are in better agreement with direct PCMRI measurements than CFD alone, and converges faster, while closely satisfying the fluid dynamics equations. For the implementation that provided the best results, the signaltoerror ratio (with respect to the PCMRI measurements) in the phantom experiment was 6.56 dB higher than that of conventional CFD; in the in vivo experiment, it was 2.15 dB higher.
Conclusions
The proposed approach allows partial or complete measurements to be incorporated into a modified CFD solver, for improving the accuracy of the resulting flow fields estimates. This can be used for reducing scan time, increasing the spatial resolution, and/or denoising the PCMRI measurements.
Background
Knowledge of blood flow patterns in the human body is a critical component in cardiovascular disease research and diagnosis. Two different approaches to 3D flow assessment are currently available to the researcher and clinician: direct, modelindependent velocity mapping using phase contrast magnetic resonance imaging (PCMRI) [1–3] or Doppler ultrasound, and modelbased computational fluid dynamics (CFD) calculations [4–16]. Among the direct methods, PCMRI has gained prominence in recent years due to its unrestricted 3D anatomical coverage and minimal operator dependence [1, 17–19]. The connection between MRIbased complex blood flow analysis (such as, turbulence [20] and helical blood flow [21]) and MRIbased biomarkers (such as, wall shear stress [22–24] and pressure gradients [25–29]) with disease progression and diagnosis are active and promising areas of research. However, PCMRI provides limited spatial and temporal resolutions, which inevitably impacts the accuracy of MRIbased hemodynamic parameter estimates [30].
CFD is an alternative that has been used to predict flow patterns in various vascular geometries, including intracranial aneurysms [10], the thoracic aorta [11], and the carotid bifurcation, both in models [12–15] and in vivo [16]. The equations describing Newtonian fluid flow are solved numerically for specified boundary and initial condition data. Such approach provides arbitrarily high spatial and temporal resolution, and is in principle capable of estimating flow fields for arbitrarily complex vessel geometries. Absolute hemodynamic parameter estimates can be obtained directly from the highresolution flow fields produced by CFD, obviating the need for data smoothing or interpolation schemes.
The accuracy of conventional CFD routines hinges on many modeling assumptions that are not strictly true for in vivo vascular flow, including rigid vessel walls and uniform blood viscosity. Indeed, the proper choice of the underlying physics model is itself an open research question [10]. CFD predictions have so far shown variable agreement with PCMRI measurements [10, 11, 13, 15], and the applicability of CFD to robust flow estimation is still being debated.
The use of fluid mechanics techniques for improving PCMRI data is an active research field. Several algorithms from the literature use regularizations based on curl and divergence of the velocity field, which are associated with the irrotationality and incompressibility characteristics of the fluid flow, respectively. These include algorithms capable of improving streamlines [31, 32], and also algorithms for denoising the PCMRI data [33–36].
However, the solutions found by those methods do not necessarily satisfy the NavierStokes momentum equation. Other authors integrated pointbased measurements within a CFD solver, by adding a “force” term to the momentum equations that is proportional to the difference between predicted and measured velocities for a given grid point [37–39]. They used such an approach to integrate, respectively, Dopplerultrasound velocity measurements and cerebral aneurysm blood flow MRI data into a CFD solver.
More recently, a method to accelerate 4D cardiac flow MRI using CFD simulations was proposed. The image model was generated by integrating numerical blood flow simulations (calculated using openFOAM [40]) into the MRI imagereconstruction algorithm [41]. However, this approach can not guarantee that the fluid physics model (momentum and continuity equations) is satisfied by the reconstructed velocity map.
Conventional CFD uses medical imaging data only to specify the vessel geometry and the flow at the inlet and outlet boundaries, or other previously known initial and boundary condition data, and uses the assumed fluid physics model to find the solution in the interior of the calculation domain [28]. The goal of this work is to develop a more general, flexible and easy to implement numerical framework for harnessing additional PCMRI velocity measurements to construct a more robust and potentially more accurate CFDbased solution, considering PCMRI data as ground truth. The proposed method is able to make use of full (or incomplete) PCMRI measurements of one or more velocity components within the entire 3D volume. This is achieved through generalized Tikhonov regularization [42], obtaining a numerical solution that is close enough to the measured flow data; satisfies the fluid physics equations; reduces noise; and, in the clinical environment, can be used to reduce scan time.
Finally, this work is presented as a proof of concept of the CFD–MRI combined solver. All simulations herein were made using the finite volume method and SIMPLER algorithm in Cartesian grids, with unrealistic assumptions about the blood flow model, such as rigid walls and Newtonian viscosity. Nevertheless, the optimal numerical solution proposed in this work is general enough to be implemented: for any type of discretization method, such as finite differences, finite volume or finite element; for steady or unsteady flow; and, for any realistic physics model, such as nonNewtonian viscosity, elastic vessel walls, or slightly compressible flow. Even in the most realistic model, the discretization (finite differences, finite volume or finite elements) of the nonlinear set of differential equations produces a large and sparse system of linear equations, that forms the basis of the proposed numerical solution. The feasibility of the proposed approach is demonstrated on data from the carotid bifurcation of one healthy volunteer, and from a pulsatile carotid flow phantom. Two implementations of the regularized computational solution were evaluated and compared: one using only the PCMRI data corresponding to the main velocity component (z axis); and another, using PCMRI data corresponding to all three velocity components.
Theory
Blood flow model
The general model for fluid motion in 3D Euclidian space is given by the Navier–Stokes momentum equation [43]:
where \(\rho\) is the fluid density, \(\vec {\nu }=(u,v,w)\) is the flow velocity vector (u, v, and w are the velocity components associated with spatial axes x, y, and z, respectively), t is time, \(\nabla\) is the gradient operator, p is the hydrodynamic pressure, \(\hat{\tau }\) is the stress tensor, and \(\vec {b}\) represents the body forces acting on the fluid during the flow. The stress tensor represents the momentum transferred in virtue of the molecular motions and interactions within the fluid. It is a function of the scalar invariants of the strain rate tensor \(\hat{e}= (1/2) \left[ \nabla \vec {\nu }+(\nabla \vec {\nu })^{T} \right]\), where \(^T\) denotes the transpose operation. For an incompressible fluid, it can be written as \(\hat{\tau } = \mu (\hat{e}) \, \hat{e}\), where the scalar \(\mu (\hat{e})\) is the generalized Newtonian viscosity for a given \(\hat{e}\).
In this work, blood is modeled as a Newtonian, incompressible and isothermal fluid, with constant viscosity \(\mu\) and constant density \(\rho\). We are also assuming that there are no body forces acting on the blood flow. Then, the simplification of the general momentum equation, Eq. (1), provides our blood model equation [43]:
where \(\Delta\) is the Laplacian operator.
Since there are no sources of blood inside an artery, the flow field must also satisfy mass conservation [43], which is expressed by the continuity equation:
SIMPLER algorithm
Equations (2) and (3) must be solved for the unknown scalar field variables u, v, w, and p. Those equations are nonlinear and coupled, and attempting to solve them directly in one step is a formidable, if not impossible, task.
The semiimplicit method for pressurelinked equations revised (SIMPLER) algorithm [44] is a wellknown and established numerical routine for solving the momentum and continuity equations, subject to given boundary and initial conditions. It belongs to a class of algorithms capable of solving the nonlinear coupled fluid dynamic equations, which also includes the SIMPLE, SIMPLEC, and PISO algorithms [45]. For our purposes, SIMPLER’s major advantage is that it does not require an initial guess for the pressure field; instead an initial estimate of the velocity field is used.
The discretization of the momentum equation, Eq. (2), forms the basis of the iterative CFD routine, yielding three linear systems. Let N be the total number of grid points in the discrete 3D calculation domain, i.e., in a rectangular grid, \(N=N_{x}\cdot N_{y}\cdot N_{z}\), where \(N_{x}\), \(N_{y}\), and \(N_{z}\) represent the number of grid points along the x, y, and z axes, respectively. Then, for the nth iteration, we have:
where \({\mathbf {S}}_{u,n1}\), \({\mathbf {S}}_{v,n1}\), and \({\mathbf {S}}_{w,n1}\) are \(N\times N\) square heptadiagonal sparse matrices, each containing previous iteration information about all three velocity components, as well as the values of the density and viscosity constants; the three \(N\times 1\) column vectors \({\mathbf {u}}_{n}\), \({\mathbf {v}}_{n}\), and \({\mathbf {w}}_{n}\) store the current iteration of the u, v, and w velocity component values, respectively, associated with all grid points in the 3D calculation domain (Fig. 1); and each of the constant \(N\times 1\) column vectors \({\mathbf {f}}_{u,n1}\), \({\mathbf {f}}_{v,n1}\), and \({\mathbf {f}}_{w,n1}\) contains previous iteration information about all three velocity components, current iteration pressure difference values, and the physical parameters \(\rho\) and \(\mu\).
The main difficulty when attempting to predict the flow of an incompressible fluid is that there is no equation for pressure. Hence, a discretized pressure equation is deduced applying Eqs. (4), (5), and (6) in the discretized continuity equation, giving rise to a discretization of the following Poisson equation:
where \(\delta t\) is the time step, i.e. the time increment between iterations (this is further discussed in the next section). For a given time instant \(t=t_i\), convergence is achieved when \(\nabla \cdot \vec {\nu }=0\). The pressure field is updated at each iteration based on the current velocity field estimate, and hence does not appear explicitly in Eqs. (4), (5), and (6). It is important to note that u, v, and w values must be defined on regular grids, staggered by half a grid spacing (along the three directions) with respect to the grid on which pressure values are defined. This is to avoid nonphysical wiggle solutions for the pressure and velocity fields [44].
The steps of the SIMPLER algorithm, for a given time instant \(t = t_i\), are summarized in Algorithm 1. A detailed explanation of the discretization of the Navier–Stokes equations is provided in Refs. [44, 45].
Methods
Algorithm implementation
On constructing the numerical solution of the unsteady Navier–Stokes equation, we assume that a velocity field and the boundary conditions at a given time instant \(t = t_0\) are known. For this initial set of data, the numerical solution for the next time step \(t_1=t_0+\delta t\) is constructed, and it converges toward the solution when the continuity equation is satisfied (step 7 in Algorithm 1). Starting with the solution for \(t_1\), the same iterative procedure is repeated to obtain the solution for \(t_2 = t_1 + \delta t\), and so forth. In this manner, a timedependent flow field is computed.
In unsteady flow simulations of the Navier–Stokes equations by implicit numerical routines, the nature of the transient SIMPLER iterative procedures is equivalent to steadystate SIMPLER calculations applied, until convergence, for each time instant [45]. In other words, solving transient problems using SIMPLER is equivalent to solving successive steadystate problems. Moreover, steadystate calculations may be interpreted as pseudotransient solutions with spatiallyvarying time steps [45]. In other words, the steadystate solutions are, in practice, unsteady solutions, considering a virtual time step \(\delta t\) with fixed boundary conditions and initial data. This approach has been used by many recent authors [13, 16, 37–39]. While this was the numerical strategy adopted in this work, the approach proposed here can also be used for unsteady flow predictions.
The steadystate solution \(\vec {\nu }_{\infty }\) corresponding to a given cardiac phase (a temporal frame within the cardiac cycle) is calculated using the MRImeasured inlet and outlet velocities for that cardiac phase. CFD calculations begin with an initial guess for \(\vec {\nu }\), and simulations are carried forward in time until convergence, i.e.:
given a suficiently small \(\varepsilon\) (\(\left\ \cdot \right\\) denotes vector magnitude) and suficiently small time step \(\delta t\). Note that, here, time t is a simulationonly parameter, and is unrelated to time instants within the cardiac cycle. It is also not related with the iteration steps (n) of the SIMPLER algorithm (Algorithm 1), since the entire algorithm—with multiple iterations until convergence criterion \(\nabla \cdot \vec {\nu }=0\) is satisfied—is performed at each time instant t, until the convergence criterion shown in Eq. 8 is satisfied. At this point, \(\vec {\nu }_{\infty }\) is obtained. If multiple cardiac phases were to be reconstructed, then \(\vec {\nu }_{\infty }\) would be independently calculated for each cardiac phase.
Our implementation of the SIMPLER algorithm was validated with the bidimensional liddriven cavity flow problem, known in the literature as a benchmark for testing CFD algorithms [46–48]. All algorithms were implemented in Matlab (The MathWorks, Inc., Natick, MA, USA). Linear systems were solved using the biconjugate gradients stabilized method.
Proposed numerical solution
In this paper, we solve for a simulated velocity field, \(\vec {\nu }_{\infty } = (u_{\infty },v_{\infty },w_{\infty })\), that is close enough to the MRImeasured vector field \(\vec {\nu }_{\mathrm {mri}} = (u_{\mathrm {mri}},v_{\mathrm {mri}},w_{\mathrm {mri}})\), and satisfies the fluid dynamics equations, Eqs. (2) and (3).
Let M be the total number of voxels in the reconstructed \(\vec {\nu }_{\mathrm {mri}}\) 3D velocity field, i.e., \(M=M_{x}\cdot M_{y}\cdot M_{z}\), where \(M_{x}\), \(M_{y}\), and \(M_{z}\) represent the number of voxels along the x, y, and z axes, respectively. Consider \({\mathbf {u}}_{\mathrm {mri}}\), \({\mathbf {v}}_{\mathrm {mri}}\), and \({\mathbf {w}}_{\mathrm {mri}}\) as the stacked \(M \times 1\) column vectors with the PCMRI measurements. Since the numerical solution of the Navier–Stokes continuityequation system is based on the solution of linear systems, we propose that our numerical optimal solution is obtained by minimizing, for each velocity component, at iteration n, the following equations:
The first term on the right hand side of Eqs. (9)–(11) is related to the numerical solution of the Navier–Stokes continuity equations, and the second term is related to the comparison between the numerical solution and the PCMRI velocity field. The \(\mathbf {S}\) matrices and \(\mathbf {f}\) vectors are defined in Eqs. (4)–(6), but note that we dropped the “\(n1\)” subscripts for simplicity and clarity; these are updated by velocity and pressure values calculated in the previous iteration. Coefficients \(\lambda _u\), \(\lambda _v\), and \(\lambda _w\) are regularization factors, which weight consistance with PCMRI data against conformance with the momentum equations. Matrices \(\Gamma _u\), \(\Gamma _v\), and \(\Gamma _w\) are of size \(M\times N\), and model the blurring effects due to finite kspace coverage in PCMRI (this is further discussed below), while adjusting the number of points on the CFD grid in order to allow a comparison between \(\vec {\nu }_n\) and \(\vec {\nu }_\mathrm {mri}\). In this approach, the number of grid points in the CFD and MRI grids are not necessarily the same; we can use a finer grid in CFD than in MRI, for example. The optimal solutions for Eqs. (9)–(11) are straightforward [42], and given by
To understand the construction of the \(\Gamma\) matrices, consider that in the absence of noise, artifacts, and distortions, the MRImeasured vector field \(\vec {\nu }_{\mathrm {mri}} = (u_{\mathrm {mri}},v_{\mathrm {mri}},w_{\mathrm {mri}})\) is a blurred version of the true vector field \(\vec {\nu } = (u,v,w)\). For the u component, for example, we can write:
where \(*\) denotes convolution, and blurring kernel \(\psi _u(x,y,z)\) is the pointspread function associated with the kspace coverage that was used when measuring \(u_{\mathrm {mri}}\). Similarly, we can write \(v_{\mathrm {mri}} = v *\psi _v\) and \(w_{\mathrm {mri}} = w *\psi _w\). If all three PCMRI velocity components are measured using the same kspace coverage, then \(\psi _u = \psi _v = \psi _w\). For a 3DFT acquisition, these spatial blurring kernels are equal to
where \(\delta x\), \(\delta y\) and \(\delta z\) are the spatial resolutions of \(\vec {\nu }_{\mathrm {mri}}\) along the x, y, and z axis, respectively.
We want the CFDestimated vector field \(\vec {\nu }_{\infty } = (u_{\infty },v_{\infty },w_{\infty })\) to be an accurate representation of the true vector field \(\vec {\nu }\). If this is so, then we should expect \(u_{\mathrm {mri}} \approx u_{\infty } *\psi _u\), \(v_{\mathrm {mri}} \approx v_{\infty } *\psi _v\), and \(w_{\mathrm {mri}} \approx w_{\infty } *\psi _w\). The discretization of these equations yields three linear systems. Then, using the same notation introduced earlier, for the nth iteration of the CFD algorithm, we can write:
The coefficients of \(\Gamma _u\), \(\Gamma _v\), and \(\Gamma _w\) are calculated from \(\psi _u\), \(\psi _v\), and \(\psi _w\), respectively. If all three PCMRI velocity components are measured using the same kspace coverage, and reconstructed onto identical grids, then \(\Gamma _u = \Gamma _v = \Gamma _w\).
The MRIguided CFD estimate corresponding to one cardiac phase was calculated as a steadystate solution \(\vec {\nu }_{\infty }\). All three components of the PCMRI velocity field \(\vec {\nu }_\mathrm {mri}\) measured at the z positions at the boundaries of the calculation domain were used as inlet and outlet boundary conditions for that cardiac phase. Note that this steady state solution \(\vec {\nu }_{\infty }\) is the closest fit in the leastsquares sense to the direct PCMRI measurements that satisfy both momentum equation (Eq. 2) and continuity equation (Eq. 3). This is guaranteed by the fact that the optimal solutions Eqs. (12)–(14) are solved in each iteration of the SIMPLER algorithm (steps 2 and 4, in Algorithm 1), and by the convergence criterion (step 7).
In each of our experiments, all three PCMRI velocity components were measured using the same kspace coverage, and reconstructed onto identical grids. In the phantom experiments, we used the same grid size for both \(\vec {\nu }_{\mathrm {mri}}\) and \(\vec {\nu }_{\infty }\), because the phantom data were measured with high spatial resolution. In these experiments, \(\vec {\nu }_{\mathrm {mri}}\) was reconstructed without zeropadding, i.e., onto \(\delta x \times \delta y \times \delta z\) voxels, and the CFD grid points were defined at the center of each of \(\vec {\nu }_{\mathrm {mri}}\)’s voxels. Hence, \(\Gamma _u = \Gamma _v = \Gamma _w\) was defined as an identity matrix. In the in vivo experiments, \(\vec {\nu }_{\mathrm {mri}}\) was reconstructed using 2fold zeropadding along each of the spatial axes, since the data was acquired with low spatial resolution. Then, \(\Gamma\) was an \(N\times N\) symmetric matrix, with coefficients calculated from the point spread function \(\psi (x,y,z)\), defined in Eq. (16). This infinite support point spread function was truncated by multiplication with the box function
where \(\text{ rect }(w)=1\), if \(1\le w\le 1\), and \(\text{ rect }(w)=0\), otherwise.
Experimental setup: phantom demonstration
PCMRI data of a pulsatile carotid flow phantom (Phantoms by Design, Inc., Bothell, WA) (Fig. 2) were obtained with high spatial resolution and high signal–to–noise ratio, from four timeresolved 3DFT FGRE image volumes (three acquired each with a velocity encoding bipolar gradient on one of the three axes, and one without a bipolar gradient). The scan parameters were: \(0.5\times 0.5\times 1.0\) mm\(^3\) spatial resolution; FOV \(4.0\times 3.5\times 5.0\) cm\(^3\); TR 11.4 ms; flip angle 8.5\(^\circ\); temporal resolution 91.2 ms; VENC 50 cm/s; 40 min per scan; 9 NEX. The data were acquired on a GE Discovery MR750 3T system (50 mT/m and 200 T/m/s max gradient amplitude and slew rate), with a 32channel receiveonly head coil array (Nova Medical, Inc., Wilmington, MA, USA). The throughslab (z) axis was oriented along the S/I direction. The phantom’s pulse cycle was set to 60 bpm.
Only the temporal frame corresponding to peak flow was reconstructed. PCMRI velocity component maps \(u_{\mathrm {mri}},\) \(v_{\mathrm {mri}}\) and \(w_{\mathrm {mri}}\) were calculated using data from all channels of the receive coil array. The lumen was segmented by manually outlining the vessel borders from a stack of 2D axial images, obtained from the reconstructed 3D volume. A few voxels presented phasewrap artifacts; these voxels were manually identified and their velocities were corrected by adding \(2\pi\) to their values.
The combined solver calculations assumed fluid viscosity of \(\mu\) = 0.005 Pa s and density of \(\rho\) = 1100 kg/m\(^3\) (these values were provided by the phantom manufacturer). Calculations were performed with time step \(\delta t\) = 0.1 ms on a Cartesian grid of \(0.5\times 0.5 \times 1.0\) mm\(^3\) voxel size.
The CFD simulation domain was rectangular, of size \(32.5\times 9.0\times 41.0\,\,\mathrm {mm}^{3}\). Each iteration required about 10 seconds of computation time on an Intel Core i7 processor running at 2.8 GHz.
Three simulated steadystate velocity fields \(\vec {\nu }_{\infty }\) were obtained:

1.
using the conventional SIMPLER algorithm, i.e., not using the PCMRI data to constrain the CFD solution (\(\vec {\nu }_\mathrm {mri}\) was used only as inlet and outlet velocities for the geometry);

2.
using the velocity component associated with the main flow axis (z) measured with PCMRI (\(w_{\mathrm {mri}}\)) to constrain the CFD solution (u and v components were determined solely from the fluid physics model); and

3.
using all three velocity components measured with PCMRI (\(u_{\mathrm {mri}}\), \(v_{\mathrm {mri}}\), and \(w_{\mathrm {mri}}\)) to constrain the CFD solution.
The first approach is equivalent to making \(\lambda _w=\lambda _u=\lambda _v=0\); in the second approach, we used \(\lambda _w=1\) and \(\lambda _u=\lambda _v=0\); in the third approach, we used \(\lambda _u=\lambda _v=\lambda _w=1\). All three approaches used all three components of the PCMRI velocity field \(\vec {\nu }_\mathrm {mri}\) measured at the z positions at the boundaries of the calculation domain as inlet and outlet boundary conditions. The number of iterations until convergence for the above simulations was 89, 40 and 5 iterations, respectively.
Experimental setup: in vivo demonstration
PCMRI data of the carotid bifurcation of one healthy volunteer were obtained from four timeresolved 3DFT FGRE image volumes (three acquired each with a velocity encoding bipolar gradient on one of the three axes, and one without a bipolar gradient). The scan parameters were: \(1.0\times 1.0\times 2.5\) mm\(^3\) spatial resolution; FOV \(7.5\times 12.0\times 36.0\) cm\(^3\); TR 7.0 ms; flip angle 15\(^\circ\); temporal resolution 56 ms; VENC 160 cm/s; 7 min per scan; 1 NEX. The data were acquired on a GE Signa 3T EXCITE HD system (40 mT/m and 150 T/m/s max gradient amplitude and slew rate), with a 4channel neck receive coil array. The throughslab (z) axis was oriented along the S/I direction. The institutional review board of the University of Southern California approved the imaging protocols. The subject was screened for MRI risk factors and provided informed consent in accordance with institutional policy.
Only the cardiac phase corresponding to peak flow was reconstructed. PCMRI velocity component maps \(u_{\mathrm {mri}},\) \(v_{\mathrm {mri}}\) and \(w_{\mathrm {mri}}\) were calculated using data from only one channel of the receive coil array. Residual linear velocity offsets in each velocity component map (e.g., due to eddycurrents) were removed by performing a linear fit to manually defined 3D regions containing only stationary tissue. The lumen was segmented by manually outlining the vessel borders from a stack of 2D axial images, obtained from the reconstructed 3D volume.
The combined solver calculations assumed blood viscosity \(\mu\) = 0.0032 Pa s and density of \(\rho\) = 1060 kg/m\(^3\) [49]. Calculations were performed with time step \(\delta t\) = 0.25 ms on a Cartesian grid of \(0.50\times 0.50 \times 1.25\) mm\(^3\) voxel size. The CFD simulation domain was rectangular, of size \(30\times 37\times 125\) mm\(^3\) (the PCMRI data was cropped to match this grid size). Each iteration required about 180 s of computation time on an Intel Core i7 processor running at 2.8 GHz.
Three simulated steadystate velocity fields \(\vec {\nu }_{\infty }\) were obtained, using the same three approaches used in the phantom experiment. The number of iterations until convergence for the simulations was 1058, 190 and 6, respectively.
Quantitative evaluation
The CFDsimulated velocity fields were quantitatively compared with the PCMRI measurements by means of the signaltoerror ratio (SER). The SER measures the ratio between the energy of the signal and the energy of the estimation error. We considered the PCMRI velocity field, \(\vec {\nu }_{\mathrm {mri}} = (u_{\mathrm {mri}},v_{\mathrm {mri}},w_{\mathrm {mri}})\), as our groundtruth “signal”; consequently, the estimation error is the vector difference between the CFDestimated velocity field, \(\vec {\nu }_{\infty } = (u_{\infty },v_{\infty },w_{\infty })\), and the groundtruth field, \(\vec {\nu }_{\mathrm {mri}}\). Thus, the SER is calculated (in decibels) as:
where integers i, j, and k represent gridpoint indexes along the x, y, and z axes, respectively. Similarly, the SER was also calculated individually for each of the velocity components, as:
Using these SER values, the three CFD approaches—pure CFD, CFD driven by one PCMRI velocity component, and CFD driven by all three PCMRI velocity components—were quantitatively evaluated and compared.
Evaluation of denoising properties
Under our hypothesis, CFD simulations provide a smooth, noisefree flow field. Therefore, we expect that the proposed approach can be used as a denoising mechanism for PCMRI flow assessment. In order to verify the denoising effects of the combined solver, we added zeromean Gaussian noise with standard deviation 8 cm/s to the phantom’s measured velocity field, \(\vec {\nu }_{\mathrm {mri}}\).
This noisy flow field was used to constrain the CFD calculations, using the approach in which all three velocity components measured with PCMRI are used. In this experiment, we used \(\lambda _u=\lambda _v=\lambda _w=\lambda\); and four different values of \(\lambda\) were evaluated: \(5\times 10^{9}\), \(5\times 10^{8}\), \(5\times 10^{7}\), and \(5\times 10^{6}\). The SER between the proposed approach and the original PCMRI measurements was calculated, and compared with the SER of the noisy flow field. The pure CFD approach, in which the noisy PCMRI data are used only as inlet and outlet velocities for the geometry, was also evaluated (this is equivalent to making \(\lambda = 0\)).
The phantom data was used in this denoising experiment, because it was acquired using 9 NEX—which results in high signaltonoise ratio (SNR), while the in vivo data was acquired using only 1 NEX. The noise levels on the phantom’s measured velocity components—estimated as the standard deviation in regions of uniform mean velocity—are lower than 3 cm/s; while the SNR of the magnitude images exceeds 26 dB. The velocitytonoise ratio (VNR) [50, 51] for the u and w components reach 28 and 31 dB, respectively (the VNR for the v component was not calculated, because v is approximately null over the entire geometry).
Finally, in order to justify our denoising experiment, we analyze the noise distribution in PCMRI images. We note that, from a maximum likelihood perspective, Eqs. (9)–(11) assume that the PCMRI data is degraded by Gaussian noise. Under certain conditions, one can prove that velocity field noise in PCMRI satisfies a zeromean Gaussian distribution [52]. Therefore, the additive noise acting on the velocity fields can be assumed to be Gaussian distributed [27, 52]. Hence, the proposed minimization is wellsuited in terms of the MR noise distribution.
Results
Phantom demonstration
Figure 3 shows a qualitative velocitymap comparison between the PCMRI phantom measurements and the three simulated results. The PCMRI velocity field does not satisfy the continuity equation, since its divergence is nonzero within the lumen (Fig. 3a). The pure CFD solution produced a velocity field that satisfies the physical model, but is considerably smoother than the PCMRI measurements (Fig. 3b). Using one MRImeasured velocity component (\(w_{\mathrm {mri}}\)) to guide the CFD simulation resulted in a solution that is qualitatively more similar to the MRImeasured field, while still satisfying the continuity and momentum equations (Fig. 3c). Even better agreement was achieved when all three MRImeasured velocity components (\(u_{\mathrm {mri}}\), \(v_{\mathrm {mri}}\), and \(w_{\mathrm {mri}}\)) were used to guide the CFD simulation (Fig. 3d). These improvements can be also appreciated on a vector field visualization of the flow field over the entire tridimensional volume (Fig. 4).
Table 1 shows the SER between the phantom experiment results from each of the three CFD approaches, relative to the MRImeasured velocity field. The approach using only one PCMRI velocity component (\(w_{\mathrm {mri}}\)) to drive the CFD calculations (labeled “CFD + 1D”) provided a 11.09 dB improvement in SER relative to pure CFD (labeled simply “CFD”) with respect to the w component; however, the improvement was only 1.81 dB when considering all three components. The approach using all three PCMRI velocity components (\(u_{\mathrm {mri}}\), \(v_{\mathrm {mri}}\), and \(w_{\mathrm {mri}}\)) to drive the CFD calculations (labeled “CFD + 3D”) provided a 6.56 dB improvement in SER relative to pure CFD when considering all three components, while still providing a 8.02 dB improvement when considering only the w component.
Note that, for each of the CFD approaches, the SER was lower for the u and v components than it was for w (Table 1). This can be explained by the fact that the same VENC was used for measuring all three PCMRI velocity components (which implies similar noise levels), but the velocities along the z axis are considerably higher than those along the x and y axes (the energy of \(w_{\mathrm {mri}}\) was 11.42 dB higher than that of \(u_{\mathrm {mri}}\), and 15.72 dB higher than that of \(v_{\mathrm {mri}}\)). As a consequence, the SNR of \(w_{\mathrm {mri}}\) is substantially higher than that of \(u_{\mathrm {mri}}\) and \(v_{\mathrm {mri}}\). Thus, even if the energy of the absolute error between the CFDestimated velocities and the MRImeasured velocities was the same for the three components, \({\mathrm {SER}}_{w}\) would be higher than \({\mathrm {SER}}_{u}\) and \({\mathrm {SER}}_{v}\). Also, note that CFD approaches provide a smooth, (ideally) noisefree velocity field, but the SER was calculated with respect to a noisy PCMRI field. This means that the denoising properties of the CFD approaches could actually hurt the SER, especially if noise levels are relatively high when compared to the velocity values (this is the case for the u and v components). However, denoising is a desirable feature, so this does not necessarily indicate an unwanted result.
Figure 5 and Table 2 show the results of the experiment in which the denoising properties of the proposed approach were evaluated. CFD calculations were constrained by all three PCMRI velocity components, with added Gaussian noise. The combined solver improved the SER of each individual velocity component, for all different weight parameters we evaluated (Table 2). The CFD solution constrained by PCMRI using \(\lambda =5\times 10^{8}\) (Fig. 5e) provides a velocity field that is less noisy and visually more similar to the measured PCMRI velocity field (Fig. 5a) than the pure CFD solution obtained using the noisy PCMRI velocity field as boundary data (Fig. 5c). Using smaller values of \(\lambda\) (Fig. 5d) results in solutions that are closer to the pure CFD solution (Fig. 5c). Using larger values of \(\lambda\) (Figs. 5f, g) results in solutions that are closer to the noisy MRI data used to constrain the solution (Fig. 5b). These results can be appreciated quantitatively on Table 2. Relative to the noisy PCMRI velocity field, the overall (\(\vec {\nu }\)) SER gain was 1.23 dB for \(\lambda = 5\times 10^{9}\); 2.47 dB for \(\lambda = 5\times 10^{8}\); 2.14 dB for \(\lambda = 5\times 10^{7}\); and 1.23 dB for \(\lambda = 5\times 10^{6}\). Moreover, all constrained CFD solutions presented better quantitative results than pure CFD. The overall (\(\vec {\nu }\)) SER gain, relative to pure CFD, was 0.55 dB for \(\lambda = 5\times 10^{9}\); 1.79 dB for \(\lambda = 5\times 10^{8}\); 1.46 dB for \(\lambda = 5\times 10^{7}\); and 0.55 dB for \(\lambda = 5\times 10^{6}\).
These results illustrate the potential of the proposed numerical framework also as a denoising technique for PCMRI.
In vivo demonstration
Figure 6 provides a qualitative velocitymap comparison between the PCMRI in vivo measurements and the three simulated results. As in the phantom experiment, the PCMRI velocity field does not satisfy the continuity equation, since its divergence is nonzero within the lumen (Fig. 6a). The pure CFD solution produced a velocity field that satisfies the physical model, but differs considerably from the PCMRI measurements (Fig. 6b). Using one MRImeasured velocity component (\(w_{\mathrm {mri}}\)) to guide the CFD simulation resulted in a solution that is qualitatively more similar to the MRImeasured field, while still satisfying the continuity and momentum equations (Fig. 6c). Even better agreement was achieved when all three MRImeasured velocity components (\(u_{\mathrm {mri}}\), \(v_{\mathrm {mri}}\), and \(w_{\mathrm {mri}}\)) were used to guide the CFD simulation (Fig. 6d). These improvements can be also appreciated on a vector field visualization of the flow field over the entire tridimensional volume (Fig. 7).
Table 3 shows the SER between the in vivo experiment results from each of the three CFD approaches, relative to the MRImeasured velocity field. The approach using only one PCMRI velocity component (\(w_{\mathrm {mri}}\)) to drive the CFD calculations (labeled “CFD + 1D”) provided a 9.79 dB improvement in SER relative to pure CFD (labeled simply “CFD”) with respect to the w component; however, the SER was only 0.76 dB higher when considering all three components, since the agreement with PCMRI for the v component was made worst. The approach using all three PCMRI velocity components (\(u_{\mathrm {mri}}\), \(v_{\mathrm {mri}}\), and \(w_{\mathrm {mri}}\)) to drive the CFD calculations (labeled “CFD + 3D”) improved the SER for all three components, by 0.24, 0.59, and 6.00 dB, respectively; as a result, there was a 2.15 dB overall improvement in SER relative to pure CFD, and a 1.39 dB improvement relative to CFD guided only by \(w_{\mathrm {mri}}\). As in the phantom experiments, the SER was lower for the u and v components than it was for w, for all three CFD approaches. This can again be explained by the fact that the same VENC was used for all three velocity components, while the velocities along the z axis are considerably higher than those along the x and y axes (the energy of \(w_{\mathrm {mri}}\) was 29.76 dB higher than that of \(u_{\mathrm {mri}}\), and 27.09 dB higher than that of \(v_{\mathrm {mri}}\)), and is a direct (and desirable) consequence of the denoising properties of the proposed approach, as previously discussed.
Discussion
The proposed methodology uses the threedimensional SIMPLER algorithm with Cartesian uniform meshes to perform blood flow simulations under the influence of threedimensional MRImeasured velocity profiles. The combined MRI–CFD methodology attempts to correct the MRImeasured flow field, forcing it to satisfy the fluid mechanics equations. The choice of the SIMPER algorithm and Cartesian discretization was performed in order to facilitate implementation of the algorithm. We showed that the proposed technique provides better agreement with the PCMRI measurements than pure CFD simulations. We also showed that this MRIguided CFD approach can be used as a means of reducing noise in the PCMRI measurements. It can also be used to reduce computation time: when all three MRImeasured velocity components are used to guide the CFD simulations, only a few iterations are required for convergence, i.e., for finding the flow field that is the most similar to the PCMRI measurements (in the least squares sense) while satisfying the fluid mechanics equations.
We believe that the proposed method can also be used as a technique for reducing scan time. The proposed methodology allows the use of only part of the velocity measurements obtained with MRI to guide computational solutions by appropriately choosing the \(\Gamma\) matrices in Eqs. (12)–(14). In this paper, we used identical grids for MRI and CFD. Therefore, the \(\Gamma\) matrices were square with size \(N\times N\). In each of the following suggested approaches, at least one of the \(\Gamma _u\), \(\Gamma _v\) or \(\Gamma _w\) matrices may be of size \(M\times N\), with \(M<N\). Possible approaches for reducing scan time include: (1) acquiring the MRI data with reduced spatial resolution, and using the MRIguided CFD simulation to improve the spatial resolution; (2) measuring portion(s) of the volume (e.g., carotid inlet and outlets) with full spatial resolution, while measuring the rest of the volume with reduced spatial resolution; (3) acquiring only a few slices along the bifurcation, and using the MRIguided CFD simulation to fill in the gaps; (4) measuring one or two velocity components with full spatial resolution, while measuring the other component(s) with reduced spatial resolution; (5) measuring one or two velocity components across the entire volume, while measuring the other component(s) for only a few slices; and (6) any combination of these approaches. We plan to explore these ideas in future studies.
It is well known that blood is a nonNewtonian fluid, therefore its viscosity is not uniform. While there exists many constitutive models and studies in the literature regarding the nonNewtonian rheological properties of blood [43, 53, 54], no goldstandard constitutive model exists, and the assumption of constant whole blood viscosity is common practice [13, 16, 33, 37, 55, 56]. Therefore, we used the Newtonian blood flow model (constant whole blood viscosity), which greatly simplified the implementation.
It is often desirable to estimate the wall shear rate at the carotid bifurcation. Determining wall shear rates from CFD simulation results would require using highly refined meshes around the neighborhood of the vessel wall. Our implementation of the SIMPLER algorithm does not allow using spatiallyvarying grid spacing, and the mesh is uniform all over the integration domain. Using a very fine grid over the entire integration domain is possible, in principle. However, this would drastically increase the computational complexity, and could make the proposed methodology impractical in a clinical environment, for example.
Finally, the proposed approach does not take the effects of vessel wall elasticity into consideration. While the pulsatile carotid flow phantom we used has rigid tube walls, the human carotid vessel wall is generally elastic. While there exists blood flow CFD simulation methods that incorporate elastic wall effects [57, 58], the assumption of rigid vessel walls is another common practice [13, 16, 33, 37, 55]. The SIMPLER algorithm implemented in this work uses a Cartesian uniform mesh, and does not allow the use of the elastic wall models.
Both wall shear stress calculations and elastic vessel walls could be properly addressed with an implementation using finite elements—which would allow an easier adaptation of the mesh near the wall and also allow simulating the effects of fluid–structure interaction. Despite the fact that the implementation of a finiteelement solver would be substantially more complex than our implementation, all the proposed methodology described in this proof of concept paper can still be applied in the same fashion, since the problem of solving the set of differential equations in a finite element discretization is replaced by a sparse system of linear equations similar to the ones obtained in this work.
In future works, this MRIguided CFD methodology for unsteady flow will be implemented on FreeFem++,^{Footnote 1} a partial differential equation solver capable of solving the Navier–Stokes equations using finite elements [59]. FreeFem++ allows access to the linear systems that can be modified in order to make use of the principles introduced here. With this approach, we expect more general, higherquality solutions, allowing the calculation of biomarkers, such as wall shear stress, since FreeFem++ can handle different types of triangular finite elements, large varieties of linear system solvers, and automatic mesh adaptation. Note that there are other free software that could be used to reproduce this methodology, such as openFOAM,^{Footnote 2} which uses finite volume discretization. Both FreeFem++ and openFOAM allow modifications of the linear systems, and have CFD solvers already implemented, which could facilitate the implementation of the methodology proposed in this study [40, 59].
All the assumptions and simplifications disscussed above (Newtonian blood flow, imprecise geometry, noncompliant walls) contribute nonlinearly to the differences observed between CFD solutions and MRImeasured velocity fields. Using the MRImeasured velocity field to constrain the CFD solution indirectly addresses these simplifications, and provides a more realistic CFD solution. However, we are still unable to clearly identify the main factors responsible for the disagreement between MRImeasured and CFDcomputed velocity fields. Nevertheless, an implementation using a more robust solver, as proposed above, could improve on these limitations and further reduce the gap between MRImeasured and CFDcomputed results.
Conclusion
We have proposed a framework for obtaining flow field estimates that are influenced by both PCMRI measurements and a fluid physics model. The results showed that the proposed technique provides better agreement with the PCMRI measurements than pure CFD simulations, and has reduced computation time (faster convergence). MRIguided CFD can be used to correct the MRImeasured flow field, forcing it to satisfy the fluid mechanics equations. It can also be used as a means of reducing noise in the PCMRI measurements, and has potential as a method for reducing scan time.
The proposed framework offers a general approach to in vivo blood flow assessment, that is complementary to improvements in PCMRI acquisition and reconstruction techniques, and can be applied to the study and diagnosis of a broad range of cardiovascular flow mapping applications.
Notes
References
 1.
Pelc NJ, Herfkens RJ, Shimakawa A, Enzmann DR. Phasecontrast cine magnetic resonance imaging. Magn Reson Q. 1991;7:229–54.
 2.
Gonzalez EG, Carvalho JLA. Does phase contrast MRI provide the mean velocity of the spins within a voxel? In: Proceedings of the International Society for Magnetic Resonance in Medicine, vol 22. Milan; 2014. p. 2480.
 3.
Sigfridsson A, Petersson S, Carlhäll CJ, Ebbers T. Fourdimensional flow MRI using spiral acquisition. Magn Reson Med. 2012;68:1065–73. doi:10.1002/mrm.23297.
 4.
Steinman DA, Taylor CA. Flow imaging and computing: large artery hemodynamics. Ann Biomed Eng. 2005;33:1704–9. doi:10.1007/s1043900587722.
 5.
Giddens DP, Zarins CK, Glagov S. The role of fluid mechanics in the localization and detection of atherosclerosis. J Biomech Eng. 1993;115:588–94. doi:10.1115/1.2895545.
 6.
Steinman DA. Imagebased computational fluid dynamics: a new paradigm for monitoring hemodynamics and atherosclerosis. Curr Drug Targets Cardiovasc Hematol Disord. 2004;4:183–97. doi:10.2174/1568006043336302.
 7.
Lei M, Kleinstreuer C, Truskey GA. Numerical investigation and prediction of atherogenic sites in branching arteries. J Biomech Eng. 1995;117:350–7. doi:10.1115/1.2794191.
 8.
Hyun S, Kleinstreuer C, Archie JP. Hemodynamics analyses of arterial expansions with implications to thrombosis and restenosis. Med Eng Phys. 2000;22:13–27. doi:10.1016/S13504533(00)000060.
 9.
Longest PW, Kleinstreuer C. Comparison of blood particle deposition models for nonparallel flow domains. J Biomech. 2003;36:421–30. doi:10.1016/S00219290(02)004347.
 10.
Boussel L, Rayz V, Martin A, AcevedoBolton G, Lawton MT, Higashida R, Smith WS, Young WL, Saloner D. Phasecontrast magnetic resonance imaging measurements in intracranial aneurysms in vivo of flow patterns, velocity fields, and wall shear stress: comparison with computational fluid dynamics. Magn Reson Med. 2009;61(2):409–17. doi:10.1002/mrm.21861.
 11.
Canstein C, Cachot P, Faust A, Stalder AF, Bock J, Frydrychowicz A, Kuffer J, Hennig J, Markl M. 3D MR flow analysis in realistic rapidprototyping model systems of the thoracic aorta: comparison with in vivo data and computational fluid dynamics in identical vessel geometries. Magn Reson Med. 2008;59(3):535–46. doi:10.1002/mrm.21331.
 12.
Milner JS, Moore JA, Rutt BK, Steinman DA. Hemodynamics of human carotid artery bifurcations: computational studies with models reconstructed from magnetic resonance imaging of normal subjects. J Vasc Surg. 1998;28:143–56. doi:10.1016/S07415214(98)702101.
 13.
Papathanasopoulou P, Zhao S, Köhler U, Robertson MB, Long Q, Hoskins P, Xu XY, Marshall I. MRI measurement of timeresolved wall shear stress vectors in a carotid bifurcation model, and comparison with CFD predictions. J Magn Reson Imaging. 2003;17:153–62. doi:10.1002/jmri.1220.
 14.
Thomas JB, Milner JS, Rutt BK, Steinman DA. Reproducibility of imagebased computational fluid dynamics models of the human carotid bifurcation. Ann Biomed Eng. 2003;31:132–41. doi:10.1114/1.1540102.
 15.
Marshall I, Zhao S, Papathanasopoulou P, Hoskins P, Xu Y. MRI and CFD studies of pulsatile flow in healthy and stenosed carotid bifurcation models. J Biomech. 2004;37:679–87. doi:10.1016/j.jbiomech.2003.09.032.
 16.
Long Q, Xu XY, Ariff B, Thom SA, Hughes AD, Stanton AV. Reconstruction of blood flow patterns in a humancarotid bifurcation: a combined CFD and MRI study. J Magn Reson Imaging. 2000;11:299–311.
 17.
Pelc NJ, Sommer FG, Li KC, Brosnan TJ, Herfkens RJ, Enzmann DR. Quantitative magnetic resonance flow imaging. Magn Reson Q. 1994;10:125–47.
 18.
Markl M, Chan FP, Alley MT, Wedding KL, Draney MT, Elkins CJ, Parker DW, Wicker R, Taylor CA, Herfkens RJ, Pelc NJ. Timeresolved threedimensional phasecontrast MRI. J Magn Reson Imaging. 2003;17:499–506. doi:10.1002/jmri.10272.
 19.
Harloff A, Albrecht F, Spreer J, Stalder AF, Bock J, Frydrychowicz A, Schllhorn J, Hetzel A, Schumacher M, Hennig J, Markl M. 3D blood flow characteristics in the carotid artery bifurcation assessed by flowsensitive 4D MRI at 3T. Magn Reson Med. 2009;61:65–74. doi:10.1002/mrm.21774.
 20.
Dyverfeldt P, Gårdhagen R, Sigfridsson A, Karlsson M, Ebbers T. On MRI turbulence quantification. Magn Reson Imaging. 2009;27:913–22. doi:10.1016/j.mri.2009.05.004.
 21.
Morbiducci U, Ponzini R, Rizzo G, Cadioli M, Esposito A, Montevecchi FM, Redaelli A. Mechanistic insight into the physiological relevance of helical blood flow in the human aorta: an in vivo study. Biomech Model Mechanobiol. 2011;10:339–55. doi:10.1007/s1023701002382.
 22.
Tsuji T, Suzuki J, Shimamoto R, Yamazaki T, Nakajima T, Nagai R, Komatsu S, Ohtomo K, ToyoOka T, Omata M. Vector analysis of the wall shear rate at the human aortoiliac bifurcation using cine MR velocity mapping. Am J Roentgenol. 2002;178:995–9. doi:10.2214/ajr.178.4.1780995.
 23.
Stalder AF, Russe MF, Frydrychowicz A, Bock J, Markl JHM. Quantitative 2D and 3D phase contrast MRI: optimized analysis of blood flow and vessel wall parameters. Magn Reson Med. 2008;60:1218–31. doi:10.1002/mrm.21778.
 24.
Dyverfeldt P, Sigfridsson A, Knutsson H, Ebbers T. A novel MRI framework for the quantification of any moment of arbitrary velocity distributions. Magn Reson Med. 2011;65(3):725–31. doi:10.1002/mrm.22649.
 25.
Ebbers T, Wigstrm L, Bolger AF, Wranne B, Karlsson M. Noninvasive measurement of timevarying threedimensional relative pressure fields within the human heart. J Biomech Eng. 2002;124:288–93. doi:10.1115/1.1468866.
 26.
Thompson RB, McVeigh ER. Fast measurement of intracardiac pressure differences with 2D breathhold phasecontrast MRI. Magn Reson Med. 2003;49:1056–66. doi:10.1002/mrm.10486.
 27.
Wang Y, Amini AA. Integrable pressure gradients via harmonicsbased orthogonal projection. In: Christensen GE, Sonka M, editors. Information processing in medical imaging: 19th international conference. Glenwood Springs; 2005. p. 431–42.
 28.
Nayak KS, Nielsen JF, Bernstein MA, Markl M, Botnar RM, Gatehouse PD, Saloner D, Lorenz C, Wen H, Hu BS, Epstein FH, Oshinski JN, Raman SV. Cardiovascular magnetic resonance phase contrast imaging. J Cardiovasc Magn Reson. 2015;17:71. doi:10.1186/s1296801501727.
 29.
Donati F, Figueroa CA, Smith NP, Lamata P, Nordslettena DA. Noninvasive pressure difference estimation from PCMRI using the workenergy equation. Med Image Anal. 2015;26:159–72. doi:10.1016/j.media.2015.08.012.
 30.
Harloff A, Nussbaumer A, Bauer S, Stalder AF, Frydrychowicz A, Weiller C, Hennig J, Markl M. In vivo assessment of wall shear stress in the atherosclerotic aorta using flowsensitive 4D MRI. Magn Reson Med. 2010;63(6):1529–36. doi:10.1002/mrm.22383.
 31.
Buonocore M. Algorithms for improving calculated streamlines in 3D phase contrast angiography. Magn Reson Med. 1994;31(1):22–30. doi:10.1002/mrm.1910310104.
 32.
Fatouraee N, Amini AA. Regularization of flow streamlines in multislice phasecontrast MR imaging. IEEE Trans Med Imaging. 2003;22(6):699–709. doi:10.1109/TMI.2003.814786.
 33.
Song SM, Napel S, Glover GH, Pelc NJ. Noise reduction in threedimensional phase contrast MR velocity measurements. J Magn Reson Imaging. 1993;3:587–96. doi:10.1002/jmri.1880030407.
 34.
Verdugo AMP, Mura J, Uribe S. Enforcing divergence free to velocity data from 4D flow MR images. In: Proceedings of International Society for Magnetic Resonance in Medicine, vol 22. Milan;2014. p. 2485.
 35.
Ong F, Uecker M, Tariq U, Hsiao A, Alley M, Vasanawala S, Lustig M. Compressed sensing 4D flow reconstruction using divergencefree wavelet transform. In: Proceedings of international society for magnetic resonance in medicine, vol 22. Milan;2014. p. 0326.
 36.
Bostan E, Lefkimmiatis S, Vardoulis O, Stergiopulos N, Unser M. Improved variational denoising of flow fields with application to phasecontrast MRI data. IEEE Signal Process Lett. 2014;22(6):762–6. doi:10.1109/LSP.2014.2369212.
 37.
Liu L, Funamoto K, Hayase T. Numerical experiment for ultrasonicmeasurementintegrated simulation of developed laminar pipe flow using axisymmetric model. J Biomech Sci Eng. 2008;3:101–15. doi:10.1007/s1043900895197.
 38.
Kato T, Funamoto K, Hayase T, Sone S, Kadowaki H, Shimazaki T, Jibiki T, Miyama K, Liu L. Development and feasibility study of a twodimensional ultrasonicmeasurementintegrated blood flow analysis system for hemodynamics in carotid arteries. Med Biol Eng Comput. 2014;52(11):933–43. doi:10.1007/s1151701411933.
 39.
Funamoto K, Suzuki Y, Hayase T, Kosugi T, Isoda H. Numerical validation of MRmeasurementintegrated simulation of blood flow in a cerebral aneurysm. Ann Biomed Eng. 2009;37(6):1105–16. doi:10.1007/s104390099689y.
 40.
Weller HG, Tabor G, Jasak H, Fureby C. A tensorial approach to computational continuum mechanics using objectoriented techniques. Computers Phys. 1998;12(6):620–31. doi:10.1063/1.168744.
 41.
Christodoulou AG, Ramb R, Menza M, Hennig J, Liang ZP. 4d flow imaging incorporating a fluid dynamics model. In: Proceedings of the international society for magnetic resonance in medicine, Toronto, vol 23. 2015. p. 2735.
 42.
Seo JK, Woo EJ. Nonlinear inverse problems in imaging. 1st ed. Chichester: Willey; 2013.
 43.
Bird R, Armstrong R, Hassager O. Dynamics of polymeric liquids: fluid mechanics, vol. 1, 2nd ed. New York: Willey; 1987.
 44.
Patankar SV. Numerical heat transfer and fluid flow. 1st ed. New York: Hemisphere Publishing Corporation; 1980.
 45.
Versteeg H, Malalasekera W. An introduction to computational fluid dynamics: the finite. 2nd ed. Glasgow: PrenticeHall; 2007.
 46.
Griebel M, Dornseifer T, Neunhoeffer T. Numerical simulation in fluid dynamics. 1st ed. Philadelphia: Society for Industrial and Applied Mathematics; 1998.
 47.
Ghia U, Ghia K, Shin C. Highresolutions for incompressible flow using the Navier–Stokes equations and a multigrid method. J Comput Phys. 1982;48:387–411. doi:10.1016/00219991(82)900584.
 48.
Gupta M, Kalita J. A new paradigm for solving Navier–Stokes equations: streamfunctionvelocity formulation. J Comput Phys. 2005;207:52–68. doi:10.1016/j.jcp.2005.01.002.
 49.
Westerhof N, Stergiopulos N, Noble M. Snapshots of hemodynamics: an aid for clinical research and graduate education. 1st ed. New York: Springer; 2005.
 50.
Papaharilaou Y, Doorly DJ, Sherwin SJ. Assessing the accuracy of twodimensional phasecontrast MRI measurements of complex unsteady flows. J Magn Reson Imaging. 2001;14:714–23. doi:10.1002/jmri.10008.
 51.
Ha H, Kim GB, Kweon J, Kim YH, Kim N, Yang DH, Lee SJ. MultiVENC acquisition of fourdimensional phasecontrast MRI to improve precision of velocity field measurement. Magn Reson Med. 2015. doi:10.1002/mrm.25715.
 52.
Gudbjartsson H, Patz S. The rician distribution of noisy MRI data. Magn Reson Med. 1995;34(6):910–4. doi:10.1002/mrm.1910340618.
 53.
Kim S. A study of nonnewtonian viscosity and yield stress of blood in a scanning capillarytube rheometer. PhD thesis, Drexel University. 2002.
 54.
Kim S, Namgung B, Ong P, Cho Y, Chun K, Lim D. Determination of rheological properties of whole blood with a scanning capillarytube rheometer using constitutive models. J Mech Sci Technol. 2009;23(6):1718–26. doi:10.1007/s1220600904206.
 55.
Steinman DA, Thomas JB, Ladak HM, Milner JS, Rutt BK, Spence JD. Reconstruction of carotid bifurcation hemodynamics and wall thickness using computational fluid dynamics and MRI. Magn Reson Med. 2002;47(1):149–59. doi:10.1002/mrm.10025.
 56.
Carvalho JLA, Nielsen JF, Nayak KS. Feasibility of in vivo measurement of carotid wall shear rate using spiral Fourier velocity encoded MRI. Magn Reson Med. 2010;63:1537–47. doi:10.1002/mrm.22325.
 57.
Perktold K, Rappitsch G. Computer simulation of local blood flow and vessel mechanics in a compliant carotid artery bifurcation model. J Biomech. 1995;28(7):845–56. doi:10.1016/00219290(95)952738.
 58.
Figueroa CA, VignonClementela IE, Jansenc KE, Hughesd TJR, Taylorb CA. A coupled momentum method for modeling blood flow in threedimensional deformable arteries. Computer Methods Appl Mech Eng. 2006;195:5685–706. doi:10.1016/j.cma.2005.11.011.
 59.
Hecht F. New development in freefem++. J Numer Math. 2012;20(3–4):251–65. doi:10.1515/jnum20120013.
Authors’ contributions
VCR: CFD–MRI code implementation, phantom data acquisition, analyzed and interpreted the data, drafted the manuscript. JFN: assisted CFD–MRI code implementation, in vivo and phantom data aquisition/reconstruction, drafted and critically revised the manuscript, supervised the research. KSN: in vivo data aquisition/reconstruction, critically revised the manuscript. JLC: design of the study, analyzed and interpreted the data, critically revised the manuscript, supervised the research. All authors read and approved the final manuscript.
Acknowledgements
This work was supported in part by the Brazilian Coordination for the Improvement of Higher Education Personnel (CAPES), scholarship Grant number 8195137.
Competing interests
The authors declare that they have no competing interests.
Author information
Affiliations
Corresponding author
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Rispoli, V.C., Nielsen, J.F., Nayak, K.S. et al. Computational fluid dynamics simulations of blood flow regularized by 3D phase contrast MRI. BioMed Eng OnLine 14, 110 (2015). https://doi.org/10.1186/s1293801501047
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s1293801501047
Keywords
 Computational fluid dynamics
 Phase contrast
 Magnetic resonance imaging