# Computational fluid dynamics simulations of blood flow regularized by 3D phase contrast MRI

- Vinicius C. Rispoli
^{1, 2}Email author, - Jon F. Nielsen
^{3}, - Krishna S. Nayak
^{4}and - Joao L. A. Carvalho
^{1}

**14**:110

https://doi.org/10.1186/s12938-015-0104-7

© Rispoli et al. 2015

**Received: **25 August 2015

**Accepted: **16 November 2015

**Published: **26 November 2015

## Abstract

### Background

Phase contrast magnetic resonance imaging (PC-MRI) is used clinically for quantitative assessment of cardiovascular flow and function, as it is capable of providing directly-measured 3D velocity maps. Alternatively, vascular flow can be estimated from model-based computation fluid dynamics (CFD) calculations. CFD provides arbitrarily high resolution, but its accuracy hinges on model assumptions, while velocity fields measured with PC-MRI 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 proof-of-concept numerical procedure for constructing a simulated flow field that is influenced by both direct PC-MRI measurements and a fluid physics model, thereby taking advantage of both the accuracy of PC-MRI 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 PC-MRI 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 PC-MRI 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 PC-MRI measurements than CFD alone, and converges faster, while closely satisfying the fluid dynamics equations. For the implementation that provided the best results, the signal-to-error ratio (with respect to the PC-MRI 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 PC-MRI measurements.

## Keywords

## 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, model-independent velocity mapping using phase contrast magnetic resonance imaging (PC-MRI) [1–3] or Doppler ultrasound, and model-based computational fluid dynamics (CFD) calculations [4–16]. Among the direct methods, PC-MRI has gained prominence in recent years due to its unrestricted 3D anatomical coverage and minimal operator dependence [1, 17–19]. The connection between MRI-based complex blood flow analysis (such as, turbulence [20] and helical blood flow [21]) and MRI-based 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, PC-MRI provides limited spatial and temporal resolutions, which inevitably impacts the accuracy of MRI-based 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 high-resolution 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 PC-MRI 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 PC-MRI 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 PC-MRI data [33–36].

However, the solutions found by those methods do not necessarily satisfy the Navier-Stokes momentum equation. Other authors integrated point-based 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, Doppler-ultrasound 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 image-reconstruction 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 PC-MRI velocity measurements to construct a more robust and potentially more accurate CFD-based solution, considering PC-MRI data as ground truth. The proposed method is able to make use of full (or incomplete) PC-MRI 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 non-Newtonian 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 PC-MRI data corresponding to the main velocity component (*z* axis); and another, using PC-MRI data corresponding to all three velocity components.

## Theory

### Blood flow model

*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}\).

### SIMPLER algorithm

Equations (2) and (3) must be solved for the unknown scalar field variables *u*, *v*, *w*, and *p*. Those equations are non-linear and coupled, and attempting to solve them directly in one step is a formidable, if not impossible, task.

The semi-implicit method for pressure-linked equations revised (SIMPLER) algorithm [44] is a well-known 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 non-linear 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.

*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

*n*-th iteration, we have:

*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,n-1}\), \({\mathbf {f}}_{v,n-1}\), and \({\mathbf {f}}_{w,n-1}\) contains previous iteration information about all three velocity components, current iteration pressure difference values, and the physical parameters \(\rho\) and \(\mu\).

*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 non-physical wiggle solutions for the pressure and velocity fields [44].

## 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 time-dependent 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 steady-state SIMPLER calculations applied, until convergence, for each time instant [45]. In other words, solving transient problems using SIMPLER is equivalent to solving successive steady-state problems. Moreover, steady-state calculations may be interpreted as pseudo-transient solutions with spatially-varying time steps [45]. In other words, the steady-state 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.

*t*is a simulation-only 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 lid-driven 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 MRI-measured 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).

*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 PC-MRI measurements. Since the numerical solution of the Navier–Stokes continuity-equation 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:

*u*component, for example, we can write:

*x*,

*y*, and

*z*axis, respectively.

*n*th iteration of the CFD algorithm, we can write:

The MRI-guided CFD estimate corresponding to one cardiac phase was calculated as a steady-state solution \(\vec {\nu }_{\infty }\). All three components of the PC-MRI 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 least-squares sense to the direct PC-MRI 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).

### Experimental setup: phantom demonstration

*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. PC-MRI 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 phase-wrap 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.

- 1.
using the conventional SIMPLER algorithm, i.e., not using the PC-MRI 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 PC-MRI (\(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 PC-MRI (\(u_{\mathrm {mri}}\), \(v_{\mathrm {mri}}\), and \(w_{\mathrm {mri}}\)) to constrain the CFD solution.

*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

PC-MRI data of the carotid bifurcation of one healthy volunteer were obtained from four time-resolved 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 4-channel neck receive coil array. The through-slab (*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. PC-MRI 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 eddy-currents) 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 PC-MRI 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 steady-state 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

*i*,

*j*, and

*k*represent grid-point indexes along the

*x*,

*y*, and

*z*axes, respectively. Similarly, the SER was also calculated individually for each of the velocity components, as:

### Evaluation of denoising properties

Under our hypothesis, CFD simulations provide a smooth, noise-free flow field. Therefore, we expect that the proposed approach can be used as a denoising mechanism for PC-MRI flow assessment. In order to verify the denoising effects of the combined solver, we added zero-mean 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 PC-MRI 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 PC-MRI measurements was calculated, and compared with the SER of the noisy flow field. The pure CFD approach, in which the noisy PC-MRI 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 signal-to-noise 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 velocity-to-noise 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 PC-MRI images. We note that, from a maximum likelihood perspective, Eqs. (9)–(11) assume that the PC-MRI data is degraded by Gaussian noise. Under certain conditions, one can prove that velocity field noise in PC-MRI satisfies a zero-mean 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 well-suited in terms of the MR noise distribution.

## Results

### Phantom demonstration

Signal-to-error ratio (in dB) between the phantom experiment results from each of the three CFD approaches—pure CFD (labeled “CFD”); CFD driven by one PC-MRI velocity component (labeled “CFD + 1D”); and CFD driven by all three PC-MRI velocity components (labeled “CFD + 3D”), relative to the MRI-measured velocity field

CFD | CFD + 1D | CFD + 3D | |
---|---|---|---|

\({\mathrm {SER}}_{u}\) | 2.97 | 4.16 (\(\uparrow\)) | 6.74 (\(\uparrow\)) |

\({\mathrm {SER}}_{v}\) | −0.25 | \(-\)0.30 (\(\approx\)) | 2.03 (\(\uparrow\)) |

\({\mathrm {SER}}_{w}\) | 5.44 | 16.53 (\(\uparrow \uparrow\)) | 13.46 (\(\uparrow\)) |

\({\mathrm {SER}}_{\vec {\nu }}\) | 6.57 | 8.38 (\(\uparrow\)) | 13.13 (\(\uparrow \uparrow\)) |

Table 1 shows the SER between the phantom experiment results from each of the three CFD approaches, relative to the MRI-measured velocity field. The approach using only one PC-MRI 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 PC-MRI 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 PC-MRI 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 CFD-estimated velocities and the MRI-measured 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) noise-free velocity field, but the SER was calculated with respect to a noisy PC-MRI 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.

Signal-to-error ratio (in dB) between noisy and original PC-MRI measurements; and between the MRI-guided CFD estimates and the original PC-MRI measurements

Noisy | CFD | CFD+3D | CFD+3D | CFD+3D | CFD+3D | |
---|---|---|---|---|---|---|

PC-MRI | \({\lambda =0}\) | \({\lambda =5\times 10^{-9}}\) | \({\lambda =5\times 10^{-8}}\) | \({\lambda =5\times 10^{-7}}\) | \({\lambda =5\times 10^{-6}}\) | |

\({\mathrm {SER}}_{u}\) | 0.80 | 2.66 | 2.87 | 3.22 | 2.34 | 1.75 |

\({\mathrm {SER}}_{v}\) | −2.46 | −0.50 | −0.41 | −0.39 | −1.18 | −1.67 |

\({\mathrm {SER}}_{w}\) | 4.64 | 5.03 | 5.82 | 7.20 | 5.88 | 5.06 |

\({\mathrm {SER}}_{\vec {\nu }}\) | 5.71 | 6.39 | 6.94 | 8.18 | 7.85 | 6.94 |

These results illustrate the potential of the proposed numerical framework also as a denoising technique for PC-MRI.

### In vivo demonstration

*w*component; however, the SER was only 0.76 dB higher when considering all three components, since the agreement with PC-MRI for the

*v*component was made worst. The approach using all three PC-MRI 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.

Signal-to-error ratio (in dB) between the in vivo experiment results from each of the three CFD approaches—pure CFD (labeled “CFD”); CFD driven by one PC-MRI velocity component (labeled “CFD + 1D”); and CFD driven by all three PC-MRI velocity components (labeled “CFD + 3D”), relative to the MRI-measured velocity field

CFD | CFD + 1D | CFD + 3D | |
---|---|---|---|

\({\mathrm {SER}}_{u}\) | 0.51 | 0.50 (\(\approx\)) | 0.75 (\(\uparrow\)) |

\({\mathrm {SER}}_{v}\) | −1.51 | \(-\)1.87 (\(\downarrow\)) | \(-\)0.92 (\(\uparrow\)) |

\({\mathrm {SER}}_{w}\) | 4.74 | 14.53 (\(\uparrow \uparrow\)) | 10.74 (\(\uparrow\)) |

\({\mathrm {SER}}_{\vec {\nu }}\) | 4.13 | 4.89 (\(\uparrow\)) | 6.28 (\(\uparrow \uparrow\)) |

## Discussion

The proposed methodology uses the three-dimensional SIMPLER algorithm with Cartesian uniform meshes to perform blood flow simulations under the influence of three-dimensional MRI-measured velocity profiles. The combined MRI–CFD methodology attempts to correct the MRI-measured 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 PC-MRI measurements than pure CFD simulations. We also showed that this MRI-guided CFD approach can be used as a means of reducing noise in the PC-MRI measurements. It can also be used to reduce computation time: when all three MRI-measured 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 PC-MRI 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 MRI-guided 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 MRI-guided 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 non-Newtonian fluid, therefore its viscosity is not uniform. While there exists many constitutive models and studies in the literature regarding the non-Newtonian rheological properties of blood [43, 53, 54], no gold-standard 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 spatially-varying 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 finite-element 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 MRI-guided CFD methodology for unsteady flow will be implemented on FreeFem++,^{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, higher-quality 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,^{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, non-compliant walls) contribute non-linearly to the differences observed between CFD solutions and MRI-measured velocity fields. Using the MRI-measured 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 MRI-measured and CFD-computed velocity fields. Nevertheless, an implementation using a more robust solver, as proposed above, could improve on these limitations and further reduce the gap between MRI-measured and CFD-computed results.

## Conclusion

We have proposed a framework for obtaining flow field estimates that are influenced by both PC-MRI measurements and a fluid physics model. The results showed that the proposed technique provides better agreement with the PC-MRI measurements than pure CFD simulations, and has reduced computation time (faster convergence). MRI-guided CFD can be used to correct the MRI-measured flow field, forcing it to satisfy the fluid mechanics equations. It can also be used as a means of reducing noise in the PC-MRI 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 PC-MRI acquisition and reconstruction techniques, and can be applied to the study and diagnosis of a broad range of cardiovascular flow mapping applications.

## Declarations

### 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 8195-13-7.

### Competing interests

The authors declare that they have no competing interests.

**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.

## Authors’ Affiliations

## References

- Pelc NJ, Herfkens RJ, Shimakawa A, Enzmann DR. Phase-contrast cine magnetic resonance imaging. Magn Reson Q. 1991;7:229–54.Google Scholar
- 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.Google Scholar
- Sigfridsson A, Petersson S, Carlhäll C-J, Ebbers T. Four-dimensional flow MRI using spiral acquisition. Magn Reson Med. 2012;68:1065–73. doi:10.1002/mrm.23297.View ArticleGoogle Scholar
- Steinman DA, Taylor CA. Flow imaging and computing: large artery hemodynamics. Ann Biomed Eng. 2005;33:1704–9. doi:10.1007/s10439-005-8772-2.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- Steinman DA. Image-based 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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/S1350-4533(00)00006-0.View ArticleGoogle Scholar
- Longest PW, Kleinstreuer C. Comparison of blood particle deposition models for non-parallel flow domains. J Biomech. 2003;36:421–30. doi:10.1016/S0021-9290(02)00434-7.View ArticleGoogle Scholar
- Boussel L, Rayz V, Martin A, Acevedo-Bolton G, Lawton MT, Higashida R, Smith WS, Young WL, Saloner D. Phase-contrast 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.View ArticleGoogle Scholar
- Canstein C, Cachot P, Faust A, Stalder AF, Bock J, Frydrychowicz A, Kuffer J, Hennig J, Markl M. 3D MR flow analysis in realistic rapid-prototyping 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.View ArticleGoogle Scholar
- 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/S0741-5214(98)70210-1.View ArticleGoogle Scholar
- Papathanasopoulou P, Zhao S, Köhler U, Robertson MB, Long Q, Hoskins P, Xu XY, Marshall I. MRI measurement of time-resolved 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.View ArticleGoogle Scholar
- Thomas JB, Milner JS, Rutt BK, Steinman DA. Reproducibility of image-based computational fluid dynamics models of the human carotid bifurcation. Ann Biomed Eng. 2003;31:132–41. doi:10.1114/1.1540102.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- Pelc NJ, Sommer FG, Li KC, Brosnan TJ, Herfkens RJ, Enzmann DR. Quantitative magnetic resonance flow imaging. Magn Reson Q. 1994;10:125–47.Google Scholar
- Markl M, Chan FP, Alley MT, Wedding KL, Draney MT, Elkins CJ, Parker DW, Wicker R, Taylor CA, Herfkens RJ, Pelc NJ. Time-resolved three-dimensional phase-contrast MRI. J Magn Reson Imaging. 2003;17:499–506. doi:10.1002/jmri.10272.View ArticleGoogle Scholar
- 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 flow-sensitive 4D MRI at 3T. Magn Reson Med. 2009;61:65–74. doi:10.1002/mrm.21774.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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/s10237-010-0238-2.View ArticleGoogle Scholar
- Tsuji T, Suzuki J, Shimamoto R, Yamazaki T, Nakajima T, Nagai R, Komatsu S, Ohtomo K, Toyo-Oka 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- Ebbers T, Wigstrm L, Bolger AF, Wranne B, Karlsson M. Non-invasive measurement of time-varying three-dimensional relative pressure fields within the human heart. J Biomech Eng. 2002;124:288–93. doi:10.1115/1.1468866.View ArticleGoogle Scholar
- Thompson RB, McVeigh ER. Fast measurement of intracardiac pressure differences with 2D breath-hold phase-contrast MRI. Magn Reson Med. 2003;49:1056–66. doi:10.1002/mrm.10486.View ArticleGoogle Scholar
- Wang Y, Amini AA. Integrable pressure gradients via harmonics-based orthogonal projection. In: Christensen GE, Sonka M, editors. Information processing in medical imaging: 19th international conference. Glenwood Springs; 2005. p. 431–42.Google Scholar
- 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/s12968-015-0172-7.View ArticleGoogle Scholar
- Donati F, Figueroa CA, Smith NP, Lamata P, Nordslettena DA. Non-invasive pressure difference estimation from PC-MRI using the work-energy equation. Med Image Anal. 2015;26:159–72. doi:10.1016/j.media.2015.08.012.View ArticleGoogle Scholar
- 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 flow-sensitive 4D MRI. Magn Reson Med. 2010;63(6):1529–36. doi:10.1002/mrm.22383.View ArticleGoogle Scholar
- Buonocore M. Algorithms for improving calculated streamlines in 3-D phase contrast angiography. Magn Reson Med. 1994;31(1):22–30. doi:10.1002/mrm.1910310104.View ArticleGoogle Scholar
- Fatouraee N, Amini AA. Regularization of flow streamlines in multislice phase-contrast MR imaging. IEEE Trans Med Imaging. 2003;22(6):699–709. doi:10.1109/TMI.2003.814786.View ArticleGoogle Scholar
- Song SM, Napel S, Glover GH, Pelc NJ. Noise reduction in three-dimensional phase contrast MR velocity measurements. J Magn Reson Imaging. 1993;3:587–96. doi:10.1002/jmri.1880030407.View ArticleGoogle Scholar
- 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.Google Scholar
- Ong F, Uecker M, Tariq U, Hsiao A, Alley M, Vasanawala S, Lustig M. Compressed sensing 4D flow reconstruction using divergence-free wavelet transform. In: Proceedings of international society for magnetic resonance in medicine, vol 22. Milan;2014. p. 0326.Google Scholar
- Bostan E, Lefkimmiatis S, Vardoulis O, Stergiopulos N, Unser M. Improved variational denoising of flow fields with application to phase-contrast MRI data. IEEE Signal Process Lett. 2014;22(6):762–6. doi:10.1109/LSP.2014.2369212.View ArticleGoogle Scholar
- Liu L, Funamoto K, Hayase T. Numerical experiment for ultrasonic-measurement-integrated simulation of developed laminar pipe flow using axisymmetric model. J Biomech Sci Eng. 2008;3:101–15. doi:10.1007/s10439-008-9519-7.View ArticleGoogle Scholar
- Kato T, Funamoto K, Hayase T, Sone S, Kadowaki H, Shimazaki T, Jibiki T, Miyama K, Liu L. Development and feasibility study of a two-dimensional ultrasonic-measurement-integrated blood flow analysis system for hemodynamics in carotid arteries. Med Biol Eng Comput. 2014;52(11):933–43. doi:10.1007/s11517-014-1193-3.View ArticleGoogle Scholar
- Funamoto K, Suzuki Y, Hayase T, Kosugi T, Isoda H. Numerical validation of MR-measurement-integrated simulation of blood flow in a cerebral aneurysm. Ann Biomed Eng. 2009;37(6):1105–16. doi:10.1007/s10439-009-9689-y.View ArticleGoogle Scholar
- Weller HG, Tabor G, Jasak H, Fureby C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Computers Phys. 1998;12(6):620–31. doi:10.1063/1.168744.View ArticleGoogle Scholar
- 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.Google Scholar
- Seo JK, Woo EJ. Nonlinear inverse problems in imaging. 1st ed. Chichester: Willey; 2013.View ArticleGoogle Scholar
- Bird R, Armstrong R, Hassager O. Dynamics of polymeric liquids: fluid mechanics, vol. 1, 2nd ed. New York: Willey; 1987.Google Scholar
- Patankar SV. Numerical heat transfer and fluid flow. 1st ed. New York: Hemisphere Publishing Corporation; 1980.MATHGoogle Scholar
- Versteeg H, Malalasekera W. An introduction to computational fluid dynamics: the finite. 2nd ed. Glasgow: Prentice-Hall; 2007.Google Scholar
- Griebel M, Dornseifer T, Neunhoeffer T. Numerical simulation in fluid dynamics. 1st ed. Philadelphia: Society for Industrial and Applied Mathematics; 1998.View ArticleGoogle Scholar
- Ghia U, Ghia K, Shin C. High-resolutions for incompressible flow using the Navier–Stokes equations and a multigrid method. J Comput Phys. 1982;48:387–411. doi:10.1016/0021-9991(82)90058-4.MATHView ArticleGoogle Scholar
- Gupta M, Kalita J. A new paradigm for solving Navier–Stokes equations: streamfunction-velocity formulation. J Comput Phys. 2005;207:52–68. doi:10.1016/j.jcp.2005.01.002.MATHMathSciNetView ArticleGoogle Scholar
- Westerhof N, Stergiopulos N, Noble M. Snapshots of hemodynamics: an aid for clinical research and graduate education. 1st ed. New York: Springer; 2005.Google Scholar
- Papaharilaou Y, Doorly DJ, Sherwin SJ. Assessing the accuracy of two-dimensional phase-contrast MRI measurements of complex unsteady flows. J Magn Reson Imaging. 2001;14:714–23. doi:10.1002/jmri.10008.View ArticleGoogle Scholar
- Ha H, Kim GB, Kweon J, Kim Y-H, Kim N, Yang DH, Lee SJ. Multi-VENC acquisition of four-dimensional phase-contrast MRI to improve precision of velocity field measurement. Magn Reson Med. 2015. doi:10.1002/mrm.25715.
- Gudbjartsson H, Patz S. The rician distribution of noisy MRI data. Magn Reson Med. 1995;34(6):910–4. doi:10.1002/mrm.1910340618.View ArticleGoogle Scholar
- Kim S. A study of non-newtonian viscosity and yield stress of blood in a scanning capillary-tube rheometer. PhD thesis, Drexel University. 2002.Google Scholar
- Kim S, Namgung B, Ong P, Cho Y, Chun K, Lim D. Determination of rheological properties of whole blood with a scanning capillary-tube rheometer using constitutive models. J Mech Sci Technol. 2009;23(6):1718–26. doi:10.1007/s12206-009-0420-6.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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.View ArticleGoogle Scholar
- 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/0021-9290(95)95273-8.View ArticleGoogle Scholar
- Figueroa CA, Vignon-Clementela IE, Jansenc KE, Hughesd TJR, Taylorb CA. A coupled momentum method for modeling blood flow in three-dimensional deformable arteries. Computer Methods Appl Mech Eng. 2006;195:5685–706. doi:10.1016/j.cma.2005.11.011.MATHView ArticleGoogle Scholar
- Hecht F. New development in freefem++. J Numer Math. 2012;20(3–4):251–65. doi:10.1515/jnum-2012-0013.MATHMathSciNetGoogle Scholar