Computational Fluid Dynamics and Its Impact on Flow Measurements Using Phase-Contrast MR-Angiography

Rationale and Objectives: Computational fluid dynamic (CFD) simulations are discussed with respect to their potential for quality assurance of flow quantification using commercial software for the evaluation of magnetic resonance phase contrast angiography (PCA) data. Materials and Methods: Magnetic resonance phase contrast angiography data was evaluated with the Nova software. CFD simulations were performed on that part of the vessel system where the flow behavior was unexpected or non-reliable. The CFD simulations were performed with in-house written software. Results: The numerical CFD calculations demonstrated that under reasonable boundary conditions, defined by the PCA velocity values, the flow behavior within the critical parts of the vessel system can be correctly reproduced. Conclusion: CFD simulations are an important extension to commercial flow quantification tools with regard to quality assurance.


Introduction
Magnetic resonance phase contrast angiography meanwhile is an established technique to measure the velocity and the direction of flow [1].The major drawbacks of this approach are often related to misaligned encoding gradients (venc), a too low SNR or high susceptibility gradients nearby (small) vessels.In each case the relative phase are not unique and the absolute velocity values are wrong.At this point the correct answer to the flow in the critical regions can only be derived on the base of physical models and computational flow dynamics (CFD).Almost all CFD-problems are associated with the solution of the Navier-Stokes equations [2,3], named after Claude-Louis Navier and George Gabriel Stokes, which describe the motion of fluid substances.The calculations were performed on incompressible fluids within human blood vessels.The idea is to selectively model the relevant vessel system and to include the known velocity values from the PCA measurements within the big vessels as boundary conditions.These velocity conditions are sufficient in most cases because the velocities tune in to the given pressure conditions.The intrinsic symmetry of the vessels furthermore can be used to reduce the dimensionality of the problem from three to two, which is necessary to use the software on commercial computers with limited working memory and CPU power-only the relative orientation of the vessels and their diameters (e.g.due to stenoses) have to be considered.
The spectrum of potential applications span the simulation of pressure and force distributions near the neck of aneurysms, the velocity and direction of flow within complex vessel systems with or without stenoses and pressure and vorticity distributions near atherosclerotic plaques.
In this technical note we will concentrate on a few representative examples in order to demonstrate the importance of quality assurance in addition to commercial flow evaluation software.

Flow Measurement
Magnetic resonance phase contrast angiography measurements of flow were performed on a 3 Tesla scanner.The FLASH sequence parameters were as follows: Base resolution 256, Phase resolution 60%, Flip angle 25 deg, TR 113.20 ms, TE 4.66 ms, Voxel size: 0.9 × 0.5 × 4.0 mm, TA: 1:06 min, 1st Signal/Mode Pulse/Retro, calculated phases 12, segments 5, slices 1.Each slice was planned on the base of a time-of-flight (TOF) measurement using the Nova software (http://www.vassolinc.com),which guarantees the exact orthogonal positioning relative to the course of the vessel.The flow measurements were performed for a number of selected vessels that define a part of the simulated vessel system.Figure 1 shows a typical Nova rendered 3D-presentation of the vessel system of a TOF measurement.The velocities are measured at the positions given by the small slice cartoons perpendicularly oriented to the vessel.The most important aspects which are essential for the choice of these locations are the distance to bifurcations and stenoses because often turbulences may occur in their immediate vicinity.Nova currently provides no option for the separate presentation of systolic and diastolic phases but merely shows the mean volumetric flow rate values on the diagrams (Figure 1(b)).The mean velocity values, necessary for the simulations, are not shown on the Nova diagram, but are stored in a text file.

Preparation of CFD
2D-projections from individual TOF maximum intensity projection maps are co-registered to a template (see e.g.http://commons/Wikimedia.org/wiki/Template:Other_ver sions/Circle_of_Willis, binary mask in Figure 2(a)) while preserving the individual properties (e.g.stenoses, malformations, relative orientation).The final mask was used to define the coordinates and nodes essential for the mesh2d Matlab program, which calculates the nodes and elements (Figure 2(b)) necessary for the numerical simulations.

Numerical Algorithms
The numerical code for the flow simulations is an in-house development, written in Matlab.The program simulates the motion of incompressible fluids via numerical solutions of the 2D, unsteady Navier-Stokes equations.The CFD-program is a vertex centered Finite Volume (FV) code that uses the median dual mesh to form the control volumes about each vertex by connecting the centroids of each cell to their edge midpoints.The time and spatial linearizations are dealt with separately.Fully 2nd order piecewise linear reconstructions are used to develop sparse operators, which approximate the continuous differential terms.These methods can be found in detail in [4,5] described by terms like Green's theorem and its application to face averaged gradients (Laplacian, unweighted least squares) or upwind piece-wise linear extrapolation with the HQUICK slope limiter (non-linear part).
The time integration of the equations is based on the well-known fractional step method [5,6].The non-staggered pressure/velocity arrangement was shown to support the common mesh scale pressure oscillation and hence the standard 2nd order incremental pressure correction method [7,8] was modified to incorporate Rhie-Chow stabilization [9].This reduced the order of the time integration to 1st/2nd.Unlike many other Navier-Stokes solvers the CFD-program makes use of a complete LU factorizetion for the Poisson pressure equation.Instead of using standard iterative approaches like preconditioned biconjugate gradients stabilized method (bicgstab (Matlab)) we decided on including the fill reducing UMFPACK algorithm [10], an approach which was shown to be far more efficient in this case.The momentum equations are solved in a semi-implicit manner.The non-linear terms are evaluated explicitly via a 2nd order TVD Runge-Kutta method [11].In order to solve the linear systems that arise due to the implicit treatment of the viscous terms, the bicgstab method was favored.Minimization of errors associated with the numerical solution was managed according to the concepts of Roache et al. [12].Copyright © 2012 SciRes.OJMI

CPU Parameters
The CFD simulations were performed on a personal computer with Intel ® Xeon processor, 2.4 GHz, 12 GB RAM, 4 CPUs, NVIDIA Quadro FX 580, Video RAM 512 MB, GPU frequency 450 MHz, Ubuntu 10.04.

Program Features
The mesh2d program of Matlab was used for the 2D mesh (triangles) generation.This type of mesh allows the automatic treatment of complex geometries.The mesh2d algorithm optimizes the mesh quality which is important to guarantee correct results.The CFD program solves the Navier-Stokes equations for the velocity and pressure fields of incompressible fluids.The time step size and the stability of the integration are adjustable.The velocity values measured by PCA were used as the boundary conditions (at the nodes).The program includes a module that calculates the [x,y] forces applied to a boundary due to pressure and skin friction effects and can be used in this example to calculate the lift and drag forces on an aneurysm neck.

Patient
High grade stenosis of the left ICA and occlusion of the right ICA.

Results
The total time needed for the preparation and simulation is in the range of 3 min which is actually feasible for emergency cases.The computer, however, should at least have 4 GB working memory and processors of the new generation (Pentium IV, Xeon) with at least 2.4 GHz. Figure 2(a) shows a typical bitmap of the circle-ofwillis which is then transformed to a 2D-mesh with nodes and elements (Figure 2(b)) using the mesh2d Matlab program.As a typical example for the application of the CFD-simulations, the data of a patient before and after stenting of the left ICA was used.The velocity values at the boundaries were given by the PCA measurements.The simulations of the corresponding vessel system before (Figures 3(a)-(c)) and after stenting (Figures 4(a)-(c 1.
In order to demonstrate the potential of the CFD simulation technique, an example for such an (artificial) aneurysm is given in Figure 5, for which the velocity (Figure 5   Copyright © 2012 SciRes.OJMI   The absolute values are not relevant in this case but merely the relative pressure and force enhancements at the neck of the aneurysm which finally may lead to a disrupture.The step from the TOF-MIP to a two-dimensional bitmap is currently not fully automated but necessitates manual intervention by an experienced radiologist.

Discussion
In future developments, a separate evaluation of different heart phases will be possible as well.This deficit, however, does not lower the impact of the results, especially because the flow directions before and after stenting do not depend on the heart phase.
In summarize, the intention of this work was to demonstrate, that the simulation of model based fluid dynamics provides important clinical information which is currently not available in most commercial software packages for flow quantification.

Figure 1 .
Figure 1.A typical 3D-plot (a) of a NOVA flow planning system, illustrating the rendered vessel system according to the TOF measurement.The velocities are measured at the positions given as small slice cartoons perpendicularly oriented to the vessel (here the right ACA).In (b) one possible schematic Nova diagram of a part of the circle-of-willis system is shown, which contains information about the vessel names, the direction of flow (green arrow) and volumetric flow rate values (ml/min).
)) confirm the change in direction within the LPCOM vessel (Figure 4(b)) with respect to the situation before stenting (Figure 3(b)): it is opposite to the flow within the RPCOM (Figures 3(c) and 4(c)).The associated Nova diagrams are shown in Figure 3(d) (before stent) and Figure 4(d) (after stent).The measured mean velocity values are shown in Table

Figure 2 .Figure 3 .
Figure 2. The bitmap of the circle-of-willis vessel system (a) and the mesh nodes and elements for a part of it (b), calculated by the mesh2d Matlab program.The mesh serves as the base for the CFD simulations.Potential boundaries which may have non-zero in-or outflow are marked with red cycles.

Figure 4 .
Figure 4. Simulated CFD data of the circle-of-willis vessel system, shown in Figure 2, if the boundary conditions are given by the velocity values of a real PCA measurement after stenting: note the opposite flow directions in the LPCOM (b) and RPCOM (c) and the change in direction within the LPCOM vessel (Figure 4(b)) with respect to the situation before stenting (Figure 3(b)) which coincides with the results of the Nova data (d).

Figure 5 .
Figure 5. CFD simulations of an artificial aneurysm.(a) The pressure distribution (b) and the forces (in dependence of time [sec]) (c) at the marked position (red arrow in a)) demonstrate, that under special conditions the local load can achieve critical values, followed by a disrupture of the aneurysm neck.
are shown if physiologically reasonable boundary conditions are presumed.At the border of the area, indicated by the red arrow, the forces (Figure 5(c)) in x-(blue) and y-direction (red) are shown in dependence of time.
CFD simulations were used to verify unexpected flow patterns that occur after stenting in the LPCOM vessel and to demonstrate the critical pressure and force distributions at the neck of an artificial aneurysm.The data shows that CFD simulations are feasible on comer-cially available PCs if the simulations are restricted to the essential vessel systems and if intrinsic symmetries are used to reduce the dimensionality of the problem.CFD data may also provide clinically useful information regarding post-stenotic vessel wall alteration and progresssion of atherosclerotic disease.More generally, wall shear stress is believed to play an important role in the evolution of atherosclerosis and other arterial diseases, such as aneurysms.At the current stage only vorticity but no turbulence model is implemented.Further CFD program options are different boundary conditions such as to extrapolate the boundary values of the pressure from the interior nodes or just to fix the value of the pressure at the outflow to zero.