**Assessment of Carotid Flow Using Magnetic Resonance Imaging and Computational Fluid Dynamics**

Vinicius C. Rispoli1, Joao L. A. Carvalho1, Jon F. Nielsen2 and Krishna S. Nayak3 <sup>1</sup>*Universidade de Brasília* <sup>2</sup>*University of Michigan* <sup>3</sup>*University of Southern California* <sup>1</sup>*Brazil* 2,3*USA*

#### **1. Introduction**

512 Fluid Dynamics, Computational Modeling and Applications

Shah, P.K. (1997). Inflammation, Metalloproteinases, and Increase Proteolysis: an Emmerging Paradigm in Aortic Aneurysm. *Circulation*, Vol. 96, pp. 2115-2117. Shojima, M., Oshima, M., Takagi, K., Torii, R., Hayakawa, M., Katada, K., Morita, A. &

Shojima, M., Oshima, M., Takagi, K., Torii, R., Nagata, K., Shirouzu, I., Morita, A. & Kirino,

Tang, B.T., Cheng, C.P., Draney, M.T., Wilson, N.M., Tsao, P.S., Herfkens, R.J. & Taylor, C.A.

Taylor, C.A. & Figueroa, C.A. (2009). Patient-Specific Modeling of Cardiovascular Mechanics. *The Annual Review of Biomedical Engineering*, Vol. 11, pp. 109-134. Thiriet, M. (2008). *Biology and Mechanics of Blood Flows: Part II – Mechanics and Medical* 

Thubrikar, M. J. (2007). *Vascular Mechanics and Pathology*, Springer, ISBN-10: 0-387-33816-0,

Toth, M., Nadasy, G.L., Nyary, I., Kerenyi, T., Orosz, M., Molnárka, G. & Monos, E. (1998).

Utter, G. & Rossman, J.S. (2007). Numerical Simulation of Saccular Aneurysm

Vorp DA, Raghavan ML, Webster MW.(1998 ) Mechanical wall stress in abdominal aortic aneurysm: influence of diameter and asymmetry. *J Vasc Surg*;27:632-639. Vorp, D.A. (2007). Biomechanics of Abdominal Aortic Aneurysm. *Journal of Biomechanics*,

Wernig, F. & Xu, Q. (2002). Mechanical Stress-Induced Apoptosis in the Cardiovascular

Zamir, M. (2005). *The Physics of Coronary Blood Flow*, Springer, ISBN 978-0387-25297-1, New

Zhang, B., Fugleholm, K., Day, L.B., Ye, S., Weller, R.O. & Day, I.N. (2003). Molecular

Sterically Inhomogenenous Viscoelastic Behavior of Human Saccular Cerebral

Hemodynamics: Influence of Morphology on Rupture Risk. *Journal of Biomechancis*,

System. *Progress in Biophysisc & Molecular Biology*, Vol. 78, Issues 2-3, (February-

Pathogenesis of Subarachnoid Haemorrhage. *The International Journal of* 

*Aspects*, Springer, ISBN 978-0-387-74848-1, New York, USA.

Aneurysms. *Journal of Vascular Reserach*, vol. 35, pp. 345-355.

*Biochemistry & Cell Biology*, Vol. 35, pp. 1341-1360.

*Stroke*, Vol. 35, 2500-2505.

pp. H668-H676.

New York, USA.

Vol. 40, 2716-2722.

Vol. 40, pp. 1887-1902.

April), pp. 105-137.

York, USA.

*Journal of Epidemiology*, Vol. 154, pp. 236-244.

Kirino, T. (2004). Magnitude and Role of Wall Shear Stress on Cerebral Aneurysm: Computational Fluid Dynamic Study of 20 Middle Cerebral Artery Aneurysms.

T. (2005). Role of the Bloodstream Impacting Force and the Local Pressure Elevation in the Rupture of Cerebral Aneurysms. *Stroke*, Vol. 36, pp. 1933-1938. Singh, K., Bonaa, K.H., Jacobsen, B.K., Bjork, L. & Solberg, S. (2001). Prevalence of and Risk

Factors for Abdominal Aortic Aneurysms in a Population-Based Study. *American* 

(2006). Abdominal Aortic Hemodynamics in Young Healthy Adults at Rest and During Lower Limb Exercise: Quantification Using Image-Based Computer Modeling. *American Journal of Physiology. Heart and Circulatory Physiology*, Vol. 291,

> Knowledge of blood flow patterns in the human body is a critical component in cardiovascular disease research and diagnosis. Carotid atherosclerosis, for example, refers to the narrowing of the carotid arteries. One symptom of atherosclerosis is abnormal flow. The carotid arteries supply blood to the brain, so early detection of carotid stenosis may prevent thrombotic stroke. The current clinical gold standard for cardiovascular flow measurement is Doppler ultrasound. However, evaluation by ultrasound is inadequate when there is fat, air, bone, or surgical scar in the acoustic path, and flow measurement is inaccurate when the ultrasound beam cannot be properly aligned with the axis of flow (Hoskins, 1996; Winkler & Wu, 1995). Two alternative approaches to 3D flow assessment are currently available to the researcher and clinician: (i) direct, model-independent velocity mapping using flow-encoded magnetic resonance imaging (MRI), and (ii) model-based computational fluid dynamics (CFD) calculations.

> MRI is potentially the most appropriate technique for addressing all aspects of cardiovascular disease examination. MRI overcomes the acoustical window limitations of ultrasound, potentially allowing flow measurements to be obtained along any direction, and for any vessel in the cardiovascular system. MRI measurements are also less operator-dependent than those of Doppler ultrasound, and the true direction of flow can generally be precisely measured. The current gold standard for MRI flow quantification is phase contrast (PC) (O'Donnell, 1985). However, PC suffers from partial-volume effects when a wide distribution of velocities is contained within a single voxel (Tang et al., 1993). This is particularly problematic when flow is turbulent and/or complex (e.g., flow jets due to stenosis) or at the interface between blood and vessel wall (viscous sublayer). This issue is typically addressed by increasing the spatial resolution, which dramatically affects the signal-to-noise ratio (SNR) and increases the scan time. Therefore, PC may be inadequate for estimating the peak velocity of stenotic flow jets and for assessing wall shear rate.

a million) will point north, creating a small net magnetization *M*0, which is strong enough to be detected. The net magnetization is proportional to the strength of the *B*<sup>0</sup> field, so MRI

<sup>515</sup> Assessment of Carotid Flow Using Magnetic

Other important components of the scanner are the gradient coils. There are typically three gradient coils (*Gx*, *Gy*, and *Gz*), that produce an intentional perturbation in the *B*<sup>0</sup> field when turned on ("played"). This perturbation varies linearly along each spatial direction (*x*, *y* and *z*), such that no perturbation exists at the iso-center of the magnet when these gradients are used. In the presence of an external magnetic field, the spins rotate about the axis of that field. *B*<sup>0</sup> is (approximately) spatially uniform, so all spins initially rotate at the same frequency (the Larmor frequency), *ω* = *γB*0, where *γ* is the gyromagnetic ratio (*γ* = 42.6 MHz/Tesla for hydrogen protons). However, when any of the gradients is played, the magnetic field becomes spatially varying, and so does the rotation frequency of the spins. Therefore, *Gx*, *Gy*, and *Gz* are used to frequency-encode (or phase-encode) spatial position along the *x*, *y* and *z* directions,

The final major component of the MR scanner is the radio-frequency (RF) coil. This is used to transmit an RF "excitation" pulse to the body, and also to receive the frequency-encoded signal from the "excited" portion of the body. In practice, independent coils may be used for transmission and reception. The RF pulse is typically modulated to the Larmor frequency. While *B*<sup>0</sup> is aligned with the *z*-axis (by definition), *B*1, which is a very weak magnetic field associated with the RF pulse, is aligned with the *x*-axis (also by definition). When the RF pulse is played, some of the spins which are in resonance with the RF pulse (i.e., rotating at the RF pulse's frequency) will now begin to rotate around the *x*-axis (thus the name magnetic resonance). This tilts the net magnetization towards the *x*-*y* plane, and the net magnetization

The RF pulse is typically designed to have a somewhat rectangular profile in Fourier domain, centered at the modulation frequency (e.g., a modulated windowed sinc). This implies that the RF pulse in fact contains a certain range of frequencies, thus all spins rotating within that range become "excited", or tilted towards the *x*-axis. So, by playing gradient(s) of appropriate amplitude, and designing the RF pulse accordingly, one can excite only a thin slice of the body, which correspond to the region containing all spins that are in resonance with the RF pulse's range of frequencies. Excitation profiles other than "slices" may also be obtained (e.g., a pencil beam, or cylindrical excitation (Hu et al., 1993)), by designing an appropriate gradient/pulse

When the RF pulse is turned off, *Mxy* begins to rotate (at the Larmor frequency) around the *z*-axis, as the net magnetization begins to realign with *B*0. This rotating magnetization generates an oscillating signal, which can be detected by the receive coil. The frequency content of the received signal can be used to obtain spatial information about the excited portion of the body. In order to frequency-encode spatial information, gradients are also played during signal acquisition. These are called readout gradients. For imaging a slice perpendicular to the *z*-axis (an axial image), *Gz* is played during excitation (for slice-selection), and *Gx* and *Gy* are played during acquisition. These can be switched, for acquiring sagittal or coronal images, or all three gradients may be used during both excitation and acquisition to

When the readout gradients are played, the acquired signal at a particular time instant corresponds to the sum of sinusoidal signals generated by spins located at different regions of the body, each rotating at different frequencies corresponding to their spatial locations. If

scanners with stronger magnetic fields (e.g., 3 Tesla) provide higher SNR.

Resonance Imaging and Computational Fluid Dynamics

will now have a component in the *x*-*y* plane,*Mxy*.

respectively.

combination.

image oblique planes.

An alternative approach for cardiovascular flow assessment is to reconstruct a complex flow field via CFD simulation, using vascular geometries and input/output functions derived from non-invasive imaging data (Kim et al., 2006; Papathanasopoulou et al., 2003; Steinman et al., 2002). 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. Furthermore, conventional CFD routines typically employ a structured non-Cartesian finite-element mesh that conforms to the vessel geometry, which leads to relatively complex algorithms that incur significant computational costs.

This chapter proposes two new approaches for MRI assessment of carotid flow.

First, we introduce a rapid MRI method for fully-localized measurement of cardiovascular blood flow. The proposed method consists of combining slice-selective spiral imaging with Fourier velocity-encoded (FVE) imaging. The "spiral FVE" method provides intrinsically higher SNR than PC. Furthermore, it is not affected by partial-volume effects, as it measures the velocity distribution within each voxel. This makes this method particularly useful for assessment of flow in stenotic vessels. We show that spiral FVE measurements of carotid flow show good agreement with Doppler ultrasound measurements. We also introduce a mathematical model for deriving spiral FVE velocity distributions from PC velocity maps, and we use this model to show that spiral FVE provides good quantitative agreement with PC measurements, in an experiment with a pulsatile carotid flow phantom. Then, we show that it is possible to accurately estimate *in vivo* wall shear rate (WSR) and oscillatory shear index (OSI) at the carotid bifurcation using spiral FVE.

Finally, we propose a hybrid MRI/CFD approach that integrates low-resolution PC flow measurements directly into the CFD solver. The feasibility of this MRI-driven CFD approach is demonstrated in the carotid bifurcation of a healthy volunteer. We show that MRI-driven CFD has a regularizing effect on the flow fields obtained with MRI alone, and produces flow patterns that are in better agreement with direct MRI measurements than CFD alone. This approach may provide improved estimates of clinically significant blood flow patterns and derived hemodynamic parameters, such as wall shear stress (WSS).

#### **2. Magnetic resonance flow imaging**

#### **2.1 Basic principles of MRI**

MRI is a modality uniquely capable of imaging all aspects of cardiovascular disease, and is a potential "one-stop shop" for cardiovascular health assessment. MRI can generate cross-sectional images in any plane (including oblique planes), and can also measure blood flow. The image acquisition is based on using strong magnetic fields and non-ionizing radiation in the radio frequency range, which are harmless to the patient.

The main component of a MRI scanner is a strong magnetic field, called the *B*<sup>0</sup> field. This magnetic field is always on, even when the scanner is not being used. Typically, MR is used to image hydrogen nuclei, because of its abundance in the human body. Spinning charged particles (or "spins"), such as hydrogen nuclei, act like a tiny bar magnet, presenting a very small magnetic field, emanating from the south pole to the north pole. In normal conditions, each nucleus points to a random direction, resulting in a null net magnetization. However, in the presence of an external magnetic field (such as the *B*<sup>0</sup> field), they will line up with that field. However, they will not all line up in the same direction. Approximately half will point north, and half will point south. Slightly more than half of these spins (about one in 2 Will-be-set-by-IN-TECH

An alternative approach for cardiovascular flow assessment is to reconstruct a complex flow field via CFD simulation, using vascular geometries and input/output functions derived from non-invasive imaging data (Kim et al., 2006; Papathanasopoulou et al., 2003; Steinman et al., 2002). 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. Furthermore, conventional CFD routines typically employ a structured non-Cartesian finite-element mesh that conforms to the vessel geometry, which leads to

First, we introduce a rapid MRI method for fully-localized measurement of cardiovascular blood flow. The proposed method consists of combining slice-selective spiral imaging with Fourier velocity-encoded (FVE) imaging. The "spiral FVE" method provides intrinsically higher SNR than PC. Furthermore, it is not affected by partial-volume effects, as it measures the velocity distribution within each voxel. This makes this method particularly useful for assessment of flow in stenotic vessels. We show that spiral FVE measurements of carotid flow show good agreement with Doppler ultrasound measurements. We also introduce a mathematical model for deriving spiral FVE velocity distributions from PC velocity maps, and we use this model to show that spiral FVE provides good quantitative agreement with PC measurements, in an experiment with a pulsatile carotid flow phantom. Then, we show that it is possible to accurately estimate *in vivo* wall shear rate (WSR) and oscillatory shear

Finally, we propose a hybrid MRI/CFD approach that integrates low-resolution PC flow measurements directly into the CFD solver. The feasibility of this MRI-driven CFD approach is demonstrated in the carotid bifurcation of a healthy volunteer. We show that MRI-driven CFD has a regularizing effect on the flow fields obtained with MRI alone, and produces flow patterns that are in better agreement with direct MRI measurements than CFD alone. This approach may provide improved estimates of clinically significant blood flow patterns and

MRI is a modality uniquely capable of imaging all aspects of cardiovascular disease, and is a potential "one-stop shop" for cardiovascular health assessment. MRI can generate cross-sectional images in any plane (including oblique planes), and can also measure blood flow. The image acquisition is based on using strong magnetic fields and non-ionizing

The main component of a MRI scanner is a strong magnetic field, called the *B*<sup>0</sup> field. This magnetic field is always on, even when the scanner is not being used. Typically, MR is used to image hydrogen nuclei, because of its abundance in the human body. Spinning charged particles (or "spins"), such as hydrogen nuclei, act like a tiny bar magnet, presenting a very small magnetic field, emanating from the south pole to the north pole. In normal conditions, each nucleus points to a random direction, resulting in a null net magnetization. However, in the presence of an external magnetic field (such as the *B*<sup>0</sup> field), they will line up with that field. However, they will not all line up in the same direction. Approximately half will point north, and half will point south. Slightly more than half of these spins (about one in

relatively complex algorithms that incur significant computational costs. This chapter proposes two new approaches for MRI assessment of carotid flow.

index (OSI) at the carotid bifurcation using spiral FVE.

**2. Magnetic resonance flow imaging**

**2.1 Basic principles of MRI**

derived hemodynamic parameters, such as wall shear stress (WSS).

radiation in the radio frequency range, which are harmless to the patient.

a million) will point north, creating a small net magnetization *M*0, which is strong enough to be detected. The net magnetization is proportional to the strength of the *B*<sup>0</sup> field, so MRI scanners with stronger magnetic fields (e.g., 3 Tesla) provide higher SNR.

Other important components of the scanner are the gradient coils. There are typically three gradient coils (*Gx*, *Gy*, and *Gz*), that produce an intentional perturbation in the *B*<sup>0</sup> field when turned on ("played"). This perturbation varies linearly along each spatial direction (*x*, *y* and *z*), such that no perturbation exists at the iso-center of the magnet when these gradients are used. In the presence of an external magnetic field, the spins rotate about the axis of that field. *B*<sup>0</sup> is (approximately) spatially uniform, so all spins initially rotate at the same frequency (the Larmor frequency), *ω* = *γB*0, where *γ* is the gyromagnetic ratio (*γ* = 42.6 MHz/Tesla for hydrogen protons). However, when any of the gradients is played, the magnetic field becomes spatially varying, and so does the rotation frequency of the spins. Therefore, *Gx*, *Gy*, and *Gz* are used to frequency-encode (or phase-encode) spatial position along the *x*, *y* and *z* directions, respectively.

The final major component of the MR scanner is the radio-frequency (RF) coil. This is used to transmit an RF "excitation" pulse to the body, and also to receive the frequency-encoded signal from the "excited" portion of the body. In practice, independent coils may be used for transmission and reception. The RF pulse is typically modulated to the Larmor frequency. While *B*<sup>0</sup> is aligned with the *z*-axis (by definition), *B*1, which is a very weak magnetic field associated with the RF pulse, is aligned with the *x*-axis (also by definition). When the RF pulse is played, some of the spins which are in resonance with the RF pulse (i.e., rotating at the RF pulse's frequency) will now begin to rotate around the *x*-axis (thus the name magnetic resonance). This tilts the net magnetization towards the *x*-*y* plane, and the net magnetization will now have a component in the *x*-*y* plane,*Mxy*.

The RF pulse is typically designed to have a somewhat rectangular profile in Fourier domain, centered at the modulation frequency (e.g., a modulated windowed sinc). This implies that the RF pulse in fact contains a certain range of frequencies, thus all spins rotating within that range become "excited", or tilted towards the *x*-axis. So, by playing gradient(s) of appropriate amplitude, and designing the RF pulse accordingly, one can excite only a thin slice of the body, which correspond to the region containing all spins that are in resonance with the RF pulse's range of frequencies. Excitation profiles other than "slices" may also be obtained (e.g., a pencil beam, or cylindrical excitation (Hu et al., 1993)), by designing an appropriate gradient/pulse combination.

When the RF pulse is turned off, *Mxy* begins to rotate (at the Larmor frequency) around the *z*-axis, as the net magnetization begins to realign with *B*0. This rotating magnetization generates an oscillating signal, which can be detected by the receive coil. The frequency content of the received signal can be used to obtain spatial information about the excited portion of the body. In order to frequency-encode spatial information, gradients are also played during signal acquisition. These are called readout gradients. For imaging a slice perpendicular to the *z*-axis (an axial image), *Gz* is played during excitation (for slice-selection), and *Gx* and *Gy* are played during acquisition. These can be switched, for acquiring sagittal or coronal images, or all three gradients may be used during both excitation and acquisition to image oblique planes.

When the readout gradients are played, the acquired signal at a particular time instant corresponds to the sum of sinusoidal signals generated by spins located at different regions of the body, each rotating at different frequencies corresponding to their spatial locations. If

**2.2 Mathematical formalism**

readout gradients *Gx* and *Gy*:

*m*(*x*, *y*):

gradients:

phase is:

intensity is realized.

As discussed on section 2.1, the acquired MR signal *s*(*t*) at a particular time instant corresponds to a sample of the Fourier transform *M*(*kx*, *ky*) of the excited magnetization

<sup>517</sup> Assessment of Carotid Flow Using Magnetic

*m*(*x*, *y*)*e*

The Fourier coordinates *kx* and *ky* vary with time, according to the zeroth moment of the

 *t* 0

 *t* 0

These equations explain how the gradients can be used to "move" along k-space, as discussed in section 2.1. This formalism can be generalized for any combination of *Gx*, *Gy*, and *Gz*

*m*(*r*) · *e*

 *t* 0

where *G <sup>r</sup>* is the oblique gradient resulting from the combination of the *Gx*, *Gy* and *Gz* gradients, and*r* is its corresponding axis along which the linear variation in magnetic field

Given a spatial position function*r*(*t*) and a magnetic field gradient *G <sup>r</sup>*(*t*), the magnetization

 *t* 0

The basic principles of quantitative flow measurement using magnetic resonance were first proposed by Singer (1959) and Hahn (1960) in the late 1950s. However, clinical applications of MR flow quantification were not reported until the early 1980s (Moran et al., 1985; Nayler et al., 1986; Singer & Crooks, 1983; van Dijk, 1984). Current MR flow imaging methods are based on the fact that spins moving at a constant velocity accrue a phase proportional to the

 *t* 0

2*π*

2*π*

*r*

2*π*

<sup>−</sup>*j*2*π*(*kx <sup>x</sup>*+*ky <sup>y</sup>*) *dx dy*. (1)

*Gx*(*τ*) *dτ* (2)

*Gy*(*τ*) *dτ*. (3)

<sup>−</sup>*j*2*<sup>π</sup>kr*·*<sup>r</sup> <sup>d</sup><sup>r</sup>* (4)

*G <sup>r</sup>*(*τ*) *dτ*, (5)

*<sup>G</sup> <sup>r</sup>*(*τ*) ·*<sup>r</sup>*(*τ*) *<sup>d</sup>τ*, (6)

*G <sup>r</sup>*(*τ*) *dτ* (7)

<sup>=</sup> <sup>2</sup>*<sup>π</sup>kr* ·*<sup>r</sup>*, (8)

*<sup>M</sup>*(*kx*, *ky*) =

Resonance Imaging and Computational Fluid Dynamics

*x y*

*kx*(*t*) = *<sup>γ</sup>*

*ky*(*t*) = *<sup>γ</sup>*

*kr*) =

*kr*(*t*) = *<sup>γ</sup>*

*φ*(*r*, *t*) = *γ*

*φ* = *γr* ·

For static spins,*r*(*t*) is constant (*r*), and this becomes:

as in the exponential in equation 4.

**2.3 Principles of MR flow imaging**

*M*(

an axial slice is being acquired, for example, the demodulated signal value is equivalent to a sample of the Fourier transform *M*(*kx*, *ky*) of the cross-sectional image *m*(*x*, *y*). In this case, by changing the amplitudes of *Gx* and *Gy* during acquisition, one may acquire different samples of *M*(*kx*, *ky*). In fact, by playing *Gx* and/or *Gy*, one can move along the *kx*-*ky* plane (which is known in MRI as "k-space"), collecting samples of *M*(*kx*, *ky*). When enough samples of *M*(*kx*, *ky*) have been collected, an inverse Fourier transform produces *m*(*x*, *y*).

The required coverage of k-space, and the number of samples, depend on the specified spatial resolution and field-of-view. For low spatial resolution imaging, only the central portion of *kx*-*ky* needs to be sampled. For higher spatial resolution, the periphery of k-space must also be covered. The field of view is associated with the spacing between samples. For a larger field-of-view, k-space needs to be more densely sampled, requiring an increased number of samples. If k-space is not sufficiently sampled, and the resulting field-of-view is not large enough to cover the entire object, overlap in spatial domain (aliasing) is observed.

Because signal amplitude is lost as the net magnetization realigns with *B*<sup>0</sup> (this is called relaxation), multiple acquisitions (excitation + readout) may be needed in order to cover the entire k-space. Some trajectories are more efficient in covering k-space than others. For example, spiral imaging, which uses oscillating gradients to achieve spiral k-space trajectories (Figure 1b), are generally faster than 2DFT imaging, i.e., require fewer acquisitions. In 2DFT imaging, each acquisition readout acquires a single line of k-space, sampling *kx*-*ky* in a Cartesian fashion (Figure 1a). This is generally slower, but may be advantageous in some applications with respect to the nature of associated image artifacts. The full sequence of RF pulses and gradients is called a "pulse sequence". The time between acquisitions is called the pulse repetition time, or TR.

Fig. 1. Timing diagram (left) and corresponding k-space trajectories (right) for (a) 2DFT, and (b) spiral acquisitions.

#### **2.2 Mathematical formalism**

4 Will-be-set-by-IN-TECH

an axial slice is being acquired, for example, the demodulated signal value is equivalent to a sample of the Fourier transform *M*(*kx*, *ky*) of the cross-sectional image *m*(*x*, *y*). In this case, by changing the amplitudes of *Gx* and *Gy* during acquisition, one may acquire different samples of *M*(*kx*, *ky*). In fact, by playing *Gx* and/or *Gy*, one can move along the *kx*-*ky* plane (which is known in MRI as "k-space"), collecting samples of *M*(*kx*, *ky*). When enough samples of

The required coverage of k-space, and the number of samples, depend on the specified spatial resolution and field-of-view. For low spatial resolution imaging, only the central portion of *kx*-*ky* needs to be sampled. For higher spatial resolution, the periphery of k-space must also be covered. The field of view is associated with the spacing between samples. For a larger field-of-view, k-space needs to be more densely sampled, requiring an increased number of samples. If k-space is not sufficiently sampled, and the resulting field-of-view is not large

Because signal amplitude is lost as the net magnetization realigns with *B*<sup>0</sup> (this is called relaxation), multiple acquisitions (excitation + readout) may be needed in order to cover the entire k-space. Some trajectories are more efficient in covering k-space than others. For example, spiral imaging, which uses oscillating gradients to achieve spiral k-space trajectories (Figure 1b), are generally faster than 2DFT imaging, i.e., require fewer acquisitions. In 2DFT imaging, each acquisition readout acquires a single line of k-space, sampling *kx*-*ky* in a Cartesian fashion (Figure 1a). This is generally slower, but may be advantageous in some applications with respect to the nature of associated image artifacts. The full sequence of RF pulses and gradients is called a "pulse sequence". The time between acquisitions is called the

Fig. 1. Timing diagram (left) and corresponding k-space trajectories (right) for (a) 2DFT,

kx

kx

ky

ky

*M*(*kx*, *ky*) have been collected, an inverse Fourier transform produces *m*(*x*, *y*).

enough to cover the entire object, overlap in spatial domain (aliasing) is observed.

pulse repetition time, or TR. **a**

**b**

and (b) spiral acquisitions.

Gz

RF

Gz

RF

Gx

Gy

Gx

Gy

As discussed on section 2.1, the acquired MR signal *s*(*t*) at a particular time instant corresponds to a sample of the Fourier transform *M*(*kx*, *ky*) of the excited magnetization *m*(*x*, *y*):

$$M(k\_{\mathbf{x}}, k\_{\mathbf{y}}) = \int\_{\mathbf{x}} \int\_{\mathbf{y}} m(\mathbf{x}, \mathbf{y}) \, e^{-j2\pi(k\_{\mathbf{x}}\mathbf{x} + k\_{\mathbf{y}}\mathbf{y})} \, d\mathbf{x} \, dy. \tag{1}$$

The Fourier coordinates *kx* and *ky* vary with time, according to the zeroth moment of the readout gradients *Gx* and *Gy*:

$$k\_{\mathbf{x}}(t) = \frac{\gamma}{2\pi} \int\_{0}^{t} \mathbf{G}\_{\mathbf{x}}(\tau) \, d\tau \tag{2}$$

$$k\_{\mathcal{Y}}(t) = \frac{\gamma}{2\pi} \int\_0^t \mathcal{G}\_{\mathcal{Y}}(\tau) \, d\tau. \tag{3}$$

These equations explain how the gradients can be used to "move" along k-space, as discussed in section 2.1. This formalism can be generalized for any combination of *Gx*, *Gy*, and *Gz* gradients:

$$M(\vec{k}\_r) = \int\_{\vec{r}} m(\vec{r}) \cdot e^{-j2\pi \vec{k}\_r \cdot \vec{r}} \, d\vec{r} \tag{4}$$

$$
\vec{k}\_r(t) = \frac{\gamma}{2\pi} \int\_0^t \vec{G}\_r(\tau) \, d\tau,\tag{5}
$$

where *G <sup>r</sup>* is the oblique gradient resulting from the combination of the *Gx*, *Gy* and *Gz* gradients, and*r* is its corresponding axis along which the linear variation in magnetic field intensity is realized.

Given a spatial position function*r*(*t*) and a magnetic field gradient *G <sup>r</sup>*(*t*), the magnetization phase is:

$$\Phi(\vec{r},t) = \gamma \int\_0^t \vec{G}\_r(\tau) \cdot \vec{r}(\tau) \,d\tau,\tag{6}$$

For static spins,*r*(*t*) is constant (*r*), and this becomes:

$$
\phi = \gamma \vec{r} \cdot \int\_0^t \vec{G}\_r(\tau) \, d\tau \tag{7}
$$

$$
\vec{k} = 2\pi \vec{k}\_r \cdot \vec{r}\_r \tag{8}
$$

as in the exponential in equation 4.

#### **2.3 Principles of MR flow imaging**

The basic principles of quantitative flow measurement using magnetic resonance were first proposed by Singer (1959) and Hahn (1960) in the late 1950s. However, clinical applications of MR flow quantification were not reported until the early 1980s (Moran et al., 1985; Nayler et al., 1986; Singer & Crooks, 1983; van Dijk, 1984). Current MR flow imaging methods are based on the fact that spins moving at a constant velocity accrue a phase proportional to the

*kv* <sup>=</sup> *<sup>γ</sup>* 2*π*

four-dimensional dataset: *m*(*x*, *y*, *v*, *t*).

Resonance Imaging and Computational Fluid Dynamics

where *δ*(*v*) is the Dirac delta function.

*<sup>s</sup>*ˆ(*x*, *<sup>y</sup>*, *<sup>v</sup>*) = [*m*(*x*, *<sup>y</sup>*) <sup>×</sup> *<sup>δ</sup>*( *<sup>v</sup>* <sup>−</sup> *vo*(*x*, *<sup>y</sup>*) )] <sup>∗</sup> sinc *<sup>x</sup>*

*<sup>m</sup>*(*x*, *<sup>y</sup>*) <sup>×</sup> sinc *<sup>v</sup>* <sup>−</sup> *vo*(*x*, *<sup>y</sup>*)

**2.6 FVE signal model**

resolution, as follows:

*<sup>s</sup>*ˆ(*x*, *<sup>y</sup>*, *<sup>v</sup>*) =

simulation purposes (Carvalho et al., 2010).

equivalent to:

approximately:

Each voxel of the two-dimensional image is associated with a distribution of velocities. This three-dimensional function *m*(*x*, *y*, *v*) is associated with a three-dimensional Fourier space *M*(*kx*, *ky*, *kv*). Thus, an extra dimension is added to k-space, and multiple acquisitions are required to cover the entire *kx*-*ky*-*kv* space. In order to move along *kv*, a bipolar gradient with the appropriate amplitude (and first moment) is played before the *kx*-*ky* readout gradients, in each acquisition. Placing the bipolar gradient along the *z*-axis will encode through-plane velocities. Placing the bipolar gradient along *x* or *y* will encode in-plane velocities. Oblique flow can be encoded using a combination of bipolar gradients along the *x*, *y* and *z* axes. Each acquisition along *kv* is called a velocity encode. The number of required velocity encodes depends on the desired velocity resolution and velocity field-of-view (the maximum range of velocities measured without aliasing). For example, to obtain a 25 cm/s resolution over a 600 cm/s field-of-view, 24 velocity encodes are needed. The spatial-velocity distribution, *m*(*x*, *y*, *v*), is obtained by inverse Fourier transforming the acquired data, *M*(*kx*, *ky*, *kv*). If cine imaging (Glover & Pelc, 1988) is used, measurements are also time resolved, resulting in a

<sup>519</sup> Assessment of Carotid Flow Using Magnetic

2DFT phase contrast provides two two-dimensional functions, *m*(*x*, *y*) and *vo*(*x*, *y*), the magnitude and velocity maps, respectively. If these maps are measured with sufficiently high spatial resolution, and flow is laminar, one can assume that each voxel contains only one velocity, and therefore the spatial-velocity distribution associated with the object is

In 2DFT FVE, k-space data is truncated to a rectangular cuboid in *kx*-*ky*-*kv* space. The associated object domain spatial-velocity blurring can be modeled as a convolution of the true object distribution, *s*(*x*, *y*, *v*), with sinc(*x*/Δ*x*), sinc(*y*/Δ*y*), and sinc(*v*/Δ*v*), where Δ*x* and Δ*y* are the spatial resolutions along the *x* and *y* axes, respectively, and Δ*v* is the velocity

where *s*ˆ(*x*, *y*, *v*) is the measured object distribution, and ∗ denotes convolution. This is

Δ*v*

This approach for deriving FVE data from high-resolution velocity maps can be used for many

Δ*x* 

 ∗ 

*s*(*x*, *y*, *v*) = *m*(*x*, *y*) × *δ*( *v* − *vo*(*x*, *y*) ), (15)

<sup>∗</sup> sinc *<sup>y</sup>*

sinc *<sup>x</sup>* Δ*x* 

Δ*y* 

<sup>∗</sup> sinc *<sup>v</sup>* Δ*v* 

<sup>×</sup> sinc *<sup>y</sup>*

Δ*y*

, (16)

. (17)

*M* 1. (14)

velocity times the first moment of the gradient waveform along the direction in which they are moving.

For spins moving along the *r*-axis with a constant velocity *v*, and initial position*r*0, we can write*r*(*t*) = *r*<sup>0</sup> +*vt*. Rewriting equation 6, for *t* = *t*0:

$$
\phi = \gamma \int\_0^{t\_0} \vec{G}\_r(t) \cdot (\vec{r}\_0 + \vec{v}t) \, dt \tag{9}
$$

$$\vec{r} = \gamma \,\vec{r}\_0 \cdot \int\_0^{t\_0} \vec{G}\_r(t) \, dt + \gamma \,\vec{v} \cdot \int\_0^{t\_0} \vec{G}\_r(t) \, t \, dt \tag{10}$$

$$
\vec{r} = \gamma \vec{r}\_0 \cdot \vec{M}\_0 + \gamma \vec{v} \cdot \vec{M}\_{1\prime} \tag{11}
$$

where *M* <sup>0</sup> and *M* <sup>1</sup> are the zeroth and first moments of the*r*-gradient waveform at the time of signal acquisitions ("echo time", or "time to echo" (TE)), respectively. Thus, if a gradient with null zeroth moment is used (e.g., a bipolar gradient, aligned with *v*), the phase accrued for a constant velocity spin is *<sup>φ</sup>* <sup>=</sup> *<sup>γ</sup><sup>v</sup>* · *<sup>M</sup>* 1.

Therefore, if a bipolar gradient waveform is played between the excitation and the readout, the phase measured in a pixel of the acquired image is directly proportional to the velocity of the spins contained within its corresponding voxel. However, factors other than flow (such as inhomogeneities of the magnetic field) may cause additional phase shifts that would cause erroneous interpretation of the local velocity (Rebergen et al., 1993).

#### **2.4 Phase contrast**

The phase contrast method addresses the problem mentioned above by using two gradient-echo data acquisitions in which the first moment of the bipolar gradient waveform is varied between measurements (O'Donnell, 1985). The velocity in each voxel is measured as:

$$v(\mathbf{x}, y) = \frac{\phi\_a(\mathbf{x}, y) - \phi\_b(\mathbf{x}, y)}{\gamma (M\_1^a - M\_1^b)},\tag{12}$$

where *φa*(*x*, *y*) and *φb*(*x*, *y*) are the phase images acquired in each acquisition, and *M<sup>a</sup>* <sup>1</sup> and *<sup>M</sup><sup>b</sup>* 1 are the first moment of the bipolar gradients used in each acquisition.

#### **2.5 Fourier velocity encoding**

While phase contrast provides a single velocity measurement associated with each voxel, Fourier velocity encoding (Moran, 1982) provides a velocity histogram for each spatial location, which is a measurement of the velocity distribution within each voxel.

FVE involves phase-encoding along a velocity dimension. Instead of only two acquisitions, as in phase contrast, multiple acquisitions are performed, and a bipolar gradient with a different amplitude (and first moment) is used in each acquisition. Equation 10 can be rewritten as:

$$\phi(\vec{r}, \vec{v}, t) = 2\pi(\vec{k}\_r \cdot \vec{r} + \vec{k}\_\upsilon \cdot \vec{v}),\tag{13}$$

where*kv* is the velocity frequency variable associated with *v*, and is proportional to the first moment of *G <sup>r</sup>*(*t*):

$$
\vec{k}\_v = \frac{\gamma}{2\pi} \vec{M}\_1. \tag{14}
$$

Each voxel of the two-dimensional image is associated with a distribution of velocities. This three-dimensional function *m*(*x*, *y*, *v*) is associated with a three-dimensional Fourier space *M*(*kx*, *ky*, *kv*). Thus, an extra dimension is added to k-space, and multiple acquisitions are required to cover the entire *kx*-*ky*-*kv* space. In order to move along *kv*, a bipolar gradient with the appropriate amplitude (and first moment) is played before the *kx*-*ky* readout gradients, in each acquisition. Placing the bipolar gradient along the *z*-axis will encode through-plane velocities. Placing the bipolar gradient along *x* or *y* will encode in-plane velocities. Oblique flow can be encoded using a combination of bipolar gradients along the *x*, *y* and *z* axes.

Each acquisition along *kv* is called a velocity encode. The number of required velocity encodes depends on the desired velocity resolution and velocity field-of-view (the maximum range of velocities measured without aliasing). For example, to obtain a 25 cm/s resolution over a 600 cm/s field-of-view, 24 velocity encodes are needed. The spatial-velocity distribution, *m*(*x*, *y*, *v*), is obtained by inverse Fourier transforming the acquired data, *M*(*kx*, *ky*, *kv*). If cine imaging (Glover & Pelc, 1988) is used, measurements are also time resolved, resulting in a four-dimensional dataset: *m*(*x*, *y*, *v*, *t*).

#### **2.6 FVE signal model**

6 Will-be-set-by-IN-TECH

velocity times the first moment of the gradient waveform along the direction in which they

For spins moving along the *r*-axis with a constant velocity *v*, and initial position*r*0, we can

*<sup>G</sup> <sup>r</sup>*(*t*) *dt* <sup>+</sup> *<sup>γ</sup><sup>v</sup>* ·

where *M* <sup>0</sup> and *M* <sup>1</sup> are the zeroth and first moments of the*r*-gradient waveform at the time of signal acquisitions ("echo time", or "time to echo" (TE)), respectively. Thus, if a gradient with null zeroth moment is used (e.g., a bipolar gradient, aligned with *v*), the phase accrued for a

Therefore, if a bipolar gradient waveform is played between the excitation and the readout, the phase measured in a pixel of the acquired image is directly proportional to the velocity of the spins contained within its corresponding voxel. However, factors other than flow (such as inhomogeneities of the magnetic field) may cause additional phase shifts that would cause

The phase contrast method addresses the problem mentioned above by using two gradient-echo data acquisitions in which the first moment of the bipolar gradient waveform is varied between measurements (O'Donnell, 1985). The velocity in each voxel is measured as:

> *<sup>v</sup>*(*x*, *<sup>y</sup>*) = *<sup>φ</sup>a*(*x*, *<sup>y</sup>*) <sup>−</sup> *<sup>φ</sup>b*(*x*, *<sup>y</sup>*) *γ*(*M<sup>a</sup>*

While phase contrast provides a single velocity measurement associated with each voxel, Fourier velocity encoding (Moran, 1982) provides a velocity histogram for each spatial

FVE involves phase-encoding along a velocity dimension. Instead of only two acquisitions, as in phase contrast, multiple acquisitions are performed, and a bipolar gradient with a different amplitude (and first moment) is used in each acquisition. Equation 10 can be rewritten as:

where*kv* is the velocity frequency variable associated with *v*, and is proportional to the first

where *φa*(*x*, *y*) and *φb*(*x*, *y*) are the phase images acquired in each acquisition, and *M<sup>a</sup>*

location, which is a measurement of the velocity distribution within each voxel.

*φ*(*r*,*v*, *t*) = 2*π*(

<sup>1</sup> <sup>−</sup> *<sup>M</sup><sup>b</sup>*

*<sup>G</sup> <sup>r</sup>*(*t*) · (*<sup>r</sup>*<sup>0</sup> <sup>+</sup>*vt*) *dt* (9)

*G <sup>r</sup>*(*t*)*t dt* (10)

<sup>1</sup>) , (12)

*kr* ·*<sup>r</sup>* <sup>+</sup>*kv* ·*<sup>v</sup>*), (13)

<sup>1</sup> and *<sup>M</sup><sup>b</sup>* 1

 *t*<sup>0</sup> 0

<sup>=</sup> *<sup>γ</sup><sup>r</sup>*<sup>0</sup> · *<sup>M</sup>* <sup>0</sup> <sup>+</sup> *<sup>γ</sup><sup>v</sup>* · *<sup>M</sup>* 1, (11)

write*r*(*t*) = *r*<sup>0</sup> +*vt*. Rewriting equation 6, for *t* = *t*0:

constant velocity spin is *<sup>φ</sup>* <sup>=</sup> *<sup>γ</sup><sup>v</sup>* · *<sup>M</sup>* 1.

**2.4 Phase contrast**

**2.5 Fourier velocity encoding**

moment of *G <sup>r</sup>*(*t*):

*φ* = *γ*

= *γr*<sup>0</sup> ·

 *t*<sup>0</sup> 0

erroneous interpretation of the local velocity (Rebergen et al., 1993).

are the first moment of the bipolar gradients used in each acquisition.

 *t*<sup>0</sup> 0

are moving.

2DFT phase contrast provides two two-dimensional functions, *m*(*x*, *y*) and *vo*(*x*, *y*), the magnitude and velocity maps, respectively. If these maps are measured with sufficiently high spatial resolution, and flow is laminar, one can assume that each voxel contains only one velocity, and therefore the spatial-velocity distribution associated with the object is approximately:

$$\mathbf{s}(\mathbf{x}, \mathbf{y}, \boldsymbol{\upsilon}) = m(\mathbf{x}, \mathbf{y}) \times \delta(\boldsymbol{\upsilon} - \boldsymbol{\upsilon}\_o(\mathbf{x}, \mathbf{y})\ ),\tag{15}$$

where *δ*(*v*) is the Dirac delta function.

In 2DFT FVE, k-space data is truncated to a rectangular cuboid in *kx*-*ky*-*kv* space. The associated object domain spatial-velocity blurring can be modeled as a convolution of the true object distribution, *s*(*x*, *y*, *v*), with sinc(*x*/Δ*x*), sinc(*y*/Δ*y*), and sinc(*v*/Δ*v*), where Δ*x* and Δ*y* are the spatial resolutions along the *x* and *y* axes, respectively, and Δ*v* is the velocity resolution, as follows:

$$\mathcal{S}(\mathbf{x}, y, \upsilon) = \left[ m(\mathbf{x}, y) \times \delta(\upsilon - \upsilon\_o(\mathbf{x}, y) \mid) \right] \ast \text{sinc} \left( \frac{\mathbf{x}}{\Delta \mathbf{x}} \right) \ast \text{sinc} \left( \frac{y}{\Delta y} \right) \ast \text{sinc} \left( \frac{\upsilon}{\Delta \upsilon} \right), \tag{16}$$

where *s*ˆ(*x*, *y*, *v*) is the measured object distribution, and ∗ denotes convolution. This is equivalent to:

$$\mathbf{s}(\mathbf{x}, y, \upsilon) = \left[ m(\mathbf{x}, y) \times \text{sinc}\left(\frac{\upsilon - \upsilon\_o(\mathbf{x}, y)}{\Delta \upsilon}\right) \right] \ast \left[ \text{sinc}\left(\frac{\mathbf{x}}{\Delta x}\right) \times \text{sinc}\left(\frac{y}{\Delta y}\right) \right]. \tag{17}$$

This approach for deriving FVE data from high-resolution velocity maps can be used for many simulation purposes (Carvalho et al., 2010).

**3.2 Spiral FVE signal model**

resulting in:

**3.3** *In vitro* **validation**

As spiral FVE acquisitions follow a stack-of-spirals pattern in *kx*-*ky*-*kv* space (Figure 3), k-space data is truncated to a cylinder, i.e., a circle along *kx*-*ky* (with diameter 1/Δ*r*), and a rect function along *kv* (with width 1/Δ*v*), where Δ*r* and Δ*v* are the prescribed spatial and velocity resolutions, respectively. Using the same approach we used in section 2.6 for 2DFT FVE, the associated object domain spatial-velocity blurring in spiral FVE can be modeled as a

<sup>521</sup> Assessment of Carotid Flow Using Magnetic

Δ*v*

An *in vitro* comparison of velocity distributions measured with spiral FVE with those derived from high-resolution 2DFT phase contrast — the current MR gold standard — was performed. The signal model presented in section 3.2 was used to generate simulated FVE data based on

The validation experiments were performed using a pulsatile carotid flow phantom (Phantoms by Design, Inc., Bothell, WA). A slice perpendicular to the phantom's carotid bifurcation was prescribed, and through-plane velocities were measured. A cine gradient-echo 2DFT phase contrast sequence with high spatial resolution and high SNR (0.33 mm resolution, 10 averages, 80 cm/s Venc) was used as a reference. Cine spiral FVE data with Δ*r* = 3 mm and Δ*v* = 10 cm/s was obtained from the same scan plane. Both acquisitions were prospectively gated, and used the same TR (11.6 ms), flip angle (30◦), slice profile (3 mm), temporal resolution (23.2 ms), and pre-scan settings. The total scan time was

A simulated spiral FVE dataset was computed from the PC magnitude and velocity maps, using the convolution model described in Eq. 18. The PC-derived and FVE-measured data were registered by taking one magnitude image *m*(*x*, *y*) from each dataset, and then using the phase difference between their Fourier transforms *M*(*kx*, *ky*) to estimate the spatial shift between the images. Amplitude scaling was performed by normalizing the 2-norm of each FVE dataset. The difference between PC-derived and FVE-measured time-velocity distributions was calculated for select voxels, and the associated signal-to-error ratios were

Figure 4 shows measured and PC-derived time-velocity FVE distributions from nine representative voxels, selected around the circumference of the vessel wall of the pulsatile carotid flow phantom's bifurcation. The signal-to-error ratio between measured and PC-derived time-velocity distributions was measured to be within 9.3–11.7 dB. Imperfect registration between the datasets, combined with spatial blurring due to off-resonance in the measured spiral FVE data, may have contributed to this moderate signal-to-error ratio. Nevertheless, the two datasets show good visual agreement, and no significant spatial variation was observed in terms of accuracy. These results show that velocity distributions

computed. This was used as a quantitative assessment of spiral FVE's accuracy.

Δ*v* 

<sup>∗</sup> jinc *<sup>x</sup>*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*<sup>2</sup> Δ*r*

*x*<sup>2</sup> + *y*2/Δ*r*) and sinc(*v*/Δ*v*),

. (18)

<sup>∗</sup> jinc *<sup>x</sup>*<sup>2</sup> <sup>+</sup> *<sup>y</sup>*<sup>2</sup> Δ*r*

convolution of the true object distribution, *s*(*x*, *y*, *v*), with jinc(

Resonance Imaging and Computational Fluid Dynamics

= 

high-resolution 2DFT phase contrast data.

40 minutes for phase contrast, and 12 seconds for FVE.

*<sup>s</sup>*ˆ(*x*, *<sup>y</sup>*, *<sup>v</sup>*) = [*m*(*x*, *<sup>y</sup>*) <sup>×</sup> *<sup>δ</sup>*( *<sup>v</sup>* <sup>−</sup> *vo*(*x*, *<sup>y</sup>*) )] <sup>∗</sup> sinc *<sup>v</sup>*

*<sup>m</sup>*(*x*, *<sup>y</sup>*) <sup>×</sup> sinc *<sup>v</sup>* <sup>−</sup> *vo*(*x*, *<sup>y</sup>*)

#### **3. Slice-selective FVE with spiral readouts (spiral FVE)**

Phase contrast imaging is fast, but has limitations associated with partial-volume effects. On the other hand, 2DFT FVE addresses those limitations, but requires long scan times. Thus, we propose the use of slice-selective FVE MRI with spiral acquisitions. The proposed spiral FVE method is capable of acquiring fully localized, time-resolved velocity distributions in a short breath-hold. Scan-plane prescription is performed using classic protocols.

We present practical implementations for measuring blood flow through the carotid arteries, and comparisons with Doppler ultrasound and high-resolution 2DFT phase contrast MRI. The proposed method is demonstrated in healthy volunteers. Subjects provided informed consent, and were imaged using a protocol approved by the institutional review board of the University of Southern California.

#### **3.1 Pulse sequence**

The spiral FVE imaging pulse sequence (Figure 2) consists of a slice-selective excitation, a velocity-encoding bipolar gradient, a spiral readout, and refocusing and spoiling gradients. The dataset corresponding to each temporal frame is a stack-of-spirals in *kx*-*ky*-*kv* space (Figure 3). The bipolar gradient effectively phase-encodes in *kv*, while each spiral readout acquires one "disc" in *kx*-*ky*.

Fig. 2. Spiral FVE pulse sequence. It consists of (a) slice selective excitation, (b) velocity encoding bipolar gradient, (c) spiral readout, and (d) refocusing and spoiling gradients.

Fig. 3. Spiral FVE k-space sampling scheme. The dataset corresponding to each temporal frame is a stack-of-spirals in *kx*-*ky*-*kv* space. Each spiral acquisition corresponds to a different *kv* encode level.

#### **3.2 Spiral FVE signal model**

8 Will-be-set-by-IN-TECH

Phase contrast imaging is fast, but has limitations associated with partial-volume effects. On the other hand, 2DFT FVE addresses those limitations, but requires long scan times. Thus, we propose the use of slice-selective FVE MRI with spiral acquisitions. The proposed spiral FVE method is capable of acquiring fully localized, time-resolved velocity distributions in a short

We present practical implementations for measuring blood flow through the carotid arteries, and comparisons with Doppler ultrasound and high-resolution 2DFT phase contrast MRI. The proposed method is demonstrated in healthy volunteers. Subjects provided informed consent, and were imaged using a protocol approved by the institutional review board of the

The spiral FVE imaging pulse sequence (Figure 2) consists of a slice-selective excitation, a velocity-encoding bipolar gradient, a spiral readout, and refocusing and spoiling gradients. The dataset corresponding to each temporal frame is a stack-of-spirals in *kx*-*ky*-*kv* space (Figure 3). The bipolar gradient effectively phase-encodes in *kv*, while each spiral readout

**ab c d**

Fig. 2. Spiral FVE pulse sequence. It consists of (a) slice selective excitation, (b) velocity encoding bipolar gradient, (c) spiral readout, and (d) refocusing and spoiling gradients.

kx

24:

1:

kv

Fig. 3. Spiral FVE k-space sampling scheme. The dataset corresponding to each temporal frame is a stack-of-spirals in *kx*-*ky*-*kv* space. Each spiral acquisition corresponds to a different

ky

**3. Slice-selective FVE with spiral readouts (spiral FVE)**

University of Southern California.

acquires one "disc" in *kx*-*ky*.

**RF**

**Gz**

**Gx Gy**

**3.1 Pulse sequence**

*kv* encode level.

breath-hold. Scan-plane prescription is performed using classic protocols.

As spiral FVE acquisitions follow a stack-of-spirals pattern in *kx*-*ky*-*kv* space (Figure 3), k-space data is truncated to a cylinder, i.e., a circle along *kx*-*ky* (with diameter 1/Δ*r*), and a rect function along *kv* (with width 1/Δ*v*), where Δ*r* and Δ*v* are the prescribed spatial and velocity resolutions, respectively. Using the same approach we used in section 2.6 for 2DFT FVE, the associated object domain spatial-velocity blurring in spiral FVE can be modeled as a convolution of the true object distribution, *s*(*x*, *y*, *v*), with jinc( *x*<sup>2</sup> + *y*2/Δ*r*) and sinc(*v*/Δ*v*), resulting in:

$$\mathcal{S}(\mathbf{x}, \boldsymbol{y}, \upsilon) = \left[ m(\mathbf{x}, \boldsymbol{y}) \times \delta(\upsilon - v\_o(\mathbf{x}, \boldsymbol{y}) \,) \right] \ast \text{sinc} \left( \frac{\upsilon}{\Delta \upsilon} \right) \ast \text{jinc} \left( \frac{\sqrt{\mathbf{x}^2 + \mathbf{y}^2}}{\Delta r} \right)$$

$$= \left[ m(\mathbf{x}, \boldsymbol{y}) \times \text{sinc} \left( \frac{\upsilon - v\_o(\mathbf{x}, \boldsymbol{y})}{\Delta \upsilon} \right) \right] \ast \text{jinc} \left( \frac{\sqrt{\mathbf{x}^2 + \mathbf{y}^2}}{\Delta r} \right) . \tag{18}$$

#### **3.3** *In vitro* **validation**

An *in vitro* comparison of velocity distributions measured with spiral FVE with those derived from high-resolution 2DFT phase contrast — the current MR gold standard — was performed. The signal model presented in section 3.2 was used to generate simulated FVE data based on high-resolution 2DFT phase contrast data.

The validation experiments were performed using a pulsatile carotid flow phantom (Phantoms by Design, Inc., Bothell, WA). A slice perpendicular to the phantom's carotid bifurcation was prescribed, and through-plane velocities were measured. A cine gradient-echo 2DFT phase contrast sequence with high spatial resolution and high SNR (0.33 mm resolution, 10 averages, 80 cm/s Venc) was used as a reference. Cine spiral FVE data with Δ*r* = 3 mm and Δ*v* = 10 cm/s was obtained from the same scan plane. Both acquisitions were prospectively gated, and used the same TR (11.6 ms), flip angle (30◦), slice profile (3 mm), temporal resolution (23.2 ms), and pre-scan settings. The total scan time was 40 minutes for phase contrast, and 12 seconds for FVE.

A simulated spiral FVE dataset was computed from the PC magnitude and velocity maps, using the convolution model described in Eq. 18. The PC-derived and FVE-measured data were registered by taking one magnitude image *m*(*x*, *y*) from each dataset, and then using the phase difference between their Fourier transforms *M*(*kx*, *ky*) to estimate the spatial shift between the images. Amplitude scaling was performed by normalizing the 2-norm of each FVE dataset. The difference between PC-derived and FVE-measured time-velocity distributions was calculated for select voxels, and the associated signal-to-error ratios were computed. This was used as a quantitative assessment of spiral FVE's accuracy.

Figure 4 shows measured and PC-derived time-velocity FVE distributions from nine representative voxels, selected around the circumference of the vessel wall of the pulsatile carotid flow phantom's bifurcation. The signal-to-error ratio between measured and PC-derived time-velocity distributions was measured to be within 9.3–11.7 dB. Imperfect registration between the datasets, combined with spatial blurring due to off-resonance in the measured spiral FVE data, may have contributed to this moderate signal-to-error ratio. Nevertheless, the two datasets show good visual agreement, and no significant spatial variation was observed in terms of accuracy. These results show that velocity distributions

other. Therefore, the unaliased field-of-view is wide enough to enclose all vessels on one side of the neck, so that quantification of a vessel of interest will not be disturbed by flow from neighboring vessels (e.g. measurement of flow in the left carotid arteries will not be disturbed by flow in the left jugular vein). In order to suppress signal from the opposite side of the neck, we separately reconstruct data from the left and right neck-coil elements. This uses the receiver coil sensitivity profile to avoid spiral view-sharing artifacts. Based on these

<sup>523</sup> Assessment of Carotid Flow Using Magnetic

time (ms) 0 200 400 600

observations, we used spiral interleaf view-sharing in all subsequent carotid studies.

Fig. 5. Comparison of different view-orderings for multi-shot spiral FVE, in a healthy volunteer carotid study. When two or more *kv* levels are acquired during the same heartbeat (a), velocity distribution changes between consecutive TRs cause ghosting artifacts along the velocity axis (arrow). This artifact is not seen if, in the same heartbeat, different spiral

Representative *in vivo* results are compared with Doppler ultrasound in Figure 6. The MRI measured time-velocity histogram shows good agreement with the ultrasound measurement, as the peak velocity and the shape of the flow waveform were comparable to those observed

c

Fig. 6. Comparison of the spiral FVE method (a) with Doppler ultrasound (b), in a healthy volunteer common carotid artery study. Peak velocity and time-velocity waveforms are in

**b**

time (ms) <sup>0</sup> <sup>200</sup> <sup>400</sup> <sup>600</sup> <sup>−</sup><sup>100</sup>

time (ms)

0 200 400 600 800

interleaves, but only one *kv* encoding, are acquired (b).

**3.4.2 Comparison with Doppler ultrasound**

**a b**

Resonance Imaging and Computational Fluid Dynamics

−50

in the ultrasound studies.

**a**

velocity (cm/s)

good agreement.

−20

60

140

0

velocity (cm/s)

50

100

measured with spiral FVE agree well with those obtained with 2DFT phase contrast, the current MRI gold standard.

Fig. 4. *In vitro* evaluation of the accuracy of spiral FVE velocity histograms. Results are shown for nine representative voxels, selected around the circumference of the vessel wall of the pulsatile carotid flow phantom's bifurcation (a). For each voxel, it is shown: (b) time-velocity distribution derived from high-resolution 2DFT phase contrast; (c) time-velocity distribution measured with spiral FVE; (d) absolute difference between spiral FVE and 2DFT PC-derived histograms; (e) signal-to-error ratio.

#### **3.4** *In vivo* **evaluation**

The spiral FVE method was evaluated *in vivo*, aiming at quantifying flow through the common carotid artery of a healthy volunteer. Doppler ultrasound was used as a gold standard and was qualitatively compared with the proposed method. We also show experiments for evaluating the most appropriate view-ordering scheme for measuring carotid flow, and we demonstrate spiral FVE's ability to measure multiple flows in a single acquisition.

#### **3.4.1 View ordering**

Experiments were performed to determine the appropriate view-ordering scheme for carotid spiral FVE studies. Flow was measured through the carotid artery of a healthy volunteer, using a 4-interleaf spiral FVE acquisition. The measurement was performed twice, and in each acquisition, a different view-ordering scheme was used. In the first acquisition, one spiral interleaf was acquired per heartbeat, and two different *kv* levels were encoded throughout each R-R interval. In the second acquisition, two spiral interleaves were acquired per heartbeat, encoding one *kv* level per cardiac cycle. Reconstructed velocity histograms were compared with respect to data inconsistency artifacts related to view-ordering.

The results are shown in Figure 5. Note that the ghosting artifacts that arise in the velocity histograms when acquiring two different velocity encode levels during the same heartbeat (Figure 5a) do not appear when we used interleaf segmentation instead (Figure 5b). In contrast to view-sharing along *kv*, which causes ghosting in the velocity direction, view-sharing among spiral interleaves introduces swirling artifacts in image domain, reducing the effective unaliased spatial field-of-view by a factor of 2. However, only moving spins (flowing blood) experience these artifacts, and vessels on the same side of the neck are relatively close to each 10 Will-be-set-by-IN-TECH

measured with spiral FVE agree well with those obtained with 2DFT phase contrast, the

4

5

6

7

8

9

3

current MRI gold standard.

2 3 4

9 8 7

**3.4** *In vivo* **evaluation**

**3.4.1 View ordering**

6 5

1

**a**

**b** velocity 1

**c**

**d**

**e**

10.7 dB

FVE and 2DFT PC-derived histograms; (e) signal-to-error ratio.

spiral FVE's ability to measure multiple flows in a single acquisition.

compared with respect to data inconsistency artifacts related to view-ordering.

11.3 dB

9.3 dB

Fig. 4. *In vitro* evaluation of the accuracy of spiral FVE velocity histograms. Results are shown for nine representative voxels, selected around the circumference of the vessel wall of

time-velocity distribution measured with spiral FVE; (d) absolute difference between spiral

The spiral FVE method was evaluated *in vivo*, aiming at quantifying flow through the common carotid artery of a healthy volunteer. Doppler ultrasound was used as a gold standard and was qualitatively compared with the proposed method. We also show experiments for evaluating the most appropriate view-ordering scheme for measuring carotid flow, and we demonstrate

Experiments were performed to determine the appropriate view-ordering scheme for carotid spiral FVE studies. Flow was measured through the carotid artery of a healthy volunteer, using a 4-interleaf spiral FVE acquisition. The measurement was performed twice, and in each acquisition, a different view-ordering scheme was used. In the first acquisition, one spiral interleaf was acquired per heartbeat, and two different *kv* levels were encoded throughout each R-R interval. In the second acquisition, two spiral interleaves were acquired per heartbeat, encoding one *kv* level per cardiac cycle. Reconstructed velocity histograms were

The results are shown in Figure 5. Note that the ghosting artifacts that arise in the velocity histograms when acquiring two different velocity encode levels during the same heartbeat (Figure 5a) do not appear when we used interleaf segmentation instead (Figure 5b). In contrast to view-sharing along *kv*, which causes ghosting in the velocity direction, view-sharing among spiral interleaves introduces swirling artifacts in image domain, reducing the effective unaliased spatial field-of-view by a factor of 2. However, only moving spins (flowing blood) experience these artifacts, and vessels on the same side of the neck are relatively close to each

the pulsatile carotid flow phantom's bifurcation (a). For each voxel, it is shown: (b) time-velocity distribution derived from high-resolution 2DFT phase contrast; (c)

10.8 dB

10.3 dB

10.5 dB

10.8 dB

10.9 dB

11.7 dB

time

2

other. Therefore, the unaliased field-of-view is wide enough to enclose all vessels on one side of the neck, so that quantification of a vessel of interest will not be disturbed by flow from neighboring vessels (e.g. measurement of flow in the left carotid arteries will not be disturbed by flow in the left jugular vein). In order to suppress signal from the opposite side of the neck, we separately reconstruct data from the left and right neck-coil elements. This uses the receiver coil sensitivity profile to avoid spiral view-sharing artifacts. Based on these observations, we used spiral interleaf view-sharing in all subsequent carotid studies.

Fig. 5. Comparison of different view-orderings for multi-shot spiral FVE, in a healthy volunteer carotid study. When two or more *kv* levels are acquired during the same heartbeat (a), velocity distribution changes between consecutive TRs cause ghosting artifacts along the velocity axis (arrow). This artifact is not seen if, in the same heartbeat, different spiral interleaves, but only one *kv* encoding, are acquired (b).

#### **3.4.2 Comparison with Doppler ultrasound**

Representative *in vivo* results are compared with Doppler ultrasound in Figure 6. The MRI measured time-velocity histogram shows good agreement with the ultrasound measurement, as the peak velocity and the shape of the flow waveform were comparable to those observed in the ultrasound studies.

Fig. 6. Comparison of the spiral FVE method (a) with Doppler ultrasound (b), in a healthy volunteer common carotid artery study. Peak velocity and time-velocity waveforms are in good agreement.

**−80 −40 0 40 80**

<sup>525</sup> Assessment of Carotid Flow Using Magnetic

**<sup>0</sup> 10 20 30 −80**

**<sup>0</sup> 10 20 30 −80**

**spatial dimension (mm)**

**velocity profile**

**spatial dimension (mm)**

∑*v* |*s*(*v*)|

, the intra-voxel position *r* is incremented by a fraction

. (19)

. The volume fractions are

**velocity (cm/s)**

**velocity (cm/s)**

**velocity (cm/s)**

High spatial resolution is typically associated with low SNR (a). Averaging multiple acquisitions improves SNR (b), but also increases the total scan time, and may cause loss of effective resolution due to subject motion. Scan parameters: 0.33x0.33x3 mm3 spatial resolution, 37 ms temporal resolution, 30◦ flip angle, 80 cm/s Venc, 2-minute scan (120

Fig. 8. Illustration of *in vivo* high-resolution 2DFT phase contrast MRI, obtained at the carotid bifurcation of a healthy volunteer at peak flow: (a) single acquisition; (b) ten signal averages.

the use of spiral FVE to estimate WSR and the oscillatory shear index. This is made possible by the method proposed by Frayne & Rutt (1995) for FVE-based WSR estimation, which is

The Frayne method for FVE-based WSR estimation (Frayne & Rutt, 1995) involves obtaining the velocity distribution for a voxel spanning the blood/vessel wall interface, and then using this distribution to reconstruct the velocity profile across the voxel, with sub-voxel spatial resolution. Assuming that signal intensity is spatially and velocity invariant, the sum of the signals from all velocities within a voxel is proportional to the total volume of material within the voxel. Furthermore, two assumptions can be made about the shape of the velocity profile: (i) the fluid velocity at the vessel wall is approximately zero, and (ii) the velocity profile within a voxel is monotonically increasing or decreasing. Using these assumptions, a step-wise discrete approximation *v*˜(*r*) to the true velocity profile across a voxel can be obtained from its velocity distribution *s*(*v*) by inverting the discrete function *r*(*v*), which is constructed as

*<sup>r</sup>*(*vi*) = *<sup>r</sup>*(*vi*−<sup>1</sup> ) + <sup>Δ</sup>*<sup>r</sup>* · *<sup>h</sup>*(*vi*), where *<sup>h</sup>*(*vi*) = <sup>|</sup>*s*(*vi*)<sup>|</sup>

of the total radial extent of the voxel (Δ*r*). This fraction is proportional to *h*(*vi*), which

is the volume fraction of the voxel that has velocity *v* = *vi*

**magnitude map velocity map**

**jugular vein**

**carotid bifurcation**

Resonance Imaging and Computational Fluid Dynamics

**4.1 FVE-based WSR estimation: the Frayne method**

**a**

**b**

heartbeats) per acquisition.

discussed next.

follows:

Note that, for each velocity bin *vi*

#### **3.4.3 Measurement of multiple flows**

Figure 7 illustrates spiral FVE's ability to resolve different flows from a single dataset. A different time-velocity distribution was calculated for each voxel, and the distributions from select voxels are presented.

Fig. 7. Multiple flow distributions obtained from a single spiral FVE dataset. For each voxel in the image, a time-velocity distribution was calculated.

#### **4. Wall shear rate estimation using spiral FVE**

Arterial wall shear stress, the drag force acting on the endothelium as a result of blood flow, is widely believed to influence the formation and growth of atherosclerotic plaque, and may have prognostic value. According to Newton's law of viscosity, WSS can be estimated as the product of wall shear rate and blood viscosity (*μ*), where WSR is the radial gradient of blood flow velocity (*dv*/*dr*) at the vessel wall. Low WSS (Tang et al., 2008; Zarins et al., 1983) and highly oscillatory WSS (Ku et al., 1985) have been linked to the formation and growth of atherosclerotic plaques, and this link has been validated *in vitro* (Dai et al., 2004). High WSS has also been hypothesized as a factor responsible for the topography of atherosclerotic lesions (Thubrikar & Robicsek, 1995).

Phase contrast MRI (PC-MRI) suffers from data inconsistency, partial-volume effects (Tang et al., 1993), intravoxel phase dispersion and inadequate SNR at high spatial resolutions (Figure 8). PC-MRI is not currently capable of providing accurate absolute measurements of WSS (Boussel et al., 2009), and has been shown to underestimate blood flow velocities in the carotid bifurcation by 31–39% (Harloff et al., 2009).

As shown in section 3, spiral FVE produces accurate velocity histograms compared with those acquired with Doppler ultrasound. In addition, spiral FVE method is capable of rapid acquisition of fully localized, time-resolved velocity distributions. In this section, we propose 12 Will-be-set-by-IN-TECH

Figure 7 illustrates spiral FVE's ability to resolve different flows from a single dataset. A different time-velocity distribution was calculated for each voxel, and the distributions from

Fig. 7. Multiple flow distributions obtained from a single spiral FVE dataset. For each voxel

Arterial wall shear stress, the drag force acting on the endothelium as a result of blood flow, is widely believed to influence the formation and growth of atherosclerotic plaque, and may have prognostic value. According to Newton's law of viscosity, WSS can be estimated as the product of wall shear rate and blood viscosity (*μ*), where WSR is the radial gradient of blood flow velocity (*dv*/*dr*) at the vessel wall. Low WSS (Tang et al., 2008; Zarins et al., 1983) and highly oscillatory WSS (Ku et al., 1985) have been linked to the formation and growth of atherosclerotic plaques, and this link has been validated *in vitro* (Dai et al., 2004). High WSS has also been hypothesized as a factor responsible for the topography of atherosclerotic

Phase contrast MRI (PC-MRI) suffers from data inconsistency, partial-volume effects (Tang et al., 1993), intravoxel phase dispersion and inadequate SNR at high spatial resolutions (Figure 8). PC-MRI is not currently capable of providing accurate absolute measurements of WSS (Boussel et al., 2009), and has been shown to underestimate blood flow velocities in

As shown in section 3, spiral FVE produces accurate velocity histograms compared with those acquired with Doppler ultrasound. In addition, spiral FVE method is capable of rapid acquisition of fully localized, time-resolved velocity distributions. In this section, we propose

120

0

−120

120

−120

120

−120

120

0 500

0

0

0

0 500

0 500

0 500

−120

**3.4.3 Measurement of multiple flows**

velocity (cm/s)

120

0

0 500

−120

time (ms) <sup>120</sup>

0 500

0 500

0 500

lesions (Thubrikar & Robicsek, 1995).

in the image, a time-velocity distribution was calculated.

**4. Wall shear rate estimation using spiral FVE**

the carotid bifurcation by 31–39% (Harloff et al., 2009).

select voxels are presented.

−120

120

0

−120

120

0

−120

0

Fig. 8. Illustration of *in vivo* high-resolution 2DFT phase contrast MRI, obtained at the carotid bifurcation of a healthy volunteer at peak flow: (a) single acquisition; (b) ten signal averages. High spatial resolution is typically associated with low SNR (a). Averaging multiple acquisitions improves SNR (b), but also increases the total scan time, and may cause loss of effective resolution due to subject motion. Scan parameters: 0.33x0.33x3 mm3 spatial resolution, 37 ms temporal resolution, 30◦ flip angle, 80 cm/s Venc, 2-minute scan (120 heartbeats) per acquisition.

the use of spiral FVE to estimate WSR and the oscillatory shear index. This is made possible by the method proposed by Frayne & Rutt (1995) for FVE-based WSR estimation, which is discussed next.

#### **4.1 FVE-based WSR estimation: the Frayne method**

The Frayne method for FVE-based WSR estimation (Frayne & Rutt, 1995) involves obtaining the velocity distribution for a voxel spanning the blood/vessel wall interface, and then using this distribution to reconstruct the velocity profile across the voxel, with sub-voxel spatial resolution. Assuming that signal intensity is spatially and velocity invariant, the sum of the signals from all velocities within a voxel is proportional to the total volume of material within the voxel. Furthermore, two assumptions can be made about the shape of the velocity profile: (i) the fluid velocity at the vessel wall is approximately zero, and (ii) the velocity profile within a voxel is monotonically increasing or decreasing. Using these assumptions, a step-wise discrete approximation *v*˜(*r*) to the true velocity profile across a voxel can be obtained from its velocity distribution *s*(*v*) by inverting the discrete function *r*(*v*), which is constructed as follows:

$$r(v\_i) = r(v\_{i-1}) + \Delta r \cdot h(v\_i), \quad \text{where} \quad h(v\_i) = \frac{|s(v\_i)|}{\sum\_{v} |s(v)|}. \tag{19}$$

Note that, for each velocity bin *vi* , the intra-voxel position *r* is incremented by a fraction of the total radial extent of the voxel (Δ*r*). This fraction is proportional to *h*(*vi*), which is the volume fraction of the voxel that has velocity *v* = *vi* . The volume fractions are

**4.2** *In vivo* **WSR measurement**

Resonance Imaging and Computational Fluid Dynamics

observed a similar variation using a PC-based approach.

**4.3 Estimation of oscillatory shear index**

blood viscosity.

[*v*0,*v*1] interval.

The *in vivo* measurement of carotid WSR using spiral FVE acquisitions with Frayne's reconstruction is now demonstrated. Three healthy subjects were studied. For each volunteer, five 5 mm contiguous slices (2.5 cm coverage) were prescribed perpendicular to the left carotid bifurcation. Each slice was imaged independently (separate acquisitions). Localized gradient shimming was performed, and acquisitions were prospectively ECG-gated. Several cardiac phases were acquired, spanning the systolic portion of the cardiac cycle. A time-bandwidth product 2 RF pulse was used for excitation, and the flip angle was 30◦. Only through-plane velocities were measured, using 32 *kv* encoding steps over a 160 cm/s velocity field-of-view (5 cm/s resolution). The velocity field-of-view was shifted from the −80 to 80 cm/s range to the −40 to 120 cm/s range during reconstruction, in order to avoid aliasing and maximize field-of-view usage. Negative velocities are encoded in order to assess negative WSR values, and also to accommodate leakage and ringing due to *kv* truncation (i.e., finite velocity resolution). Spatial encoding was performed using eight 4 ms variable-density spiral interleaves (16∼4 cm field-of-view, 1.4 mm resolution). The temporal resolution was 24 ms (12 ms TR, 2 views per beat). Scan time was 128 heartbeats per slice, i.e., approximately 2 minutes per slice. The subjects provided informed consent, and were scanned using a protocol approved by the institutional review board of the University of Southern California. Figures 10 and 11 show two representative sets of *in vivo* results. The WSR values are shown for manually-segmented regions-of-interest. Figure 10 illustrates the variation in WSR along all three spatial dimensions near the carotid bifurcation of subject #1. These results correspond to the cardiac phase with the highest peak velocity. Figure 11 illustrates the temporal variation of pulsatile WSR in the common carotid artery of subject #2. The results show a circumferential variation in WSR around the wall of the common carotid artery. Markl et al. (2009) recently

<sup>527</sup> Assessment of Carotid Flow Using Magnetic

The OSI is important for the evaluation of shear stress imposed by pulsatile flow. This index describes the shear stress acting in directions other than the direction of the temporal mean shear stress vector (He & Ku, 1996; Ku et al., 1985). It is defined as the relation between the time-integral of the shear stress component acting in the direction opposite to the main direction of flow and the time-integral of the absolute shear stress. In this work, OSI was calculated as defined by He & Ku (1996), and assuming spatially and temporally invariant

Measuring negative WSR with the proposed method requires voxels to be small enough to contain only reverse flow. Voxels containing both forward and reverse flow would violate the assumption of a monotonically increasing/decreasing velocity profile within the voxel. Under sufficient spatial resolution, negative WSR may be measured simply by using a negative

Figure 12 presents a demonstration of *in vivo* OSI estimation, using the proposed method. Results are shown for eight representative voxels, selected around the circumference of the wall of the carotid bifurcation of subject #3. Measured velocity distributions, wall shear rates, and OSI values, corresponding to the systolic portion of the cardiac cycle, are shown for each voxel. High OSI values were observed at the wall corresponding to the carotid bulb. Non-zero OSI was also observed at the opposite wall. These findings are in agreement with PC-based OSI measurements recently reported by Markl et al. (2009). The results also suggest

calculated by normalizing the velocity distribution. This process is demonstrated graphically in Figure 9. Spatial variations in signal intensity due to radio-frequency saturation (i.e., inflow enhancement) and due to differences in 1H density and relaxation properties between vessel wall tissue and blood may be compensated by adjusting *s*(*v*) accordingly, prior to calculating *h*(*v*) (Frayne & Rutt, 1995).

Fig. 9. Graphical illustration of the Frayne method (Frayne & Rutt, 1995). WSR is estimated from FVE velocity distributions in voxels spanning the blood/vessel wall interface. First, a threshold is applied to the velocity histogram to reduce noise sensitivity (a). Then, the volume fraction within each velocity bin is converted into a radial position across the voxel, using Eq. 19. Finally, the velocity gradient is estimated from the reconstructed velocity profile (b), within a small velocity interval [*v*0,*v*1].

In order to reduce the effects of ringing and noise rectification due to the magnitude operation in Eq. 19, a threshold is applied to *s*(*v*) before normalization (Figure 9a). The appropriate threshold level must be determined by analyzing the signal intensities in the velocity distribution for a range of velocities outside the range of expected blood flow velocities. It is assumed that signal outside this expected range is exclusively due to rectified noise (Frayne & Rutt, 1995). In our implementation, only components that are below the specified threshold *and* outside this expected range of velocities are set to zero.

The WSR can be estimated by prescribing a velocity interval [*v*0,*v*1] and then fitting a first-order polynomial to the points of *v*˜(*r*) within this interval (Figure 9b). Ideally, *v*<sup>0</sup> = 0 and *v*<sup>1</sup> = Δ*v*, because we wish to estimate the velocity derivative at the blood-wall interface. The SNR of shear rate estimates will increase as this velocity interval becomes larger, because of averaging across multiple velocity steps. However, as the interval becomes larger, the shear rate is averaged over a larger distance within the voxel and may deviate from the true local shear at the wall (Frayne & Rutt, 1995). Therefore, it is important to prescribe a reasonable [*v*0,*v*1] interval. In our implementation, *v*<sup>0</sup> = Δ*v* and *v*<sup>1</sup> ≈ 30 cm/s are used for an initial assessment, and then the interval is manually adjusted for selected voxels of interest. The same voxel-based approach is used with respect to the noise threshold discussed above, with a fixed threshold value being used for the initial assessment. A negative [*v*0,*v*1] interval — e.g., *v*<sup>0</sup> = −5 cm/s and *v*<sup>1</sup> = −15 cm/s — is used when large volume fractions are measured on negative velocities, i.e., when there are large *h*(*v*) values for *v* < 0. This allows measurements of negative WSR values.

#### **4.2** *In vivo* **WSR measurement**

14 Will-be-set-by-IN-TECH

calculated by normalizing the velocity distribution. This process is demonstrated graphically in Figure 9. Spatial variations in signal intensity due to radio-frequency saturation (i.e., inflow enhancement) and due to differences in 1H density and relaxation properties between vessel wall tissue and blood may be compensated by adjusting *s*(*v*) accordingly, prior to calculating

threshold

**...**

Fig. 9. Graphical illustration of the Frayne method (Frayne & Rutt, 1995). WSR is estimated from FVE velocity distributions in voxels spanning the blood/vessel wall interface. First, a threshold is applied to the velocity histogram to reduce noise sensitivity (a). Then, the volume fraction within each velocity bin is converted into a radial position across the voxel, using Eq. 19. Finally, the velocity gradient is estimated from the reconstructed velocity

In order to reduce the effects of ringing and noise rectification due to the magnitude operation in Eq. 19, a threshold is applied to *s*(*v*) before normalization (Figure 9a). The appropriate threshold level must be determined by analyzing the signal intensities in the velocity distribution for a range of velocities outside the range of expected blood flow velocities. It is assumed that signal outside this expected range is exclusively due to rectified noise (Frayne & Rutt, 1995). In our implementation, only components that are below the specified threshold

The WSR can be estimated by prescribing a velocity interval [*v*0,*v*1] and then fitting a first-order polynomial to the points of *v*˜(*r*) within this interval (Figure 9b). Ideally, *v*<sup>0</sup> = 0 and *v*<sup>1</sup> = Δ*v*, because we wish to estimate the velocity derivative at the blood-wall interface. The SNR of shear rate estimates will increase as this velocity interval becomes larger, because of averaging across multiple velocity steps. However, as the interval becomes larger, the shear rate is averaged over a larger distance within the voxel and may deviate from the true local shear at the wall (Frayne & Rutt, 1995). Therefore, it is important to prescribe a reasonable [*v*0,*v*1] interval. In our implementation, *v*<sup>0</sup> = Δ*v* and *v*<sup>1</sup> ≈ 30 cm/s are used for an initial assessment, and then the interval is manually adjusted for selected voxels of interest. The same voxel-based approach is used with respect to the noise threshold discussed above, with a fixed threshold value being used for the initial assessment. A negative [*v*0,*v*1] interval — e.g., *v*<sup>0</sup> = −5 cm/s and *v*<sup>1</sup> = −15 cm/s — is used when large volume fractions are measured on negative velocities, i.e., when there are large *h*(*v*) values for *v* < 0. This allows measurements

−40

0

40

velocity (cm/s)

**b**

80

120

0 0.35 0.7 1.05 1.4 intra-voxel position (mm)

v1

v0

*h*(*v*) (Frayne & Rutt, 1995).

**a**

velocity (cm/s)

of negative WSR values.

−40

0

40

80

120

20 15 10 5 0 volume fraction (%)

profile (b), within a small velocity interval [*v*0,*v*1].

*and* outside this expected range of velocities are set to zero.

after threshold measured

The *in vivo* measurement of carotid WSR using spiral FVE acquisitions with Frayne's reconstruction is now demonstrated. Three healthy subjects were studied. For each volunteer, five 5 mm contiguous slices (2.5 cm coverage) were prescribed perpendicular to the left carotid bifurcation. Each slice was imaged independently (separate acquisitions). Localized gradient shimming was performed, and acquisitions were prospectively ECG-gated. Several cardiac phases were acquired, spanning the systolic portion of the cardiac cycle. A time-bandwidth product 2 RF pulse was used for excitation, and the flip angle was 30◦. Only through-plane velocities were measured, using 32 *kv* encoding steps over a 160 cm/s velocity field-of-view (5 cm/s resolution). The velocity field-of-view was shifted from the −80 to 80 cm/s range to the −40 to 120 cm/s range during reconstruction, in order to avoid aliasing and maximize field-of-view usage. Negative velocities are encoded in order to assess negative WSR values, and also to accommodate leakage and ringing due to *kv* truncation (i.e., finite velocity resolution). Spatial encoding was performed using eight 4 ms variable-density spiral interleaves (16∼4 cm field-of-view, 1.4 mm resolution). The temporal resolution was 24 ms (12 ms TR, 2 views per beat). Scan time was 128 heartbeats per slice, i.e., approximately 2 minutes per slice. The subjects provided informed consent, and were scanned using a protocol approved by the institutional review board of the University of Southern California. Figures 10 and 11 show two representative sets of *in vivo* results. The WSR values are shown for manually-segmented regions-of-interest. Figure 10 illustrates the variation in WSR along all three spatial dimensions near the carotid bifurcation of subject #1. These results correspond to the cardiac phase with the highest peak velocity. Figure 11 illustrates the temporal variation of pulsatile WSR in the common carotid artery of subject #2. The results show a circumferential variation in WSR around the wall of the common carotid artery. Markl et al. (2009) recently observed a similar variation using a PC-based approach.

#### **4.3 Estimation of oscillatory shear index**

The OSI is important for the evaluation of shear stress imposed by pulsatile flow. This index describes the shear stress acting in directions other than the direction of the temporal mean shear stress vector (He & Ku, 1996; Ku et al., 1985). It is defined as the relation between the time-integral of the shear stress component acting in the direction opposite to the main direction of flow and the time-integral of the absolute shear stress. In this work, OSI was calculated as defined by He & Ku (1996), and assuming spatially and temporally invariant blood viscosity.

Measuring negative WSR with the proposed method requires voxels to be small enough to contain only reverse flow. Voxels containing both forward and reverse flow would violate the assumption of a monotonically increasing/decreasing velocity profile within the voxel. Under sufficient spatial resolution, negative WSR may be measured simply by using a negative [*v*0,*v*1] interval.

Figure 12 presents a demonstration of *in vivo* OSI estimation, using the proposed method. Results are shown for eight representative voxels, selected around the circumference of the wall of the carotid bifurcation of subject #3. Measured velocity distributions, wall shear rates, and OSI values, corresponding to the systolic portion of the cardiac cycle, are shown for each voxel. High OSI values were observed at the wall corresponding to the carotid bulb. Non-zero OSI was also observed at the opposite wall. These findings are in agreement with PC-based OSI measurements recently reported by Markl et al. (2009). The results also suggest

that higher spatial resolution is needed for accurately estimating OSI in some of the voxels. Notably, voxel (h) presents both positive and negative velocity components simultaneously during post-systolic deceleration. This indicates that the spatial resolution was insufficient, and the assumption of a monotonically decreasing velocity profile within the voxel was

<sup>529</sup> Assessment of Carotid Flow Using Magnetic

b

OSI=0.45

OSI=0.00

OSI=0.00

OSI=0.24

d

h

OSI=0.00

OSI=0.00

c

g

f

OSI=0.48

Fig. 12. *In vivo* assessment of oscillatory shear index. Results are shown for eight representative voxels (a–h), selected around the circumference of the wall of the carotid bifurcation of subject #3 (see inset). Measured velocity distributions, wall shear rates, and OSI values, corresponding to the systolic portion of the cardiac cycle, are shown for each voxel.

Computational fluid dynamics methods are concerned with the approximate solution of the fluid motion as well as with the interaction of the fluid with solid bodies. Fluid dynamic problems, essentially, are based on nonlinear systems of coupled partial differential equations. Attempting to solve those systems analytically, in general, is an impracticable task. Over the years, researchers have developed methodologies, codes, schemes and algorithms to find approximate solutions, focusing on accuracy and speed of convergence, for problems

Historically, CFD started in the early 1970s, triggered by the availability of increasingly more powerful mainframes. In the beginning, the study was limited to high-technology engineering areas of aeronautics and astronautics. Nowadays, computational fluid dynamics methodologies are routinely employed in many fields such as: racing car design, ship design, meteorology, oil recovery, civil engineering, airspace engineering, and biomedical

e

**5. Assessment of carotid flow using CFD (and MRI)**

involving fluid dynamics and heat transfer.

engineering.

<sup>50</sup> <sup>175</sup> <sup>300</sup> <sup>−</sup><sup>40</sup>

OSI=0.44

<sup>50</sup> <sup>175</sup> <sup>300</sup> <sup>−</sup><sup>850</sup>

time(ms)

0 40 80

velocity (cm/s)

Resonance Imaging and Computational Fluid Dynamics

<sup>c</sup> <sup>d</sup> e

f g h

a b

shear rate (s−1)

a

0 850 1700

violated.

Fig. 10. Carotid WSR measured across the carotid bifurcation of subject #1: (a) slice prescription; (b) spiral FVE WSR measurements. Results correspond to the cardiac phase with the highest peak velocity (96 ms after the ECG trigger). The common (CCA), external (ECA), and internal (ICA) carotid arteries, and the jugular vein (JV), are indicated.

Fig. 11. Temporal variation of pulsatile WSR measured in the common carotid artery of subject #2. Results correspond to the cardiac phases acquired 48–168 ms after the ECG trigger, and to a slice prescribed 10 mm below the carotid bifurcation.

16 Will-be-set-by-IN-TECH

**b**

**a**

−10 mm

Fig. 10. Carotid WSR measured across the carotid bifurcation of subject #1: (a) slice prescription; (b) spiral FVE WSR measurements. Results correspond to the cardiac phase with the highest peak velocity (96 ms after the ECG trigger). The common (CCA), external

(ECA), and internal (ICA) carotid arteries, and the jugular vein (JV), are indicated.

48 ms 72 ms 96 ms 120 ms 144 ms 168 ms

Fig. 11. Temporal variation of pulsatile WSR measured in the common carotid artery of subject #2. Results correspond to the cardiac phases acquired 48–168 ms after the ECG

time

cardiac cycle

trigger, and to a slice prescribed 10 mm below the carotid bifurcation.

velocity

CCA

JV

ICA

−5 mm

0 mm

5 mm

10 mm

ECA

200

400

600

800

1000

wall shear rate (s

−1)

2000

1000

wall shear rate (s

−1

)

1500

500

0

1200

1400

1600

1800

2000

that higher spatial resolution is needed for accurately estimating OSI in some of the voxels. Notably, voxel (h) presents both positive and negative velocity components simultaneously during post-systolic deceleration. This indicates that the spatial resolution was insufficient, and the assumption of a monotonically decreasing velocity profile within the voxel was violated.

Fig. 12. *In vivo* assessment of oscillatory shear index. Results are shown for eight representative voxels (a–h), selected around the circumference of the wall of the carotid bifurcation of subject #3 (see inset). Measured velocity distributions, wall shear rates, and OSI values, corresponding to the systolic portion of the cardiac cycle, are shown for each voxel.

#### **5. Assessment of carotid flow using CFD (and MRI)**

Computational fluid dynamics methods are concerned with the approximate solution of the fluid motion as well as with the interaction of the fluid with solid bodies. Fluid dynamic problems, essentially, are based on nonlinear systems of coupled partial differential equations. Attempting to solve those systems analytically, in general, is an impracticable task. Over the years, researchers have developed methodologies, codes, schemes and algorithms to find approximate solutions, focusing on accuracy and speed of convergence, for problems involving fluid dynamics and heat transfer.

Historically, CFD started in the early 1970s, triggered by the availability of increasingly more powerful mainframes. In the beginning, the study was limited to high-technology engineering areas of aeronautics and astronautics. Nowadays, computational fluid dynamics methodologies are routinely employed in many fields such as: racing car design, ship design, meteorology, oil recovery, civil engineering, airspace engineering, and biomedical engineering.

volunteer. The results show that this methodology produces flow fields that are in better

<sup>531</sup> Assessment of Carotid Flow Using Magnetic

Here, blood is modeled as an incompressible Newtonian fluid with constant viscosity *μ*. This model is widely used in *in vivo* CFD analysis, and underlies the commonly held definition of wall shear stress as being proportional to *dv*/*dr*, where *v* is the velocity tangential to the vessel wall, and *r* is the perpendicular distance from the wall. The task of a CFD routine is to solve the momentum and continuity equations that the flow field — subject to the assumption of Newtonian flow — must obey. Following the control-volume formulation introduced by

where *ρ* is the fluid density, **v** = (*u*, *v*, *w*) is the velocity vector field, *p* is the pressure field, ∇ is the gradient operator, and Δ is the Laplacian operator. The flow field must also satisfy mass

Equations (20) and (21) must be solved for the unknown scalar field variables *u*, *v*, *w*, and *p*. These equations are non-linear and coupled, and attempting to solve them directly in one step

The solver was built on the SIMPLER algorithm developed by Patankar (1980), which is a well-known and established numerical routine for solving Eqs. (20) and (21). SIMPLER starts with an initial estimate for the velocity field, and updates this estimate in an iterative fashion. At each iteration *<sup>i</sup>*, the discretized equations are linearized using the velocity estimate **<sup>v</sup>***i*−<sup>1</sup> at the previous iteration, which produces a square system matrix **A** for each velocity component.

where **w** is the *z* velocity in all voxels in the 3D calculation domain (stacked to form a 1D column vector), and **b** is a constant column vector. Note that the pressure field is updated periodically based on the current velocity field estimate, and hence does not appear explicitly

The key step in our approach is to add additional rows to **A** and **b** that incorporate MRI measurements of one or more velocity components. The underlying assumption we make here is that the velocity measured with PC-MRI is equal to the average velocity within the voxel, and can hence be expressed as a linear combination of the velocities on the underlying

CFD calculation grid. For example, for the *z* (S/I) velocity component we have

where **w***MR* is the *z* velocity component measured with PC-MRI.

*Dt* <sup>=</sup> −∇*<sup>p</sup>* <sup>+</sup> *<sup>μ</sup>*Δ**v**, (20)

∇ · **v** = 0. (21)

**<sup>A</sup>***i*−1**w***<sup>i</sup>* = **<sup>b</sup>***i*−1, (22)

**w***MR* ≈ **A***MR***w**, (23)

agreement with direct PC-MRI measurements than CFD alone.

Resonance Imaging and Computational Fluid Dynamics

*ρ D***v**

**5.2.1 Proposed approach**

Patankar (1980), the momentum equation is

is a formidable, if not impossible, task.

in (22).

conservation (or continuity), which can be expressed as

For example, for the *z* velocity component, we can write

This section is dedicated to the use of CFD for evaluation of blood flow in the carotid arteries. *In vivo* 3D blood flow patterns can be either measured directly using PC-MRI, or obtained from model-based CFD calculations. PC-MRI is accurate, but suffers from long scan times and limited spatio-temporal resolution and SNR. CFD provides arbitrarily high resolution and reduced scan times, but its accuracy hinges on the model assumptions. A numerical framework for constructing a flow field that is influenced by both direct measurements and a fluid physics model is described.

#### **5.1 CFD in biomedical engineering**

Recently, CFD is playing an important role in the analysis of blood flow. This technique of flow visualization has been widely applied in problems involving arterial diseases with reconstruction of hemodynamics in realistic models based on images generated by standard *in vivo* medical visualization tools, such as MRI. CFD can be used to improve the data obtained using MRI for estimation of the flow properties of blood vessels. CFD can also be used to compare real data obtained from patients with simulated data models using realistic geometries. These geometries may even be constructed using MRI data.

Boussel et al. (2009) compared a time-dependent 3D phase-contrast MRI sequence with patient-specific CFD models for patients who had intracranial aneurysms. The evolution of intracranial aneurysms is known to be related to hemodynamic forces such as wall shear stress and maximum shear stress. Harloff et al. (2010) used CFD image-based modeling as an option to improve the accuracy of MRI-based WSS and OSI estimation, with the purpose of correlating atherogenic low WSS and high OSI with the localization of aortic plaques. Steinman et al. (2002) used a novel approach for noninvasively reconstructing artery wall thickness and local hemodynamics at the human carotid bifurcation. Three-dimensional models of the lumen and wall boundaries, from which wall thickness can be measured, were reconstructed from black blood MRI. Along with time-varying inlet/outlet flow rates measured via PC-MRI, the lumen boundary was used as input for a CFD simulation of the subject-specific flow patterns and wall shear stress. Long et al. (2003) were concerned with the reproducibility of geometry reconstruction, one of the most crucial steps in the modeling process. Canstein et al. (2008) used rapid prototyping to transform aortic geometries as measured by contrast-enhanced MR angiography into realistic vascular models with large anatomical coverage. Visualization of characteristic 3D flow patterns and quantitative comparisons of the *in vitro* experiments with *in vivo* data and CFD simulations in identical vascular geometries were performed to evaluate the accuracy of vascular model systems.

CFD can also be used with other imaging techniques such as tomography. Howell et al. (2007) used CFD to study the temporal and spatial variations in surface pressure and shear through the cardiac cycle on models of bifurcated stent-grafts derived from computed tomography in patients who had previously undergone endovascular repair of abdominal aortic aneurysm.

#### **5.2 MRI-driven CFD**

A numerical framework for constructing a flow field that is influenced by both direct measurements and a fluid physics model is now described. The PC-MRI signal is expressed as a linear function of the velocities on the underlying computational grid, and a tunable parameter controls the relative influence of direct measurements and the model assumptions. The feasibility of the proposed approach is demonstrated in the carotid bifurcation of a healthy volunteer. The results show that this methodology produces flow fields that are in better agreement with direct PC-MRI measurements than CFD alone.

#### **5.2.1 Proposed approach**

18 Will-be-set-by-IN-TECH

This section is dedicated to the use of CFD for evaluation of blood flow in the carotid arteries. *In vivo* 3D blood flow patterns can be either measured directly using PC-MRI, or obtained from model-based CFD calculations. PC-MRI is accurate, but suffers from long scan times and limited spatio-temporal resolution and SNR. CFD provides arbitrarily high resolution and reduced scan times, but its accuracy hinges on the model assumptions. A numerical framework for constructing a flow field that is influenced by both direct measurements and a

Recently, CFD is playing an important role in the analysis of blood flow. This technique of flow visualization has been widely applied in problems involving arterial diseases with reconstruction of hemodynamics in realistic models based on images generated by standard *in vivo* medical visualization tools, such as MRI. CFD can be used to improve the data obtained using MRI for estimation of the flow properties of blood vessels. CFD can also be used to compare real data obtained from patients with simulated data models using realistic

Boussel et al. (2009) compared a time-dependent 3D phase-contrast MRI sequence with patient-specific CFD models for patients who had intracranial aneurysms. The evolution of intracranial aneurysms is known to be related to hemodynamic forces such as wall shear stress and maximum shear stress. Harloff et al. (2010) used CFD image-based modeling as an option to improve the accuracy of MRI-based WSS and OSI estimation, with the purpose of correlating atherogenic low WSS and high OSI with the localization of aortic plaques. Steinman et al. (2002) used a novel approach for noninvasively reconstructing artery wall thickness and local hemodynamics at the human carotid bifurcation. Three-dimensional models of the lumen and wall boundaries, from which wall thickness can be measured, were reconstructed from black blood MRI. Along with time-varying inlet/outlet flow rates measured via PC-MRI, the lumen boundary was used as input for a CFD simulation of the subject-specific flow patterns and wall shear stress. Long et al. (2003) were concerned with the reproducibility of geometry reconstruction, one of the most crucial steps in the modeling process. Canstein et al. (2008) used rapid prototyping to transform aortic geometries as measured by contrast-enhanced MR angiography into realistic vascular models with large anatomical coverage. Visualization of characteristic 3D flow patterns and quantitative comparisons of the *in vitro* experiments with *in vivo* data and CFD simulations in identical vascular geometries were performed to evaluate the accuracy of vascular model systems. CFD can also be used with other imaging techniques such as tomography. Howell et al. (2007) used CFD to study the temporal and spatial variations in surface pressure and shear through the cardiac cycle on models of bifurcated stent-grafts derived from computed tomography in patients who had previously undergone endovascular repair of abdominal aortic aneurysm.

A numerical framework for constructing a flow field that is influenced by both direct measurements and a fluid physics model is now described. The PC-MRI signal is expressed as a linear function of the velocities on the underlying computational grid, and a tunable parameter controls the relative influence of direct measurements and the model assumptions. The feasibility of the proposed approach is demonstrated in the carotid bifurcation of a healthy

geometries. These geometries may even be constructed using MRI data.

fluid physics model is described.

**5.1 CFD in biomedical engineering**

**5.2 MRI-driven CFD**

Here, blood is modeled as an incompressible Newtonian fluid with constant viscosity *μ*. This model is widely used in *in vivo* CFD analysis, and underlies the commonly held definition of wall shear stress as being proportional to *dv*/*dr*, where *v* is the velocity tangential to the vessel wall, and *r* is the perpendicular distance from the wall. The task of a CFD routine is to solve the momentum and continuity equations that the flow field — subject to the assumption of Newtonian flow — must obey. Following the control-volume formulation introduced by Patankar (1980), the momentum equation is

$$
\rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \mu \Delta \mathbf{v},\tag{20}
$$

where *ρ* is the fluid density, **v** = (*u*, *v*, *w*) is the velocity vector field, *p* is the pressure field, ∇ is the gradient operator, and Δ is the Laplacian operator. The flow field must also satisfy mass conservation (or continuity), which can be expressed as

$$
\nabla \cdot \mathbf{v} = 0.\tag{21}
$$

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

The solver was built on the SIMPLER algorithm developed by Patankar (1980), which is a well-known and established numerical routine for solving Eqs. (20) and (21). SIMPLER starts with an initial estimate for the velocity field, and updates this estimate in an iterative fashion. At each iteration *<sup>i</sup>*, the discretized equations are linearized using the velocity estimate **<sup>v</sup>***i*−<sup>1</sup> at the previous iteration, which produces a square system matrix **A** for each velocity component. For example, for the *z* velocity component, we can write

$$\mathbf{A}\_{i-1}\mathbf{w}\_i = \mathbf{b}\_{i-1} \tag{22}$$

where **w** is the *z* velocity in all voxels in the 3D calculation domain (stacked to form a 1D column vector), and **b** is a constant column vector. Note that the pressure field is updated periodically based on the current velocity field estimate, and hence does not appear explicitly in (22).

The key step in our approach is to add additional rows to **A** and **b** that incorporate MRI measurements of one or more velocity components. The underlying assumption we make here is that the velocity measured with PC-MRI is equal to the average velocity within the voxel, and can hence be expressed as a linear combination of the velocities on the underlying CFD calculation grid. For example, for the *z* (S/I) velocity component we have

$$\mathbf{w}\_{MR} \approx \mathbf{A}\_{MR} \mathbf{w}\_{\prime} \tag{23}$$

where **w***MR* is the *z* velocity component measured with PC-MRI.

Figure 13 compares flow fields in the carotid artery obtained with PC-MRI only, CFD only (*s* = 0), and the proposed combined solver algorithm with *s* = 5. These flow fields were calculated from 4, 1, or 2 time-resolved 3D image acquisitions, corresponding to a total scan time of 28, 7, and 14 minutes, as indicated in the figure. In the common carotid artery, flow is predominantly along *z*, and all methods produce comparable flow fields. In the bifurcation, the bulk flow pattern appears qualitatively similar for all three methods, which indicates that the underlying fluid physics model makes reasonable predictions regarding the transverse velocities. However, comparison of the CFD and PC-MRI results shows that CFD underestimates the velocities in the external carotid artery. The combined solver, which strikes a compromise between CFD and PC-MRI, brings the calculated flow field in closer agreement

<sup>533</sup> Assessment of Carotid Flow Using Magnetic

Fig. 13. Blood flow in the carotid bifurcation (top) and the common carotid artery (bottom) obtained with 4-point phase contrast MRI (left), CFD (center), and the combined solver (*s* = 5, right). The lines show the path of massless particles during the course of 60 ms (top) or 20 ms (bottom), under the assumption of a constant velocity field. Pathlines are RGB color-coded according to the local velocity direction (vertical=blue; in-plane=red and green).

corresponding to total scan times of 28, 7 and 14 minutes, respectively. Compared to CFD, the combined solver produces flow fields that are in better qualitative and quantitative

In this chapter, we have introduced spiral FVE, a rapid MRI method for fully-localized measurement of cardiovascular blood flow. The proposed method was shown to be capable of measuring blood flow in the carotid arteries and estimating wall shear stress and oscillatory shear index at the carotid bifurcation. In addition, we proposed a combined CFD-MRI solver that integrates the non-linear coupled system of the fluid partial differential equations using MRI data as initial data. This methodology, which has a tunable parameter, is capable of

These velocity fields were obtained from 4 (left), 1 (center), or 2 (right) MRI scans,

with the values measured with PC-MRI.

Resonance Imaging and Computational Fluid Dynamics

agreement with PC-MRI.

**6. Conclusion**

Combining Eqs. (22) and (23), we have

$$
\begin{bmatrix}
\mathbf{A}\_{i-1} \\
\mathbf{s}\mathbf{A}\_{MR}
\end{bmatrix}\mathbf{w}\_{i} = \begin{bmatrix}
\mathbf{b}\_{i-1} \\
\mathbf{s}\mathbf{w}\_{MR}
\end{bmatrix}.
\tag{24}
$$

The weighting parameter *s* is an adjustable scalar that determines the influence of the MRI measurements on the solution. For example, *s* = 0 produces a conventional, unconstrained CFD solution. Here, we solve Eq. (24) in the least squares sense, using the conjugate gradient method on the normal equations. Hence, the solution **w** at each iteration represents a weighted least squares estimate, with relative contributions of fluid physics and direct MRI measurement being controlled by the parameters after multiple interactions. The velocity field **v** = (*u*, *v*, *w*) converges toward the solution for a particular point *t*. Time is then incremented by an amount Δ*t*, and the iterative procedure is repeated to obtain the solution at *t* + Δ*t*. This way, it is possible to predict a time-dependent flow field, as long as the true velocity field **v** at *t* = 0 is known. In this work, however, the exact velocity field is not known, since only a relatively low-resolution and noisy PC-MRI velocity field is available. Instead, we will calculate the steady flow **v***<sup>s</sup>* subject to the measured inlet and outlet velocities. We will furthermore assume that the velocity field obtained with PC-MRI at one time-point (near peak flow) is representative of the steady flow given the inlet and outlet velocities at that time-point. We obtain **v***<sup>s</sup>* by starting with an initial guess for **v**, and carrying the simulations forward in time until convergence.

#### **5.2.2** *In vivo* **demonstration**

PC-MRI data were obtained from four time-resolved 3DFT FGRE image volumes in the carotid artery in one healthy volunteer (1×1×2.5 mm3 voxel size; FOV 16×12×7.5 cm3; TR 7.0 ms; flip angle 15◦; temporal resolution 56 ms; Venc 1.6 m/s; 7 minutes per scan), on a GE Signa 3T EXCITE HD system (4 G/cm and 15 G/cm/ms maximum gradient amplitude and slew rate) with a 4-channel carotid receive coil array. The through-slab (*z*) axis was oriented along the S/I direction. The subject provided informed consent, and was scanned using a protocol approved by the institutional review board of the University of Southern California. PC-MRI velocity component maps **u***MR*, **v***MR*, and **w***MR* were calculated using data from one receive coil. 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 vessel lumen was segmented by manually outlining the vessel borders from a stack of 2D axial slices.

The solver calculations assumed a blood viscosity of 0.0027 Pa·sec and a blood density of 1060 kg/m<sup>3</sup> (Reynolds number of order 1000), and no-slip boundary conditions. Calculations were performed on a Cartesian grid of 1 mm isotropic resolution, and all algorithms were implemented in Matlab. Low-resolution (truncated to 1×1×3 mm3) measurements of the S/I velocity component **w***MR* near the time-point of peak flow were incorporated into the solver, and a steady velocity vector field **v***<sup>s</sup>* was calculated as described above. Hence, **w** was partially constrained by the measured velocity component **w***MR*, whereas **u** and **v** were determined solely from the fluid physics model. As in the original SIMPLER algorithm (Patankar, 1980), *u*, *v*, and *w* were defined on regular grids that were staggered by half a grid spacing (in different directions) with respect to the centered grid. This is done to avoid a checkerboard solution for the pressure and velocity fields (Patankar, 1980). The calculation domain was rectangular of size 20.5×16.5×25 mm3.

20 Will-be-set-by-IN-TECH

The weighting parameter *s* is an adjustable scalar that determines the influence of the MRI measurements on the solution. For example, *s* = 0 produces a conventional, unconstrained CFD solution. Here, we solve Eq. (24) in the least squares sense, using the conjugate gradient method on the normal equations. Hence, the solution **w** at each iteration represents a weighted least squares estimate, with relative contributions of fluid physics and direct MRI measurement being controlled by the parameters after multiple interactions. The velocity field **v** = (*u*, *v*, *w*) converges toward the solution for a particular point *t*. Time is then incremented by an amount Δ*t*, and the iterative procedure is repeated to obtain the solution at *t* + Δ*t*. This way, it is possible to predict a time-dependent flow field, as long as the true velocity field **v** at *t* = 0 is known. In this work, however, the exact velocity field is not known, since only a relatively low-resolution and noisy PC-MRI velocity field is available. Instead, we will calculate the steady flow **v***<sup>s</sup>* subject to the measured inlet and outlet velocities. We will furthermore assume that the velocity field obtained with PC-MRI at one time-point (near peak flow) is representative of the steady flow given the inlet and outlet velocities at that time-point. We obtain **v***<sup>s</sup>* by starting with an initial guess for **v**, and carrying the simulations forward in

PC-MRI data were obtained from four time-resolved 3DFT FGRE image volumes in the carotid artery in one healthy volunteer (1×1×2.5 mm3 voxel size; FOV 16×12×7.5 cm3; TR 7.0 ms; flip angle 15◦; temporal resolution 56 ms; Venc 1.6 m/s; 7 minutes per scan), on a GE Signa 3T EXCITE HD system (4 G/cm and 15 G/cm/ms maximum gradient amplitude and slew rate) with a 4-channel carotid receive coil array. The through-slab (*z*) axis was oriented along the S/I direction. The subject provided informed consent, and was scanned using a protocol approved by the institutional review board of the University of Southern California. PC-MRI velocity component maps **u***MR*, **v***MR*, and **w***MR* were calculated using data from one receive coil. 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 vessel lumen was segmented by manually outlining

The solver calculations assumed a blood viscosity of 0.0027 Pa·sec and a blood density of 1060 kg/m<sup>3</sup> (Reynolds number of order 1000), and no-slip boundary conditions. Calculations were performed on a Cartesian grid of 1 mm isotropic resolution, and all algorithms were implemented in Matlab. Low-resolution (truncated to 1×1×3 mm3) measurements of the S/I velocity component **w***MR* near the time-point of peak flow were incorporated into the solver, and a steady velocity vector field **v***<sup>s</sup>* was calculated as described above. Hence, **w** was partially constrained by the measured velocity component **w***MR*, whereas **u** and **v** were determined solely from the fluid physics model. As in the original SIMPLER algorithm (Patankar, 1980), *u*, *v*, and *w* were defined on regular grids that were staggered by half a grid spacing (in different directions) with respect to the centered grid. This is done to avoid a checkerboard solution for the pressure and velocity fields (Patankar, 1980). The calculation domain was rectangular of

 **<sup>b</sup>***i*−<sup>1</sup> *s***w***MR* . (24)

 **<sup>A</sup>***i*−<sup>1</sup> *s***A***MR*  **w***<sup>i</sup>* =

Combining Eqs. (22) and (23), we have

time until convergence.

**5.2.2** *In vivo* **demonstration**

size 20.5×16.5×25 mm3.

the vessel borders from a stack of 2D axial slices.

Figure 13 compares flow fields in the carotid artery obtained with PC-MRI only, CFD only (*s* = 0), and the proposed combined solver algorithm with *s* = 5. These flow fields were calculated from 4, 1, or 2 time-resolved 3D image acquisitions, corresponding to a total scan time of 28, 7, and 14 minutes, as indicated in the figure. In the common carotid artery, flow is predominantly along *z*, and all methods produce comparable flow fields. In the bifurcation, the bulk flow pattern appears qualitatively similar for all three methods, which indicates that the underlying fluid physics model makes reasonable predictions regarding the transverse velocities. However, comparison of the CFD and PC-MRI results shows that CFD underestimates the velocities in the external carotid artery. The combined solver, which strikes a compromise between CFD and PC-MRI, brings the calculated flow field in closer agreement with the values measured with PC-MRI.

Fig. 13. Blood flow in the carotid bifurcation (top) and the common carotid artery (bottom) obtained with 4-point phase contrast MRI (left), CFD (center), and the combined solver (*s* = 5, right). The lines show the path of massless particles during the course of 60 ms (top) or 20 ms (bottom), under the assumption of a constant velocity field. Pathlines are RGB color-coded according to the local velocity direction (vertical=blue; in-plane=red and green). These velocity fields were obtained from 4 (left), 1 (center), or 2 (right) MRI scans, corresponding to total scan times of 28, 7 and 14 minutes, respectively. Compared to CFD, the combined solver produces flow fields that are in better qualitative and quantitative agreement with PC-MRI.

#### **6. Conclusion**

In this chapter, we have introduced spiral FVE, a rapid MRI method for fully-localized measurement of cardiovascular blood flow. The proposed method was shown to be capable of measuring blood flow in the carotid arteries and estimating wall shear stress and oscillatory shear index at the carotid bifurcation. In addition, we proposed a combined CFD-MRI solver that integrates the non-linear coupled system of the fluid partial differential equations using MRI data as initial data. This methodology, which has a tunable parameter, is capable of

Howell, B. A., Kim, T., Cheer, A., Dwyer, H., Saloner, D. & Chuter, T. A. M. (2007).

<sup>535</sup> Assessment of Carotid Flow Using Magnetic

Hu, B. S., Pauly, J. M. & Macovski, A. (1993). Localized real-time velocity spectra

Kim, C. S., Kiris, C., Kwak, D. & David, T. (2006). Numerical simulation of local blood flow in

Ku, D. N., Giddens, D. P., Zarins, C. K. & Glagov, S. (1985). Pulsatile flow and atherosclerosis

Long, Q., Ariff, B., Zhao, S. Z., Thom, S. A., Hughes, A. D. & Xu, X. Y. (2003). Reproducibility

Markl, M., Wegent, F., Bauer, S., Stalder, A. F., Frydrychowicz, A., Weiller, C., Schumacher, M.

Moran, P. R., Moran, R. A. & Karstaedt, N. (1985). Verification and evaluation of internal

Nayler, G. L., Firmin, D. N. & Longmore, D. B. (1986). Blood flow imaging by cine magnetic

O'Donnell, M. (1985). NMR blood flow imaging using multiecho, phase contrast sequences,

Papathanasopoulou, P., Zhao, S., Köhler, U., Robertson, M. B., Long, Q., Hoskins, P., Xu, X. Y.

Patankar, S. V. (1980). *Numerical Heat Transfer and Fluid Flow*, Hemisphere Publishing

Rebergen, S. A., van der Wall, E. E., Doornbos, J. & de Roos, A. (1993). Magnetic resonance

Singer, J. R. (1959). Blood flow rates by nuclear magnetic resonance measurements, *Science*

Singer, J. R. & Crooks, L. E. (1983). Nuclear magnetic resonance blood flow measurements in

Steinman, D. A., Thomas, J. B., Ladak, H. M., Milner, J. S., Rutt, B. K. & Spence, J. D.

computational fluid dynamics and MRI, *Magn Reson Med* 47(1): 149–159. Tang, C., Blatter, D. D. & Parker, D. L. (1993). Accuracy of phase-contrast flow measurements in the presence of partial-volume effects, *J Magn Reson Imaging* 3(2): 377–385. Tang, D., Yang, C., Mondal, S., Liu, F., Canton, G., Hatsukami, T. & Yuan, C. (2008). A negative

wall stress: in vivo MRI-based 2D/3D FSI models, *J Biomech* 41(4): 727–736.

low oscillating shear stress, *Arterioscler Thromb Vasc Biol* 5: 293–302.

magnetic resonance images, *Magn Reson Med* 49(4): 665–674.

*Endovasc Ther* 14(2): 138–143.

Resonance Imaging and Computational Fluid Dynamics

*Magn Reson Imaging* 1(4): 197–203.

method, *Radiology* 154(2): 433–441.

*Med Phys* 12(1): 59–64.

*Imaging* 17(2): 153–162.

130(3389): 1652–1653.

Corporation.

resonance, *J Comput Assist Tomogr* 10(5): 715–722.

applications, *Am Heart J* 126(6): 1439–1456.

the human brain, *Science* 221(4611): 654–656.

determination, *Magn Reson Med* 30(3): 393–398.

Computational fluid dynamics within bifurcated abdominal aortic stent-grafts, *J*

the carotid and cerebral arteries under altered gravity, *J Biomech Eng* 128(2): 194–202.

in the human carotid bifurcation. Positive correlation between plaque location and

study of 3D geometrical reconstruction of the human carotid bifurcation from

& Harloff, A. (2009). Three-dimensional assessment of wall shear stress distribution in the carotid bifurcation, *Proc, ISMRM, 17th Annual Meeting*, Honolulu, p. 323. Moran, P. R. (1982). A flow velocity zeugmatographic interlace for NMR imaging in humans,

flow and motion. true magnetic resonance imaging by the phase gradient modulation

& Marshall, I. (2003). MRI measurement of time-resolved wall shear stress vectors in a carotid bifurcation model, and comparison with CFD predictions, *J Magn Reson*

measurement of velocity and flow: technique, validation, and cardiovascular

(2002). Reconstruction of carotid bifurcation hemodynamics and wall thickness using

correlation between human carotid atherosclerotic plaque progression and plaque

producing flow fields that are in better agreement with direct PC-MRI measurements than CFD alone.

MRI is potentially the most appropriate technique for addressing all aspects of a complete cardiovascular disease examination. The evaluation of carotid flow will be a necessary capability in such examination. The methodologies presented in this chapter can improve the quality of the MRI measured data, reducing scan time, and improving the SNR. These techniques should improve the diagnosis and understanding of carotid diseases.

#### **7. References**


22 Will-be-set-by-IN-TECH

producing flow fields that are in better agreement with direct PC-MRI measurements than

MRI is potentially the most appropriate technique for addressing all aspects of a complete cardiovascular disease examination. The evaluation of carotid flow will be a necessary capability in such examination. The methodologies presented in this chapter can improve the quality of the MRI measured data, reducing scan time, and improving the SNR. These

Boussel, L., Rayz, V., Martin, A., Acevedo-Bolton, G., Lawton, M. T., Higashida, R., Smith,

Canstein, C., Cachot, P., Faust, A., Stalder, A. F., Bock, J., Frydrychowicz, A., Kuffer, J., Hennig,

dynamics in identical vessel geometries, *Magn Reson Med* 59(3): 535–546. Carvalho, J. L. A., Nielsen, J. F. & Nayak, K. S. (2010). Feasibility of in vivo measurement of

Dai, G., Kaazempur-Mofrad, M. R., Natarajan, S., Zhang, Y., Vaughn, S., Blackman, B. R.,

Frayne, R. & Rutt, B. K. (1995). Measurement of fluid-shear rate by Fourier-encoded velocity

Glover, G. H. & Pelc, N. J. (1988). A rapid-gated cine MRI technique, *Magn Reson Annu*

Hahn, E. L. (1960). Detection of sea-water motion by nuclear precession, *J Geophys Res*

Harloff, A., Albrecht, F., Spreer, J., Stalder, A. F., Bock, J., Frydrychowicz, A., Schollhorn,

Harloff, A., Nussbaumer, A., Bauer, S., Stalder, A. F., Frydrychowicz, A., Weiller, C., Hennig,

Hoskins, P. R. (1996). Accuracy of maximum velocity estimates made using Doppler

aorta using flow-sensitive 4D MRI, *Magn Reson Med* 63(6): 1529–1536. He, X. & Ku, D. N. (1996). Pulsatile flow in the human left coronary artery bifurcation: average

J., Hetzel, A., Schumacher, M., Hennig, J. & Markl, M. (2009). 3D blood flow characteristics in the carotid artery bifurcation assessed by flow-sensitive 4D MRI

J. & Markl, M. (2010). In vivo assessment of wall shear stress in the atherosclerotic

*Sciences of the United States of America* 101(41): 14871–14876.

imaging, *Magn Reson Med* 34(3): 378–387.

at 3T, *Magn Reson Med* 61(1): 65–74.

conditions, *J Biomech Eng* 118(1): 74–82.

ultrasound systems, *Br J Radiol* 69(818): 172–177.

W. S., Young, W. L. & Saloner, D. (2009). 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*

J. & Markl, M. (2008). 3D MR flow analysis in realistic rapid-prototyping model systems of the thoracic aorta: comparison with in vivo data and computational fluid

carotid wall shear rate using spiral Fourier velocity encoded MRI, *Magn Reson Med*

Kamm, R. D., Garcia-Cardena, G. & Gimbrone Jr., M. A. (2004). Distinct endothelial phenotypes evoked by arterial waveforms derived from atherosclerosis-susceptible and -resistant regions of human vasculature, *Proceedings of the National Academy of*

techniques should improve the diagnosis and understanding of carotid diseases.

CFD alone.

**7. References**

61(2): 409–417.

63(6): 1537–1547.

pp. 299–333.

65(2): 776–777.


**0**

**24**

*Japan*

**Numerical Simulation for Intranasal**

Takahisa Yamamoto1, Seiichi Nakata2, Tsutomu Nakashima3

More than 10 million people in Japan suffer some of nasal diseases every year (Haruna (2003)); the paranasal sinusitis (so-called the empyema), hypertrophic rhinitis and inferior concha inflammation. Nebulizer treatment has been used for the nasal diseases. The effectiveness of the nebulizer treatment has been confirmed from clinical view points until now. However there are a few researches that evaluate the effect of the nebulizer treatment theoretically and quantitatively, i.e., the transport characteristics of medicinal droplets and their deposition on

The development of medical image processing technique is in progress and now gives us exquisitely detailed anatomic information. Some researchers calculated blood flow inside vital arteries and intranasal airflow characteristics by means of the medical image processing technique as well as Computed Fluid Dynamics (CFD) analysis. As for the CFD analysis of intranasal flow, Weinhold et al. constructed both a transparent resin model and numerical three-dimensional anatomy model with nasal cavities using a patient's CT data (Weinhold & Mlynski (2004)). They subsequently made clear airflow characteristics in the nose experimentally and numerically, and found that pressure drop was a main factor of nasal airflow. Lindemann et al. focused at a case which underwent radical sinus surgery (Lindemann et al. (2004; 2005)). In their case, both the lateral nasal wall and the turbinates, inhibiting physiological airflow, were removed by the surgery to realize the enlargements of the nasal cavity volumes and to increase the ratio between nasal cavity volume and surface area. However the researches mentioned here dealt with only a few patient-cases even though there are individual differences in the shape of human nasal cavity and in the grade of medical conditions. The past researches considered the individual differences insufficiently. The authors analyzed intranasal transport phenomena for several patient cases and compared each others in the past study (Monya et al. (2009); Yamamoto et al. (2009)). From the results characteristics of airflow and medicinal droplet transportation strongly depend on inflow conditions such as inflow angle, velocity and size of particle even if there are the individual

**1. Introduction**

the inflammation areas of nasal wall.

differences for the shape of patient's nasal cavity.

**Transport Phenomena**

and Tsuyoshi Yamamoto4

<sup>2</sup>*Fujita Health University* <sup>3</sup>*Nagoya University* <sup>4</sup>*Kyushu University*

<sup>1</sup>*Gifu National College of Technology*


## **Numerical Simulation for Intranasal Transport Phenomena**

Takahisa Yamamoto1, Seiichi Nakata2, Tsutomu Nakashima3 and Tsuyoshi Yamamoto4

> *Gifu National College of Technology Fujita Health University Nagoya University Kyushu University Japan*

#### **1. Introduction**

24 Will-be-set-by-IN-TECH

536 Fluid Dynamics, Computational Modeling and Applications

Thubrikar, M. J. & Robicsek, F. (1995). Pressure-induced arterial wall stress and

van Dijk, P. (1984). Direct cardiac NMR imaging of heart wall and blood flow velocity, *J Comput*

Winkler, A. J. & Wu, J. (1995). Correction of intrinsic spectral broadening errors in Doppler

Zarins, C. K., Giddens, D. P., Bharadvaj, B. K., Sottiurai, V. S., Mabon, R. F. & Glagov, S. (1983).

peak velocity measurements made with phased sector and linear array transducers,

Carotid bifurcation atherosclerosis. Quantitative correlation of plaque localization with flow velocity profiles and wall shear stress, *Circulation Research* 53(4): 502–514.

atherosclerosis, *Ann Thorac Surg* 59(6): 1594–1603.

*Assist Tomogr* 8(3): 429–436.

*Ultrasound Med Biol* 21(8): 1029–1035.

More than 10 million people in Japan suffer some of nasal diseases every year (Haruna (2003)); the paranasal sinusitis (so-called the empyema), hypertrophic rhinitis and inferior concha inflammation. Nebulizer treatment has been used for the nasal diseases. The effectiveness of the nebulizer treatment has been confirmed from clinical view points until now. However there are a few researches that evaluate the effect of the nebulizer treatment theoretically and quantitatively, i.e., the transport characteristics of medicinal droplets and their deposition on the inflammation areas of nasal wall.

The development of medical image processing technique is in progress and now gives us exquisitely detailed anatomic information. Some researchers calculated blood flow inside vital arteries and intranasal airflow characteristics by means of the medical image processing technique as well as Computed Fluid Dynamics (CFD) analysis. As for the CFD analysis of intranasal flow, Weinhold et al. constructed both a transparent resin model and numerical three-dimensional anatomy model with nasal cavities using a patient's CT data (Weinhold & Mlynski (2004)). They subsequently made clear airflow characteristics in the nose experimentally and numerically, and found that pressure drop was a main factor of nasal airflow. Lindemann et al. focused at a case which underwent radical sinus surgery (Lindemann et al. (2004; 2005)). In their case, both the lateral nasal wall and the turbinates, inhibiting physiological airflow, were removed by the surgery to realize the enlargements of the nasal cavity volumes and to increase the ratio between nasal cavity volume and surface area. However the researches mentioned here dealt with only a few patient-cases even though there are individual differences in the shape of human nasal cavity and in the grade of medical conditions. The past researches considered the individual differences insufficiently. The authors analyzed intranasal transport phenomena for several patient cases and compared each others in the past study (Monya et al. (2009); Yamamoto et al. (2009)). From the results characteristics of airflow and medicinal droplet transportation strongly depend on inflow conditions such as inflow angle, velocity and size of particle even if there are the individual differences for the shape of patient's nasal cavity.

Transport Phenomena 3

Numerical Simulation for Intranasal Transport Phenomena 539

Fig. 2. Front and side view of human head and nasal sinuses (adapted from Haruna (2003))

Fig. 3. Cross-sectional CT data for nasal cavity and nasal sinuses: lett) coronal and right)

Nasal septum locates on the middle of left and right nasal cavities, being almost straight stretches. Many people has a slightly-curved nasal septum as shown in the left-side figure of Fig.5. There are small size different between the left and right nostrils. The right figure of Fig.5 indicates a CT image that is the case of deviated nasal septum. The deviated nasal septum is caused by external nasal injuries and/or congenital origin. As for the cases, one of the nasal

sagittal sections

**2.2.1 Deviated nasal septum**

Fig. 1. Flowchart diagram of CFD analysis for heat and mass transport phenomena in nasal cavity

Figure 1 indicates the flowchart of the CFD analysis for intranasal transport phenomena. This chapter shows you detail features of nasal cavity and some nasal diseases, subsequently how to yield three-dimensional model of nasal cavity and calculation mesh for CFD analysis. Finally the characteristics of intranasal transport phenomena are presented in this chapter.

#### **2. Anatomy of nasal cavity**

#### **2.1 Structure of nasal cavity and nasal sinuses**

Figure 2 shows front and side views of human head, nasal cavity, nasal conchae and nasal sinuses. The nasal cavity is a large air filled space, and the pathway of respiration flow. The nasal cavity conditions the air to be received by the other areas of the respiratory tract. The nasal conchae, which has large surface area, warms and cools passing air through the nasal cavity. At the same time, the passing air is humidified and cleaned by nasal conchae and nasal hairs; most parts of dusts and particulate matters are removed there. The nasal sinuses are four caviums existing in the ossa faciei as shown in Fig.2; maxillary sinus, frontal sinus, ethmoidal sinus and sphenoidal sinus. These sinuses connect to the nasal cavity via very small ducts, ostia; both diameter and streamwise scales are several millimeters. Figures 3 represents coronal and sagittal cross-sectional CT data for nasal cavity and nasal sinuses. The white areas in the figure represent bones of the patient's head, and the gray areas indicate muscles and fat. The nasal cavity and sinuses are represented as dark area in the CT data.

In paranasal sinusitis, both the sinuses and the ostia are infected by bacillus and the virus. The details of such the diseases are stated in following subsections.

#### **2.2 Nasal diseases**

The nose often suffers some injuries such as fractures because the nose is protruded forward from the human face. Infection, epistaxis, as well as polyps are also found in the nasal diseases. Almost all people have contract rhinitis which is caused by the inflammation of nasal mucosa. Sinusitis will develop when the inflammation of the rhinitis spread to the sinuses via ostia.

2 Will-be-set-by-IN-TECH

Computer Tomography (CT) data of nasal cavity; five cases

Generation of 3-dimensional geometry model ・ Marching-Cubes method

CFD analysis grid generation of nasal cavity

・ Delaunay triangulation ・ Unstructured grid

・ CFD analysis of gas-liquid two phase flow

・ CFD analysis of gas flow in nasal cavities

Fig. 1. Flowchart diagram of CFD analysis for heat and mass transport phenomena in nasal

Figure 1 indicates the flowchart of the CFD analysis for intranasal transport phenomena. This chapter shows you detail features of nasal cavity and some nasal diseases, subsequently how to yield three-dimensional model of nasal cavity and calculation mesh for CFD analysis. Finally the characteristics of intranasal transport phenomena are presented in this chapter.

Figure 2 shows front and side views of human head, nasal cavity, nasal conchae and nasal sinuses. The nasal cavity is a large air filled space, and the pathway of respiration flow. The nasal cavity conditions the air to be received by the other areas of the respiratory tract. The nasal conchae, which has large surface area, warms and cools passing air through the nasal cavity. At the same time, the passing air is humidified and cleaned by nasal conchae and nasal hairs; most parts of dusts and particulate matters are removed there. The nasal sinuses are four caviums existing in the ossa faciei as shown in Fig.2; maxillary sinus, frontal sinus, ethmoidal sinus and sphenoidal sinus. These sinuses connect to the nasal cavity via very small ducts, ostia; both diameter and streamwise scales are several millimeters. Figures 3 represents coronal and sagittal cross-sectional CT data for nasal cavity and nasal sinuses. The white areas in the figure represent bones of the patient's head, and the gray areas indicate muscles and

In paranasal sinusitis, both the sinuses and the ostia are infected by bacillus and the virus. The

The nose often suffers some injuries such as fractures because the nose is protruded forward from the human face. Infection, epistaxis, as well as polyps are also found in the nasal diseases. Almost all people have contract rhinitis which is caused by the inflammation of nasal mucosa. Sinusitis will develop when the inflammation of the rhinitis spread to the sinuses via

fat. The nasal cavity and sinuses are represented as dark area in the CT data.

details of such the diseases are stated in following subsections.

CFD analysis

in nasal cavities

cavity

**2. Anatomy of nasal cavity**

**2.2 Nasal diseases**

ostia.

**2.1 Structure of nasal cavity and nasal sinuses**

Fig. 2. Front and side view of human head and nasal sinuses (adapted from Haruna (2003))

Fig. 3. Cross-sectional CT data for nasal cavity and nasal sinuses: lett) coronal and right) sagittal sections

#### **2.2.1 Deviated nasal septum**

Nasal septum locates on the middle of left and right nasal cavities, being almost straight stretches. Many people has a slightly-curved nasal septum as shown in the left-side figure of Fig.5. There are small size different between the left and right nostrils. The right figure of Fig.5 indicates a CT image that is the case of deviated nasal septum. The deviated nasal septum is caused by external nasal injuries and/or congenital origin. As for the cases, one of the nasal

Transport Phenomena 5

Numerical Simulation for Intranasal Transport Phenomena 541

Figure 6 shows cross-sectional CT data for both healthy maxillary sinus and sinusitis one. The sinusitis happens anywhere about four kinds of the sinus paranasalis such as the maxillary sinus, sinus ethmoidales, the master wrestlers caves, and the sphenoidal sinuses. The sinusitis is classified into the acute sinusitis (in a short term) and the chronic sinusitis (in a long term).

Fig. 6. Cross-sectional CT data for maxillary sinuses: left) healthy maxillary sinus and right)

The acute sinusitis is caused by a variety of bacteria and virus, often develops after the blockage of ostia which are opening of sinuses connecting to nasal cavity. The blockage normally occures as a result of virus infection in the upper respiratory caused by a cold. The cold brings on the swelling of the mucous membrane of nostril, and then the ostia are easily obstructed. In the blocked sinus, the air inside the cavity are absorbed to the bloodstream, and subsequently the pressure inside the sinus decreases. This pressure drop produce sever pains and the sinus fills with secretory fluid. The secretory fluid becomes a breeding ground for viruses. In order to attack the viruses, white blood cell and other matters such as the secretory fluid are aggregated in the blocked sinus, then the pressure in the sinus is increased and pain

The case that a symptom of the sinusitis continues more than 8-12 weeks is called the chronic sinusitis. The developing mechanism of the chronic sinusitis has not been clear, however it has been confirmed that the chronic sinusitis develops after the infection of the virus, severe allergies and the influences of the environmental pollution material. A genetic contributor is also regarded as one of the factors concerned with development of the chronic sinusitis.

In order to construct a three-dimensional model for nasal cavity, data conversion algorithm from two-dimensional CT/MRI data to three-dimensional geometric model as well as smoothing algorithm for the three-dimensional model are required. Furthermore, mesh generation models are needed to execute CFD analysis of intranasal transport phenomena.

Marching cubes algorithm is one of the latest models of surface construction used for viewing three-dimensional data (Lorensen & Cline (1987)). This algorithm extracts a polygonal mesh

**3. Construction of three-dimensional nasal cavity model**

This section describes the fundamental theories of these algorithms.

**3.1 Data conversion algorithm from 2-D to 3-D: marching cubes algorithm**

**2.2.3 Sinusitis**

sinusitis

becomes more sever.

Fig. 4. Series of axial cross-sectional CT data for nasal cavity and nasal sinuses (adapted from Haruna (2003))

Fig. 5. Axial cross-sectional CT data for nasal septum: left) healthy nasal septum and right) deviated nasal septum

cavity is significantly smaller than the other. Then the patient suffers nasal congestion and drying nasal cavity, consequently epistaxis. In the sever cases operative treatment is required.

#### **2.2.2 Rhinitis (nasal inflammation)**

The nasal cavity is the part that easily causes infection in the upper respiratory tracts. Rhinitis (nasal inflammation) is mainly classified into two cases; the one is acute rhinitis and the other is chronic rhinitis. The acute rhinitis is generally caused by the infection of virus and allergies of pollen and house dust. The difference point between both rhinitis is the term of these diseases; the former occurs in a short term and the latter in a long term (two weeks and more). In some cases, the chronic rhinitis occasionally accompanies the chronic sinusitis.

#### **2.2.3 Sinusitis**

4 Will-be-set-by-IN-TECH

Fig. 4. Series of axial cross-sectional CT data for nasal cavity and nasal sinuses (adapted from Haruna

Fig. 5. Axial cross-sectional CT data for nasal septum: left) healthy nasal septum and right)

cavity is significantly smaller than the other. Then the patient suffers nasal congestion and drying nasal cavity, consequently epistaxis. In the sever cases operative treatment is required.

The nasal cavity is the part that easily causes infection in the upper respiratory tracts. Rhinitis (nasal inflammation) is mainly classified into two cases; the one is acute rhinitis and the other is chronic rhinitis. The acute rhinitis is generally caused by the infection of virus and allergies of pollen and house dust. The difference point between both rhinitis is the term of these diseases; the former occurs in a short term and the latter in a long term (two weeks and more).

In some cases, the chronic rhinitis occasionally accompanies the chronic sinusitis.

(2003))

deviated nasal septum

**2.2.2 Rhinitis (nasal inflammation)**

Figure 6 shows cross-sectional CT data for both healthy maxillary sinus and sinusitis one. The sinusitis happens anywhere about four kinds of the sinus paranasalis such as the maxillary sinus, sinus ethmoidales, the master wrestlers caves, and the sphenoidal sinuses. The sinusitis is classified into the acute sinusitis (in a short term) and the chronic sinusitis (in a long term).

Fig. 6. Cross-sectional CT data for maxillary sinuses: left) healthy maxillary sinus and right) sinusitis

The acute sinusitis is caused by a variety of bacteria and virus, often develops after the blockage of ostia which are opening of sinuses connecting to nasal cavity. The blockage normally occures as a result of virus infection in the upper respiratory caused by a cold. The cold brings on the swelling of the mucous membrane of nostril, and then the ostia are easily obstructed. In the blocked sinus, the air inside the cavity are absorbed to the bloodstream, and subsequently the pressure inside the sinus decreases. This pressure drop produce sever pains and the sinus fills with secretory fluid. The secretory fluid becomes a breeding ground for viruses. In order to attack the viruses, white blood cell and other matters such as the secretory fluid are aggregated in the blocked sinus, then the pressure in the sinus is increased and pain becomes more sever.

The case that a symptom of the sinusitis continues more than 8-12 weeks is called the chronic sinusitis. The developing mechanism of the chronic sinusitis has not been clear, however it has been confirmed that the chronic sinusitis develops after the infection of the virus, severe allergies and the influences of the environmental pollution material. A genetic contributor is also regarded as one of the factors concerned with development of the chronic sinusitis.

#### **3. Construction of three-dimensional nasal cavity model**

In order to construct a three-dimensional model for nasal cavity, data conversion algorithm from two-dimensional CT/MRI data to three-dimensional geometric model as well as smoothing algorithm for the three-dimensional model are required. Furthermore, mesh generation models are needed to execute CFD analysis of intranasal transport phenomena. This section describes the fundamental theories of these algorithms.

#### **3.1 Data conversion algorithm from 2-D to 3-D: marching cubes algorithm**

Marching cubes algorithm is one of the latest models of surface construction used for viewing three-dimensional data (Lorensen & Cline (1987)). This algorithm extracts a polygonal mesh

Transport Phenomena 7

Numerical Simulation for Intranasal Transport Phenomena 543

In the image processing and image transmission, noises are physically captured and fed into the processes due to various factors. For instance, noises are generated caused by the image resolution and image slice thickness of CT and MRI in the field of medical engineering. When the noise generation mechanism is clear and mathematically modeled, we can remove the noise using the optimized filter. However, it is difficult to modeled the noise mathematically in almost all situations.f The smoothing processes are needed to reduce and minimize the influence of the noise. The noise makes patterned indented surfaces in the three-dimensional model, therefore affects the quality of the constructed nasal cavity model significantly. There are two smoothing algorithms that often uses in the medical image processing; moving-average algorithm and median-average algorithm. This subsection

**3.2 Smoothing algorithm**

presents both algorithms as follows.

**3.2.1 Moving-average algorithm**

1 1 1

1 20 1

1 1 1

1+1+1+1+20+1+1+1+1

1 1 1 1 3.1 1

Fig. 9. Procedure of Moving-Avarage Algorithm for 3 × 3 Cubes

*<sup>g</sup>*(*i*, *<sup>j</sup>*) = <sup>1</sup>

consequently smoothes the object surfaces.

**3.2.2 Median-average algorithm**

*n*2

1 1 1

(a) noise (b) asperity

Figure 9 shows the application of the moving-average algorithm for 3 × 3 cubes (Savitzky & Golay (1964); Sun (2006); Tamura (2002)). The number inside a cube expresses a scalar variable, for example the concentration of a chemical specie. The case of Fig.9(a) presents a noise included in the data and that of Fig.9(b) expresses an asperity in the data. The moving-average algorithm calculates an average of surroundings of the concentration *f*(*i*, *j*)

> *n*/2 ∑ *l*=−[*l*/2]

The moving-average algorithm feathers the edge and the boundary of the input image,

Figure 10 indicates the application of the median-average algorithm for 3 × 3 cubes (Boyle & Thomas (1988); Tamura (2002)). The cases of Fig.10(a) and (b) represent a noise and an asperity in the data, respectively. This algorithm sorts both a concentration value *f*(*i*, *j*) at a target position of the input image and its surroundings in ascending order, subsequently converts the output image *g*(*i*, *j*) using their median value. As for Fig.10(a), the concentration values in input image are rearranged like "1, 1, 1, 1, 1, 1, 1, 1, 20" in ascending order, and then the fifth value "1" is picked up as a median value and used as the concentration *g*(*i*, *j*) in the

in the input image, and converts to the concentration *g*(*i*, *j*) in the output image.

*n*/2 ∑ *k*=−[*n*/2] 1 1 20

1 20 20

<sup>9</sup> =3.1 1+1+20+1+20+20+1+1+20

1 20 1

1 1 20

1 9.4 20

<sup>9</sup> =9.4

*f*(*i* + *k*, *j* + *l*) (1)

1 20

1

Fig. 7. Configurations of fourteen unique cubes for marching cubes algorithm

Fig. 8. Triangular surface patterns for the marching cubes algorithm

of isosurface from a three-dimensional scalar field. The algorithm proceeds through the scalar field, taking eight neighbor locations at a time, then determining the polygons needed to represent the part of the isosurface that passes through this cube. The individual polygons are then fused into the desired surface as shown in Fig.7. This is done by creating an index to a pre-calculated array of 256 possible polygon configurations (28 = 256) within the cube, by treating each of the 8 scalar values as a bit in an 8-bit integer. Finally each vertex of the generated polygons is placed on the appropriate position along the cube's edge by linearly interpolating the two scalar values that are connected by that edge. The precalculated array of 256 cube configurations can be obtained by reflections and symmetrical rotations of 14 unique cases as shown in Fig.8.

The applications of this algorithm are mainly concerned with medical visualizations such as CT and MRI scan data images, and special effects on three-dimensional modeling.

#### **3.2 Smoothing algorithm**

6 Will-be-set-by-IN-TECH

**Triangular surface**

Fig. 7. Configurations of fourteen unique cubes for marching cubes algorithm

Fig. 8. Triangular surface patterns for the marching cubes algorithm

cases as shown in Fig.8.

of isosurface from a three-dimensional scalar field. The algorithm proceeds through the scalar field, taking eight neighbor locations at a time, then determining the polygons needed to represent the part of the isosurface that passes through this cube. The individual polygons are then fused into the desired surface as shown in Fig.7. This is done by creating an index to a pre-calculated array of 256 possible polygon configurations (28 = 256) within the cube, by treating each of the 8 scalar values as a bit in an 8-bit integer. Finally each vertex of the generated polygons is placed on the appropriate position along the cube's edge by linearly interpolating the two scalar values that are connected by that edge. The precalculated array of 256 cube configurations can be obtained by reflections and symmetrical rotations of 14 unique

The applications of this algorithm are mainly concerned with medical visualizations such as

CT and MRI scan data images, and special effects on three-dimensional modeling.

**Element point**

**Inside point Outside point Interpolation point** In the image processing and image transmission, noises are physically captured and fed into the processes due to various factors. For instance, noises are generated caused by the image resolution and image slice thickness of CT and MRI in the field of medical engineering. When the noise generation mechanism is clear and mathematically modeled, we can remove the noise using the optimized filter. However, it is difficult to modeled the noise mathematically in almost all situations.f The smoothing processes are needed to reduce and minimize the influence of the noise. The noise makes patterned indented surfaces in the three-dimensional model, therefore affects the quality of the constructed nasal cavity model significantly. There are two smoothing algorithms that often uses in the medical image processing; moving-average algorithm and median-average algorithm. This subsection presents both algorithms as follows.

#### **3.2.1 Moving-average algorithm**

Fig. 9. Procedure of Moving-Avarage Algorithm for 3 × 3 Cubes

Figure 9 shows the application of the moving-average algorithm for 3 × 3 cubes (Savitzky & Golay (1964); Sun (2006); Tamura (2002)). The number inside a cube expresses a scalar variable, for example the concentration of a chemical specie. The case of Fig.9(a) presents a noise included in the data and that of Fig.9(b) expresses an asperity in the data. The moving-average algorithm calculates an average of surroundings of the concentration *f*(*i*, *j*) in the input image, and converts to the concentration *g*(*i*, *j*) in the output image.

$$g(i,j) = \frac{1}{n^2} \sum\_{k=-\lceil n/2 \rceil}^{n/2} \sum\_{l=-\lceil l/2 \rceil}^{n/2} f(i+k, j+l) \tag{1}$$

The moving-average algorithm feathers the edge and the boundary of the input image, consequently smoothes the object surfaces.

#### **3.2.2 Median-average algorithm**

Figure 10 indicates the application of the median-average algorithm for 3 × 3 cubes (Boyle & Thomas (1988); Tamura (2002)). The cases of Fig.10(a) and (b) represent a noise and an asperity in the data, respectively. This algorithm sorts both a concentration value *f*(*i*, *j*) at a target position of the input image and its surroundings in ascending order, subsequently converts the output image *g*(*i*, *j*) using their median value. As for Fig.10(a), the concentration values in input image are rearranged like "1, 1, 1, 1, 1, 1, 1, 1, 20" in ascending order, and then the fifth value "1" is picked up as a median value and used as the concentration *g*(*i*, *j*) in the

Transport Phenomena 9

Numerical Simulation for Intranasal Transport Phenomena 545

to execute CFD analysis for such the transport phenomena, we have to consider multi-phase

The Reynolds number of intranasal flow is a few hundred. Therefore the flow is assumed as

here *ρ<sup>F</sup>* and **uF** are density and velocity of continuous phase. Equation 2 is the general form of the mass conservation equation and is valid for incompressible as well as compressible flows.

where *p* is the static pressure, *μ<sup>F</sup>* is viscosity of the fluid, and *ρ***g** and **F** are gravitational body force and external body forces, respectively. **F** contains interaction with the dispersed phase

The thermal energy transport equation is meant to be used for flows which are low speed and

where *cp* and *λ<sup>F</sup>* are specific heat at constant pressure and heat conductivity, respectively. *SFP* denotes a source term regarding the interaction between continuous phase and disperse

Consider a discrete particle traveling in a continuous fluid medium. The forces acting on the particle which affect the particle acceleration are due to the difference in velocity between the particle and fluid, as well as to the displacement of the fluid by the particle. The equation of motion for such a particle was derived by Basset, Boussinesq and Oseen for a rotating

**Basset** <sup>+</sup> **<sup>F</sup><sup>i</sup>**

where, superscript *i* denote the *i*-th particle, and **FDrag**, **FPressure**, **FBasset**, **FBuoyancy**, **FRotation** express drag force acting on the particle, pressure gradient force, Basset force, buoyancy force due to gravity and forces due to domain rotation, respectively. The pressure gradient force term applies on the particle due to the pressure gradient in the fluid surrounding the particle caused by fluid acceleration. The Basset force term accounts for the deviation in flow pattern

The aerodynamic drag force on a particle is propotional to the slip velocity between the

)2*ρFCD***uF** <sup>−</sup> **uP**(**uF** <sup>−</sup> **<sup>u</sup><sup>i</sup>**

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> ∇ · *<sup>ρ</sup>F***uF** <sup>=</sup> <sup>0</sup> (2)

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> **uF** · ∇*ρF***uF** <sup>=</sup> −∇*<sup>p</sup>* <sup>+</sup> *<sup>μ</sup>F*∇<sup>2</sup>**uF** <sup>+</sup> *<sup>ρ</sup>F***<sup>g</sup>** <sup>+</sup> **<sup>F</sup>** (3)

*<sup>∂</sup><sup>t</sup>* <sup>+</sup> **uF** · ∇*ρFcpTF* <sup>=</sup> ∇ · (*λF*∇*TF*) + **SFP** (4)

**Buoyancy** <sup>+</sup> **<sup>F</sup><sup>i</sup>**

**Rotation** <sup>+</sup> **<sup>F</sup><sup>i</sup>**

**Others** (5)

**<sup>P</sup>**) (6)

The equation for conservation of mass, or continuity equation, can be written as follows:

flow. This section describes the basic theory of CFD analysis for two-phase flow.

*∂ρ<sup>F</sup>*

Conservation of momentum in an inertial reference frame is described by

in multi-phase flow and other model-dependent source terms.

*∂ρFcpTF*

**Drag** <sup>+</sup> **<sup>F</sup><sup>i</sup>**

**Fi**

**Drag** <sup>=</sup> <sup>1</sup> 8 *π*(*d<sup>i</sup>*

**Pressure** <sup>+</sup> **<sup>F</sup><sup>i</sup>**

**4.2 Governing equations for disperse phase**

reference frame (Akiyama (2002)):

*mi p d***up** *dt* <sup>=</sup> **<sup>F</sup><sup>i</sup>**

from a steady state.

particle and the fluid velocity

**4.1 Governing equations for continuous phase**

*∂ρF***uF**

incompressible and laminar flow.

close to constant density.

phase.

Fig. 10. Procedure of Median-Average algorithm for 3 × 3 cubes

Fig. 11. Schematic drawing of the Delaunay triangulation model

output image. On the other hand, the concentration values in input image are rearranged like "1, 1, 1, 1, 1, 20, 20, 20, 20" in ascending order as shown in Fig.10(b). Finally, the fifth value "1" is selected as the concentration *g*(*i*, *j*) in the output image. The median-average algorithm can remove the noise effectively compared with the moving-average algorithm. At the same time, this algorithm involves a risk that the magnitude of rearranging the original CT/MRI data is larger than the moving-average algorithm, consequently could yield different three-dimensional model from the real target.

#### **3.3 Mesh generation: Delaunay triangulation model**

The Delaunay triangulation model is one of the most-used mesh generation technique. This model is applicable for a wide variety of very complicate geometries.

This model mathematically ensures that the circumcircle associated with each triangle contains no other point in its interior. This definition extends naturally to three-dimensional problems. Figure 11 shows a schematic drawing of the Delaunay triangulation model (Nakahashi (1995)). The triangles as shown in Fig.11(a) contains the edges and the edge's endpoints. A new point "P" is added inside the triangle "ABC", then the mathematical theory of the Delaunay triangulation is not satisfied in the triangle "ABC" and its surrounding triangle "BCD". The model can eliminate the edge "BC" which is a common edge for both triangles "ABC" and "BCD", subsequently reconstruct triangles using the added point "P". Consequently, as shown in Fig.11(d), the model can yield four triangles "ABP", "ACP", "BDP" and "CDP". The Delaunay triangulation constructs and reconstructs triangles in target domain simultaneously and give us fine meshes for CFD analysis.

#### **4. Theory of computational fluid dynamics**

In the nasal cavity, particulate matters such as pollens, house dust as well as medicinal droplets are transported via airflow of the respiration and the nebulizer treatment. In order 8 Will-be-set-by-IN-TECH

1, 1, 1, 1, 1, 1, 1, 1, 20 1, 1, 1, 1, 1, 20, 20, 20, 20

(a) (b)

<sup>D</sup> <sup>A</sup>

(a) Plane triangulation (b) Insertion P-point (c) Delitation BC-element (d) Plane triangulation

output image. On the other hand, the concentration values in input image are rearranged like "1, 1, 1, 1, 1, 20, 20, 20, 20" in ascending order as shown in Fig.10(b). Finally, the fifth value "1" is selected as the concentration *g*(*i*, *j*) in the output image. The median-average algorithm can remove the noise effectively compared with the moving-average algorithm. At the same time, this algorithm involves a risk that the magnitude of rearranging the original CT/MRI data is larger than the moving-average algorithm, consequently could yield different

The Delaunay triangulation model is one of the most-used mesh generation technique. This

This model mathematically ensures that the circumcircle associated with each triangle contains no other point in its interior. This definition extends naturally to three-dimensional problems. Figure 11 shows a schematic drawing of the Delaunay triangulation model (Nakahashi (1995)). The triangles as shown in Fig.11(a) contains the edges and the edge's endpoints. A new point "P" is added inside the triangle "ABC", then the mathematical theory of the Delaunay triangulation is not satisfied in the triangle "ABC" and its surrounding triangle "BCD". The model can eliminate the edge "BC" which is a common edge for both triangles "ABC" and "BCD", subsequently reconstruct triangles using the added point "P". Consequently, as shown in Fig.11(d), the model can yield four triangles "ABP", "ACP", "BDP" and "CDP". The Delaunay triangulation constructs and reconstructs triangles in target domain

In the nasal cavity, particulate matters such as pollens, house dust as well as medicinal droplets are transported via airflow of the respiration and the nebulizer treatment. In order

1 1 20

P P P

C

B

1 20 20

1 20 1

<sup>D</sup> <sup>A</sup>

1 1 20

B

D

C

1 1 20

1 20

1 1 1

1 1 1

Fig. 10. Procedure of Median-Average algorithm for 3 × 3 cubes

B

C

A

Fig. 11. Schematic drawing of the Delaunay triangulation model

model is applicable for a wide variety of very complicate geometries.

1 1 1

1

A

1 1 1

B

D

three-dimensional model from the real target.

**3.3 Mesh generation: Delaunay triangulation model**

simultaneously and give us fine meshes for CFD analysis.

**4. Theory of computational fluid dynamics**

C

1 20 1

1 1 1

to execute CFD analysis for such the transport phenomena, we have to consider multi-phase flow. This section describes the basic theory of CFD analysis for two-phase flow.

#### **4.1 Governing equations for continuous phase**

The Reynolds number of intranasal flow is a few hundred. Therefore the flow is assumed as incompressible and laminar flow.

The equation for conservation of mass, or continuity equation, can be written as follows:

$$\frac{\partial \rho\_F}{\partial t} + \nabla \cdot \rho\_F \mathbf{u}\_F = 0 \tag{2}$$

here *ρ<sup>F</sup>* and **uF** are density and velocity of continuous phase. Equation 2 is the general form of the mass conservation equation and is valid for incompressible as well as compressible flows. Conservation of momentum in an inertial reference frame is described by

$$\frac{\partial \rho\_F \mathbf{u\_F}}{\partial t} + \mathbf{u\_F} \cdot \nabla \rho\_F \mathbf{u\_F} = -\nabla p + \mu\_F \nabla^2 \mathbf{u\_F} + \rho\_F \mathbf{g} + \mathbf{F} \tag{3}$$

where *p* is the static pressure, *μ<sup>F</sup>* is viscosity of the fluid, and *ρ***g** and **F** are gravitational body force and external body forces, respectively. **F** contains interaction with the dispersed phase in multi-phase flow and other model-dependent source terms.

The thermal energy transport equation is meant to be used for flows which are low speed and close to constant density.

$$\frac{\partial \rho\_F c\_p T\_F}{\partial t} + \mathbf{u\_F} \cdot \nabla \rho\_F c\_p T\_F = \nabla \cdot (\lambda\_F \nabla T\_F) + \mathbf{S\_{FP}} \tag{4}$$

where *cp* and *λ<sup>F</sup>* are specific heat at constant pressure and heat conductivity, respectively. *SFP* denotes a source term regarding the interaction between continuous phase and disperse phase.

#### **4.2 Governing equations for disperse phase**

Consider a discrete particle traveling in a continuous fluid medium. The forces acting on the particle which affect the particle acceleration are due to the difference in velocity between the particle and fluid, as well as to the displacement of the fluid by the particle. The equation of motion for such a particle was derived by Basset, Boussinesq and Oseen for a rotating reference frame (Akiyama (2002)):

$$m\_p^i \frac{d\mathbf{u\_P}}{dt} = \mathbf{F\_{Drag}^i} + \mathbf{F\_{Pressure}^i} + \mathbf{F\_{Baseet}^i} + \mathbf{F\_{Buoyancy}^i} + \mathbf{F\_{Rotation}^i} + \mathbf{F\_{Others}^i} \tag{5}$$

where, superscript *i* denote the *i*-th particle, and **FDrag**, **FPressure**, **FBasset**, **FBuoyancy**, **FRotation** express drag force acting on the particle, pressure gradient force, Basset force, buoyancy force due to gravity and forces due to domain rotation, respectively. The pressure gradient force term applies on the particle due to the pressure gradient in the fluid surrounding the particle caused by fluid acceleration. The Basset force term accounts for the deviation in flow pattern from a steady state.

The aerodynamic drag force on a particle is propotional to the slip velocity between the particle and the fluid velocity

$$\mathbf{F\_{Drag}^{i}} = \frac{1}{8}\pi (d^i)^2 \rho\_F C\_D \mathbf{u\_F} - \mathbf{u\_P}(\mathbf{u\_F} - \mathbf{u\_P^i}) \tag{6}$$

Transport Phenomena 11

Numerical Simulation for Intranasal Transport Phenomena 547

Chronic sinusitis hypertrophic

3 � Left & Right Left

5 � Left & Right Left

can neglect both the latent heat transfer and the radiative heat transfer. The convective heat

*<sup>P</sup>* is the temperature of the *i*-th particle, and Nu is the Nusselt number given by

The convective heat transfer *QC* is calculated and used in the source term of the thermal

Three-dimensional models of nasal cavity were constructed by means of five actual patient's CT data. These cases suffer several grades of deviated nasal septum symptom, hypertrophic rhinitis, chronic sinusitis and inferior nasal concha swelling. The case data are summarized in Table 1. Since head CT data involves not only nasal cavity but also skin, fat, bone and purulent matter, it is required to eliminate them. This study used three-dimensional medical image processing software "Mimics (Materialize Inc.)" and the marching cubes algorithm to determine the surfaces of nasal wall. The moving-average algorithm is also applied to smoothing the surface. Figure 12 shows the three-dimensional nasal model of case-I constructed in this study. The anterior and the posterior of the figure show anterior nares and upper pharynx, respectively. The three-dimension model has curved nasal septum which corresponds to the curvature as shown in Fig.12. The cross-sectional areas of nasal cavities

*Nu* = 2 + 0.6*Re*0.5(*μ<sup>F</sup>*

*<sup>λ</sup>FNu*(*TF* <sup>−</sup> *<sup>T</sup><sup>i</sup>*

*Cp λF* ) 1

rhinitis

Bloating inferior

concha

*<sup>P</sup>*) (15)

<sup>3</sup> (16)

Fig. 12. Three-dimensional nasal cavity model (case-1)

2 Left & Right

1 � Left Left

4 � Left & Right

*Qi <sup>C</sup>* <sup>=</sup> *<sup>π</sup>d<sup>i</sup>*

**5. Numerical simulations of intranasal transport phenomena**

Case No. Deviated nasal septum

Table 1. Detail patient's case data

energy transport equation, Eq.4.

**5.1 Case data and three-dimensional model**

transfer *QC* is given by

where *T<sup>i</sup>*

where **uF** and **uP** are the velocities of the fluid flow and the particle. *CD* denotes the drag coefficient of the particle.

$$\mathcal{C}\_D = \frac{24}{Re^i (1 + 0.15 (Re^i)^{0.687})} \tag{7}$$

*Re<sup>i</sup>* is the Reynolds number based on the slip velocity between the fluid flow and the *i*-th particle velocities.

$$Re^{i} = \frac{d^{i}(\mathbf{u\_{F}} - \mathbf{u\_{P}^{i}})}{\nu\_{F}} \tag{8}$$

The pressure gradient force results from the local fluid pressure gradient around the particle and is defined as

$$\mathbf{F}\_{\text{pressure}}^{\text{i}} = -\frac{1}{6}\pi (d^{i})^{3} \nabla p \tag{9}$$

This force is only important if large fluid pressure gradients exist and if the particle density is smaller than or similar to the fluid density.

$$\mathbf{F}\_{\mathbf{Based}}^{\mathbf{i}} = \frac{3}{2} (d^{i})^{2} \sqrt{\pi \rho\_{F} \mu\_{F}} \int\_{t\_{0}}^{t} \frac{d(\mathbf{u}\_{\mathbf{F}} - \mathbf{u}\_{\mathbf{p}}^{\mathbf{i}})}{d\tau} (t - \tau)^{-\frac{1}{2}} d\tau \tag{10}$$

The buoyancy force is the force on a particle immersed in a fluid. The buoyancy force is equal to the weight of the displaced fluid and is given by

$$\mathbf{F}\_{\mathbf{Buoyancy}}^{\mathbf{i}} = \frac{3}{2} (d^i)^3 (\rho\_P^i - \rho\_F) \mathbf{g} \tag{11}$$

here **g** expresses the gravity vector.

In rotating frame of reference, the rotation term is an intrinsic part of the acceleration in and is the sum ob Coriolis and centripetal forces.

$$\mathbf{F\_{rotation}} = -\frac{1}{3} (d^i)^3 \rho\_P (\boldsymbol{\omega} \times \mathbf{u\_P}) - \frac{1}{6} \rho\_P \boldsymbol{\omega} \times (\boldsymbol{\omega} \times \mathbf{r\_P}) \tag{12}$$

where, *ω* denotes rotation speed of the reference.

In the intranasal flow and heat and mass transfer, the influence of both the pressure gradient term and the Basset force term on the transport phenomena in the nasal cavity is extremely small. The momentum conservation equation for the particle is shown as

$$m\_{\mathbf{P}}^i \frac{d\mathbf{u}\_{\mathbf{P}}^i}{dt} = \frac{1}{8} \pi (d^i)^2 \rho\_F \mathbb{C}\_D \mathbf{u}\_{\mathbf{F}} - \mathbf{u}\_{\mathbf{P}}^i (\mathbf{u}\_{\mathbf{F}} - \mathbf{u}\_{\mathbf{P}}^i) \tag{13}$$

When the Lagrangian method is applied to the particle movement, the location of the particle is given by

$$\frac{d\mathbf{x\_P^i}}{dt} = \mathbf{u\_P^i} \tag{14}$$

As for the almost all particulate matter transported in nasal cavity, the mass fraction of the disperse phase are extremely small compared with the continuous phase. Therefore it is only necessary to consider the one-way interaction between both phases, namely only the continuous phase influences the disperse phase.

Heat transfer concerning the disperse phase is governed by three physical processes, namely convective heat transfer, latent heat transfer associated with mass transfer, and radiative heat transfer. Since the heat transfer of intranasal flow is not so high temperature condition, we 10 Will-be-set-by-IN-TECH

where **uF** and **uP** are the velocities of the fluid flow and the particle. *CD* denotes the drag

*Re<sup>i</sup>* is the Reynolds number based on the slip velocity between the fluid flow and the *i*-th

The pressure gradient force results from the local fluid pressure gradient around the particle

This force is only important if large fluid pressure gradients exist and if the particle density is

 *t t*0

The buoyancy force is the force on a particle immersed in a fluid. The buoyancy force is equal

2 (*di* )3(*ρ<sup>i</sup>*

In rotating frame of reference, the rotation term is an intrinsic part of the acceleration in and

)3*ρP*(*<sup>ω</sup>* <sup>×</sup> **uP**) <sup>−</sup> <sup>1</sup>

In the intranasal flow and heat and mass transfer, the influence of both the pressure gradient term and the Basset force term on the transport phenomena in the nasal cavity is extremely

When the Lagrangian method is applied to the particle movement, the location of the particle

As for the almost all particulate matter transported in nasal cavity, the mass fraction of the disperse phase are extremely small compared with the continuous phase. Therefore it is only necessary to consider the one-way interaction between both phases, namely only the

Heat transfer concerning the disperse phase is governed by three physical processes, namely convective heat transfer, latent heat transfer associated with mass transfer, and radiative heat transfer. Since the heat transfer of intranasal flow is not so high temperature condition, we

*d***x<sup>i</sup> P** *dt* <sup>=</sup> **<sup>u</sup><sup>i</sup>**

)2*ρFCD***uF** <sup>−</sup> **<sup>u</sup><sup>i</sup>**

*<sup>d</sup>*(**uF** <sup>−</sup> **<sup>u</sup><sup>i</sup>**

6

**<sup>P</sup>**(**uF** <sup>−</sup> **<sup>u</sup><sup>i</sup>**

**p**) *<sup>d</sup><sup>τ</sup>* (*<sup>t</sup>* <sup>−</sup> *<sup>τ</sup>*)<sup>−</sup> <sup>1</sup>

(**uF** <sup>−</sup> **<sup>u</sup><sup>i</sup> P**)

*νF*

6 *π*(*d<sup>i</sup>*

*Rei*(<sup>1</sup> <sup>+</sup> 0.15(*Rei*)0.687) (7)

)3∇*<sup>p</sup>* (9)

*<sup>P</sup>* − *ρF*)**g** (11)

*ρPω* × (*ω* × **rP**) (12)

**<sup>P</sup>** (14)

**<sup>P</sup>**) (13)

<sup>2</sup> *dτ* (10)

(8)

*CD* <sup>=</sup> <sup>24</sup>

*Re<sup>i</sup>* <sup>=</sup> *<sup>d</sup><sup>i</sup>*

**pressure** <sup>=</sup> <sup>−</sup><sup>1</sup>

)2√*πρFμ<sup>F</sup>*

**Buoyancy** <sup>=</sup> <sup>3</sup>

**Fi**

coefficient of the particle.

particle velocities.

and is defined as

is given by

smaller than or similar to the fluid density.

**Fi**

here **g** expresses the gravity vector.

is the sum ob Coriolis and centripetal forces.

where, *ω* denotes rotation speed of the reference.

*mi P d***u<sup>i</sup> P** *dt* <sup>=</sup> <sup>1</sup> 8 *π*(*d<sup>i</sup>*

continuous phase influences the disperse phase.

**Frotation** <sup>=</sup> <sup>−</sup><sup>1</sup>

**Basset** <sup>=</sup> <sup>3</sup>

to the weight of the displaced fluid and is given by

2 (*di*

**Fi**

3 (*di*

small. The momentum conservation equation for the particle is shown as

Fig. 12. Three-dimensional nasal cavity model (case-1)


Table 1. Detail patient's case data

can neglect both the latent heat transfer and the radiative heat transfer. The convective heat transfer *QC* is given by

$$Q^i\_\mathbb{C} = \pi d^i \lambda\_F \text{Nu}(T\_F - T^i\_P) \tag{15}$$

where *T<sup>i</sup> <sup>P</sup>* is the temperature of the *i*-th particle, and Nu is the Nusselt number given by

$$Nu = 2 + 0.6 Re^{0.5} (\mu\_F \frac{\mathbf{C}\_p}{\lambda\_F})^{\frac{1}{5}} \tag{16}$$

The convective heat transfer *QC* is calculated and used in the source term of the thermal energy transport equation, Eq.4.

#### **5. Numerical simulations of intranasal transport phenomena**

#### **5.1 Case data and three-dimensional model**

Three-dimensional models of nasal cavity were constructed by means of five actual patient's CT data. These cases suffer several grades of deviated nasal septum symptom, hypertrophic rhinitis, chronic sinusitis and inferior nasal concha swelling. The case data are summarized in Table 1. Since head CT data involves not only nasal cavity but also skin, fat, bone and purulent matter, it is required to eliminate them. This study used three-dimensional medical image processing software "Mimics (Materialize Inc.)" and the marching cubes algorithm to determine the surfaces of nasal wall. The moving-average algorithm is also applied to smoothing the surface. Figure 12 shows the three-dimensional nasal model of case-I constructed in this study. The anterior and the posterior of the figure show anterior nares and upper pharynx, respectively. The three-dimension model has curved nasal septum which corresponds to the curvature as shown in Fig.12. The cross-sectional areas of nasal cavities

Transport Phenomena 13

Numerical Simulation for Intranasal Transport Phenomena 549

Fig. 15. Dependency between pressure drop and number of grid nodes in the CFD analysis

This study used CFD analysis code "CFX-11 (ANSYS, Inc.)" to analyze the intranasal transport phenomena of the nebulizer treatment. This code is able to solve partial differential equations for mass, momentum and energy transportations with appropriate boundary conditions. Since the Reynolds number of intranasal flow was several hundred, normal laminar flow analysis was used in this study. The Lagrange approach is adopted to solve transportation of the droplets of medicinal mist. Then the model employs assumptions that a droplet is sensitive with fluid flow, accelerating by velocity difference between fluid flow and the droplet. Air inlet velocity is set at 0.5m/s. This velocity condition corresponds to the condition of actual nebulizer treatment. As shown in Fig.16, four air inlet angles are considered in this study; 90, 60, 30 and 25 degree. In the nebulizer treatment, medicinal mist is atomized and transported by dry air. This study considers dry air (298K, 0.1MPa) as a fluid medium. No-slip condition, where fluid velocity on the wall surface is set at zero, is adopted to the nasal wall. Droplet diameter of medicinal mist is uniformly 50 micro meters. This value also corresponds

Figures 17-19 shows stream lines in case-I left nasal when inlet velocity is 0.5m/s and nebulizer angles are 90, 60 and 30 degree. As we can see from the stream lines, changing nebulizer angle 60 to 30 degree changes the main flow of nebulizer from the superior nasal concha (upper-area of nasal cavity) to the inferior nasal concha (lower-area). In the all

Fig. 16. Schematic figure of inlet nebulizer angles; 90, 60 and 30 degree

**5.2 Computational fluid dynamics model**

to the condition of actual nebulizer treatment.

**5.3 Results and discussion**

for case-I

Fig. 13. Summaries of cross-sectional areas of nasal cavities for all cases

Fig. 14. Schematic drawing of CFD analysis grid (unstructured grid)

which significantly affect intranasal flow are summarized in Fig.13. The lumen and the wall of three-dimensional models are meshed using commercial meshing software ICEM CFD (ANSYS, Inc.) and the delaunay triangulation algorithm as shown in Fig.14.

Figure 15 indicates a example of the influence of the mesh quality on the accuracy of CFD analysis. In this figure, pressure drop between inlet and outlet and the number of grid nodes are used as the factors of the accuracy and the mesh quality, respectively. The pressure drop significantly depends on the number of grid nodes; the pressure drop shows occasional ups and downs in the case of small number of grid nodes, on the contrary it becomes stable in the case of large number of grid nodes.

12 Will-be-set-by-IN-TECH

Fig. 13. Summaries of cross-sectional areas of nasal cavities for all cases

Fig. 14. Schematic drawing of CFD analysis grid (unstructured grid)

case of large number of grid nodes.

(ANSYS, Inc.) and the delaunay triangulation algorithm as shown in Fig.14.

which significantly affect intranasal flow are summarized in Fig.13. The lumen and the wall of three-dimensional models are meshed using commercial meshing software ICEM CFD

Figure 15 indicates a example of the influence of the mesh quality on the accuracy of CFD analysis. In this figure, pressure drop between inlet and outlet and the number of grid nodes are used as the factors of the accuracy and the mesh quality, respectively. The pressure drop significantly depends on the number of grid nodes; the pressure drop shows occasional ups and downs in the case of small number of grid nodes, on the contrary it becomes stable in the

Fig. 16. Schematic figure of inlet nebulizer angles; 90, 60 and 30 degree

#### **5.2 Computational fluid dynamics model**

This study used CFD analysis code "CFX-11 (ANSYS, Inc.)" to analyze the intranasal transport phenomena of the nebulizer treatment. This code is able to solve partial differential equations for mass, momentum and energy transportations with appropriate boundary conditions. Since the Reynolds number of intranasal flow was several hundred, normal laminar flow analysis was used in this study. The Lagrange approach is adopted to solve transportation of the droplets of medicinal mist. Then the model employs assumptions that a droplet is sensitive with fluid flow, accelerating by velocity difference between fluid flow and the droplet. Air inlet velocity is set at 0.5m/s. This velocity condition corresponds to the condition of actual nebulizer treatment. As shown in Fig.16, four air inlet angles are considered in this study; 90, 60, 30 and 25 degree. In the nebulizer treatment, medicinal mist is atomized and transported by dry air. This study considers dry air (298K, 0.1MPa) as a fluid medium. No-slip condition, where fluid velocity on the wall surface is set at zero, is adopted to the nasal wall. Droplet diameter of medicinal mist is uniformly 50 micro meters. This value also corresponds to the condition of actual nebulizer treatment.

#### **5.3 Results and discussion**

Figures 17-19 shows stream lines in case-I left nasal when inlet velocity is 0.5m/s and nebulizer angles are 90, 60 and 30 degree. As we can see from the stream lines, changing nebulizer angle 60 to 30 degree changes the main flow of nebulizer from the superior nasal concha (upper-area of nasal cavity) to the inferior nasal concha (lower-area). In the all

Transport Phenomena 15

Numerical Simulation for Intranasal Transport Phenomena 551

Fig. 20. Trajectories of nebulizer droplets of the case-III; inflow velocity of 0.5m/s and inlet

Fig. 21. Trajectories of nebulizer droplets of the case-III; inflow velocity of 0.5m/s and inlet

Fig. 22. Trajectories of nebulizer droplets of the case-III; inflow velocity of 0.5m/s and inlet

of nebulizer angle, we can reduce the amount of medicinal droplet's deposition on nares and

This study focuses and uses three parameters to organize the results of each case; pressure drop of nebulizer flow, maximum velocity of intranasal flow and deposition ratio which expresses how many droplets impinge and deposit on nasal wall. Figures 23-25 are correlations among them for all patient's cases. The correlations of each case show similar

middle nasal concha, consequently transport wide area of nasal cavity.

angle of 90 degree)

angle of 60 degree)

angle of 30 degree)

Fig. 17. Results of intranasal flow of the case-I; velocity magnitude distribution in each cross-section and streamlines (inflow velocity of 0.5m/s and inlet angle of 90 degree)

Fig. 18. Results of intranasal flow of the case-I; velocity magnitude distribution in each cross-section and streamlines (inflow velocity of 0.5m/s and inlet angle of 60 degree)

Fig. 19. Results of intranasal flow of the case-I; velocity magnitude distribution in each cross-section and streamlines (inflow velocity of 0.5m/s and inlet angle of 30 degree)

nebulizer angle conditions flow velocity increases near the nares. Especially, the increasing flow velocity is most obvious in the inflammation regions. Circulation flow is found in the lower inlet angle conditions. This circulation flow will be able to extend residence times of airflow and medicinal droplet in nasal cavity, consequently enhance the effect of the nebulizer treatment.

Figure 20-22 shows trajectories of medicinal droplets is the case-III when inlet velocity is 0.5m/s and inlet angles are 90, 60 and 30 degree. Almost all droplets impinge and deposit near the nares and middle nasal concha where the cross-section areas are narrow caused by deviated nasal septum and inflammation. From these results, the nebulizer angle is able to control the impingement and deposition characterisitics of medicinal droplets. By the control 14 Will-be-set-by-IN-TECH

Fig. 17. Results of intranasal flow of the case-I; velocity magnitude distribution in each cross-section and streamlines (inflow velocity of 0.5m/s and inlet angle of 90 degree)

Fig. 18. Results of intranasal flow of the case-I; velocity magnitude distribution in each cross-section and streamlines (inflow velocity of 0.5m/s and inlet angle of 60 degree)

Fig. 19. Results of intranasal flow of the case-I; velocity magnitude distribution in each cross-section and streamlines (inflow velocity of 0.5m/s and inlet angle of 30 degree)

treatment.

nebulizer angle conditions flow velocity increases near the nares. Especially, the increasing flow velocity is most obvious in the inflammation regions. Circulation flow is found in the lower inlet angle conditions. This circulation flow will be able to extend residence times of airflow and medicinal droplet in nasal cavity, consequently enhance the effect of the nebulizer

Figure 20-22 shows trajectories of medicinal droplets is the case-III when inlet velocity is 0.5m/s and inlet angles are 90, 60 and 30 degree. Almost all droplets impinge and deposit near the nares and middle nasal concha where the cross-section areas are narrow caused by deviated nasal septum and inflammation. From these results, the nebulizer angle is able to control the impingement and deposition characterisitics of medicinal droplets. By the control

Fig. 20. Trajectories of nebulizer droplets of the case-III; inflow velocity of 0.5m/s and inlet angle of 90 degree)

Fig. 21. Trajectories of nebulizer droplets of the case-III; inflow velocity of 0.5m/s and inlet angle of 60 degree)

Fig. 22. Trajectories of nebulizer droplets of the case-III; inflow velocity of 0.5m/s and inlet angle of 30 degree)

of nebulizer angle, we can reduce the amount of medicinal droplet's deposition on nares and middle nasal concha, consequently transport wide area of nasal cavity.

This study focuses and uses three parameters to organize the results of each case; pressure drop of nebulizer flow, maximum velocity of intranasal flow and deposition ratio which expresses how many droplets impinge and deposit on nasal wall. Figures 23-25 are correlations among them for all patient's cases. The correlations of each case show similar

Transport Phenomena 17

Numerical Simulation for Intranasal Transport Phenomena 553

tendencies even though there are the individual differences such as the shape of nasal cavity and the state of nasal disease. Therefore these parameters are useful to organize the results mentioned above. The maximum velocity has direct proportional relationship with pressure drop. Meanwhile pressure drop and maximum velocity have strong correlation with deposition ratio. These results lead us to conclude that the nebulizer conditions have

This chapter presents not only the procedure of CFD analysis for intranasal transport phenomena but also medicinal droplets transport characteristic in patient's nasal cavity using actual CT/MRI data. Furthermore, the result of correlation analysis for the nebulizer treatment condition is also presented here. As for the results of CFD analysis, nebulizer angle, inlet velocity and size of droplet significantly affect on the intranasal transport phenomena. The condition of low nebulizer angle achieves that the medicinal droplets are transported all over the nasal cavity, that is this condition will increases therapeutic response to the nasal diseases. The correlation among velocity of intranasal flow, pressure drop and deposition ratio, which are main parameters for intranasal transportation of medicinal droplets in the nebulizer treatment, indicates that the basis of the intranasal transportation show similar characteristics in all patient's cases. These results mean that there is an optimum condition of nebulizer treatment for all patients who suffer nasal diseases. Further research of the intranasal transport phenomena for many patient's cases is required to survey the best

Since the development of the medical image processing technique has been in progress, the CFD analysis for human body, such as blood flow from the heart, respiration flow in the pharynx and the lung, will be widely spread and achieve tangible results in clinical practices

Akiyama, M. (2002). *Analysis of Two-Phase Flow Dynamics: Multi-Physics Flow Analysis*, Korona

Boyle, R. & Thomas, R. (1988). *Computer Vision: A First Course*, Blackwell Scientific

Haruna, S. (2003). *An endoscope surgical operation of chronic sinus infection*, Hokendojinsha Inc.,

Lindemann, J., Keck, T., Wiesmiller, K., Sander, B., Brambs, H.-J., Rettinger, G. & Pless, D.

Lindemann, J., Keck, T., Wiesmiller, K., Sander, B., Brambs, H.-J., Rettinger, G. & Pless,

Lorensen, W. E. & Cline, H. E. (1987). Marching cubes:, *Computer Graphics: A high Resolution*

(2004). A numerical simulation of intranasal air temperature during inspiration, *The*

D. (2005). Numerical simulation of intranasal airflow after radical sinus surgery,

same-level influence on the intranasal transport phenomena for all cases.

This work was supported by the grant from JSPS KAKENHI 19760114.

*American Journal of Otolaryngology* Vol.26(No.1): 175–180.

*3D Surface Construction Algorithm* Vol.21(No.4): 163–169.

**6. Conclusion**

treatment condition.

within the next several years.

Inc., Tokyo.

Tokyo.

Publications, New York.

*Laryngoscope* Vol.114(No.1): 1037–1041.

**7. Acknowledgement**

**8. References**

Fig. 23. Correlation analysis between maximum velocity and pressure drop

Fig. 24. Correlation analysis between deposition ratio and maximum velocity

Fig. 25. Correlation analysis between deposition ratio and pressure drop

tendencies even though there are the individual differences such as the shape of nasal cavity and the state of nasal disease. Therefore these parameters are useful to organize the results mentioned above. The maximum velocity has direct proportional relationship with pressure drop. Meanwhile pressure drop and maximum velocity have strong correlation with deposition ratio. These results lead us to conclude that the nebulizer conditions have same-level influence on the intranasal transport phenomena for all cases.

#### **6. Conclusion**

16 Will-be-set-by-IN-TECH

Fig. 23. Correlation analysis between maximum velocity and pressure drop

Fig. 24. Correlation analysis between deposition ratio and maximum velocity

Fig. 25. Correlation analysis between deposition ratio and pressure drop

This chapter presents not only the procedure of CFD analysis for intranasal transport phenomena but also medicinal droplets transport characteristic in patient's nasal cavity using actual CT/MRI data. Furthermore, the result of correlation analysis for the nebulizer treatment condition is also presented here. As for the results of CFD analysis, nebulizer angle, inlet velocity and size of droplet significantly affect on the intranasal transport phenomena. The condition of low nebulizer angle achieves that the medicinal droplets are transported all over the nasal cavity, that is this condition will increases therapeutic response to the nasal diseases. The correlation among velocity of intranasal flow, pressure drop and deposition ratio, which are main parameters for intranasal transportation of medicinal droplets in the nebulizer treatment, indicates that the basis of the intranasal transportation show similar characteristics in all patient's cases. These results mean that there is an optimum condition of nebulizer treatment for all patients who suffer nasal diseases. Further research of the intranasal transport phenomena for many patient's cases is required to survey the best treatment condition.

Since the development of the medical image processing technique has been in progress, the CFD analysis for human body, such as blood flow from the heart, respiration flow in the pharynx and the lung, will be widely spread and achieve tangible results in clinical practices within the next several years.

#### **7. Acknowledgement**

This work was supported by the grant from JSPS KAKENHI 19760114.

#### **8. References**


**Part 5** 

**Additional Important Themes** 

Monya, M., Yamamoto, T., Nakata, S., Nakashima, T., Yamamoto, T. & Suzuki, T. (2009). Cfd analysis for intranasal mass transportation in deviated septum case (japanese), *JSME Journal Series B* Vol.75(No.751): 175–180.

Nakahashi, K. (1995). *Grid Generation and Computer Graphics*, Univ. of Tokyo, Tokyo.


**Part 5** 

**Additional Important Themes** 

18 Will-be-set-by-IN-TECH

554 Fluid Dynamics, Computational Modeling and Applications

Monya, M., Yamamoto, T., Nakata, S., Nakashima, T., Yamamoto, T. & Suzuki, T. (2009). Cfd

Savitzky, A. & Golay, M. J. E. (1964). Smoothing and differentiation of data by simplified least

Sun, C. (2006). Moving average algrithms for diamond, hexagon, and general polygonal shaped window operation, *Pattern Recognition Letters* Vol.27(No.6): 1627–1639.

Weinhold, I. & Mlynski, G. (2004). Numerical simulation of air flow in the human nose, *Eur*

Yamamoto, T., Nakata, S., Nakashima, T. & Yamamoto, T. (2009). Computational fluid

dynamics simulation for intranasal flow (japanese), *Foresight of Otorhinolaryngology*

Nakahashi, K. (1995). *Grid Generation and Computer Graphics*, Univ. of Tokyo, Tokyo.

squares procedures, *Anal. Chem.* Vol.36(No.8): 1627–1639.

*Journal Series B* Vol.75(No.751): 175–180.

Tamura, H. (2002). *Computer Image Processing*, Ohm Inc., Tokyo.

*Arch Otorhinolaryngol* Vol.261(No.1): 452–455.

Vol.52(No.1): 24–29.

analysis for intranasal mass transportation in deviated septum case (japanese), *JSME*

**25** 

*Italy* 

**Fluid-Dynamic Characterization and** 

**of the Hydraulic Separator Multidune** 

Floriana La Marca, Monica Moroni and Antonio Cenedese

*DICEA - Sapienza University of Rome, Rome* 

**Efficiency Analysis in Plastic Separation** 

Recovery of useable plastics from post-consumer and manufacturing waste remains a major recycling challenge. The global consumption of plastics was reported to be 230 million tonnes in 2005 (PlasticsEurope 2007a) of which 47.5 million tonnes were produced in Europe (25 European Union countries + Norway and Switzerland). Of the European production, only 22 million tonnes were reported as having been collected. Of this collected waste, 4 million tonnes were recycled as a manufacturing feedstock (18%) and 6.4 million tonnes went into energy recovery (29%), with the balance (11.6 million tonnes) probably being

The recycling of plastics is a process essential to reduce the efflux of materials to landfills and to decrease the production of raw materials. In recent years awareness of the importance of environmental protection has led to the development of different techniques for plastic recycling. One issue related to the recycling of this material is the presence in the market of many types of plastics (polymers with additives), often with similar

The separator "Multidune" is a hydraulic separator by density. Its name derives from the characteristic undulate profile of the channel where separation occurs. The channel is constructed from a sequence of closed parallel cylindrical tubes welded together in plane which are then sliced down the lateral mid-plane and the lower complex is laterally shifted relative to the upper complex. The Multidune allows solid particle separation according to

Previous investigations (De Sena et al., 2008; Moroni et al. 2008) suggested the flow within the Multidune is organized into three main patterns. Principally, a longitudinal transport flow takes place, where the velocity is high. A particle belonging to this region can move from one camera to another. The second region is the lower recirculation zone with high values of the vorticity field. Particles belonging to this region undergo the vertical impulse of the fluid. The thrust is proportional to the vertical velocity component and, in conjunction with gravity and buoyancy, determines the destiny of a particle. If the thrust is larger than the net weight of the particle, an interaction with the principal transport flow occurs and, consequently, the particle will move to the following chamber. The third region is the upper recirculation zone whose dimensions are smaller than the other recirculation zone. If a

characteristics that make them difficult to differentiate in the recovery phase.

their specific weight and the velocity field establishing within the apparatus.

**1. Introduction** 

disposed in landfills (PlasticsEurope 2007b).
