Analyzing Transient Turbulence in a Stenosed Carotid Artery by Proper Orthogonal Decomposition

Abstract

High-resolution three-dimensional simulations (involving 100 million degrees of freedom) were employed to study transient turbulent flow in a carotid arterial bifurcation with a stenosed internal carotid artery (ICA). The geometrical model was reconstructed from MRI images, and in vivo velocity measurements were incorporated in the simulations to provide inlet and outlet boundary conditions. Due to the high degree of the ICA occlusion and the variable flow rate, a transitional and intermittent flow between laminar and turbulent states was established. Time- and space-window proper orthogonal decomposition (POD) was applied to quantify the different flow regimes in the occluded artery. A simplified version of the POD analysis that utilizes 2D slices only—more appropriate in the clinical setting—was also investigated.

Appendix

Solver

For solution of Navier–Stokes equations we use the parallel numerical solver NEKTAR which is based on the high-order spectral element/hp method.17 The computational domain consists of structured or unstructured grids or a combination of both. In each element the solution is approximated in terms of hierarchical mixed-order Jacobi polynomials. NEKTAR employs up to a third-order splitting scheme,16 that decouples the velocity and pressure fields. The solution of three Helmholtz equations for the velocity components and Poisson equation for the pressure are obtained iteratively using a Preconditioned Conjugate Gradient solver; in this study a new parallel low energy basis preconditioner11,28 was implemented which has the parallel efficiency of a diagonal preconditioner but it is almost a ten-fold faster.

Geometry Reconstruction

A geometric model of the carotid artery was obtained from in vivo MRI images shown in Fig. 1a. In particular, MRI images of cross-sectional slices of the left common carotid (CCA), internal carotid (ICA) and external carotid (ECA) arteries were acquired by a Genesis Signa scanner using a 3D time-of-flight (TOF) sequence. The distance between the slices was 1.2 mm and a total of 145 images was acquired.

In order to obtain the required contours, the cross-sectional views of MRI images were processed by an in-house software package, a Matlab-based GUI which allows accurate detection of the arterial wall. The end product is a geometric model of the arterial wall saved in PLOT3D or STL format and is compatible with different mesh generation softwares. We use Gridgen—a commercial mesh generator developed by Pointwise,35 which also allows editing and construction of geometric databases. In Fig. 1a the medical image of stenosed carotid artery is presented while in Fig. 1b we show the geometric model reconstructed from MRI images.

In the current study, a mesh with 22,441 tetrahedral elements was created; in the reference element a Cartesian grid based on Gauss quadrature points is generated. We employ Gauss–Lobatto quadrature points in one direction and Gauss–Radau points in the other two directions and compute the number of degrees of freedom (DOF) in each element as a number of quadrature points required for exact integration of quadratic terms. For eighth-order polynomial expansion we have DOF = (P + 3) * (P + 2)² = 1100 per element, that is a total number of DOF per one variable in our simulation is 24,685,100. In addition, we employed the 3/2 over-integration rule for consistent integration and more accurate differentiation of the nonlinear term,18 thus the number of quadrature points in the domain was above 37M. Local mesh refinement was applied downstream of the ICA narrowing. The flat faces of the surface elements were projected on smooth curved boundaries. The information on the wall curvature was provided by patches of parametric surfaces and also by a surface smoothing technique described in Volino and Magnenat-Thalmann.33 In Fig. 1c we plot the parametric surface mesh (2D structured grid) at the carotid bifurcation and in Fig. 1d we provide an illustration of the computational grid, where the number of quadrature points inside each element is shown for a third-order polynomial approximation.

Resolution Study

Systematic resolution studies were performed. We used the h- and the p-type of refinement to minimize the error due to the spatial-discretization. The computational grid was refined in the regions with high velocity gradients; specifically, about 25% of the elements are located in a small fraction of ICA downstream the stenosis. The accuracy was verified by p-refinement. First, we performed unsteady 3D flow simulation in a pipe with Re = 350 and Ws = 4.375; a Womersley velocity profile was prescribed at the inlet, and zero pressure boundary condition at the outlet. Spectral convergence in Wall Shear Stress (WSS), presented in Fig. 15 (left), was obtained. Second, we performed steady flow simulation in the domain of carotid artery with parabolic flow profile at the inlet and zero pressure boundary conditions at both outlets. The L ₂-error in velocity field is presented in Table 1. The L ₂-error in the WSS, depicted in Fig. 15 (right), shows spectral convergence. In both cases the WSS computed with the highest resolution was treated as a reference solution. We note that the discretization error in WSS is at least one order higher than the error in the velocity field. In general, it is a very difficult task to obtain spectral convergence in the WSS in a domain with complex geometry and curved boundaries, e.g., the domain of carotid artery.

Figure 15

Verification of accuracy: convergence in the WSS. The magnitude of the (non-dimensional) WSS is of order O(1). Simulations performed with P = 4, 6, 8, 10, Δt = 0.0002. Solution with P = 10 is considered as a reference solution. Left: unsteady flow in pipe, Re = 350, Ws = 4.375. Right: steady flow in carotid artery, Re = 350 and Re = 700

Table 1 Verification of accuracy: steady flow simulations with Re = 350 and Re = 700

Boundary Conditions

The following boundary conditions were applied: At the inlet the Womersley velocity profile was imposed, consistent with the wave form of the flowrate range reported in Holdsworth.14 The main characteristic parameters are: mean Reynolds number at inlet Re _m = 350, the minimum and maximum of the Reynolds number at inlet are Re _min = 46 and Re _max = 1345; and the Womersley number Ws = 4.375.

At the outlets we used time-dependent RC-type boundary condition for the pressure supplemented with the zero Neumann boundary condition for velocity.10 This type of boundary condition provides an accurate and numerically efficient method to specify the flow rate ratio in numerical simulations of unsteady flow in a domain with multiple outlets. The pressure at outlet j is computed from

$$P_j(t) = -R_j(t)C_j\frac{dP_j}{dt} + R_j(t)Q_j(t), $$

(3)

where R _j (t), C _j and Q _j (t) are the specified resistance, capacitance and computed flow rate at the outlet j at time t. In a domain with two outlets, the required ratio Q ₁(t)/Q ₂(t) is achieved by fixing one of the resistance values, for example R ₁ = const, and computing the resistance at the second outlet from R ₂(t) = R ₁ f(t) where f(t) = Q ₁(t)/Q ₂(t). In this study we used R ₁ = R _ICA = 100, C ₁ = 0.002, C ₂ = 0.005; the capacitance C _j  is selected to effectively filter the high frequency flow rate oscillations at outlets which are due to numerical discretization and turbulence. The function f(t) was computed as a scaled ratio of velocity, measured with Doppler Ultrasound in ECA and ICA. The clinical data is presented in Fig. 16. The measurements of velocity in ICA and ECA were not synchronized, which results in the obvious difficulty on how to impose the correct phase shift between the two wave forms. In order to compensate for the missing parameters two simulations were performed. In the first simulation the systolic peak in ICA precedes the one in ECA whereas in the second simulation it lags. The waveforms of inlet and outlet flow rates for the two cases are presented in Fig. 17. The major difference between the two sets of waveforms is the phase shift between the flow rates at the two outlets; consequently the onset of turbulence and re-laminarization occurred at slightly different times with respect to the beginning of the cardiac cycle. In both simulations we observed transitional flow; window-POD analysis applied to two data sets produced similar results, hence in this paper we present results of the first simulation only.

Figure 16

(in color) Doppler ultrasound measurements: Velocity wave in internal (left) and external (right) carotid arteries. Left side of plots (a) and (c) shows the monitored velocity; the wide white line highlights its maximum. In the middle of each frame the 2D image of instantaneous flow in arterial segment is captured; it depicts the relative (to the flow) direction of the ultrasound signal, the two short lines positioned inside the vessel point to the spatial location where the velocity was monitored. (b) Geometry of the carotid artery; colors represent pressure

Figure 17

Simulation of intermittent laminar-turbulent-laminar flow in stenosed carotid artery with R(t) C-type boundary condition: flow rates in CCA, ICA and ECA. The upper plot shows the case where the peak velocity at ICA precedes the one in ECA; the lower plot shows the alternative scenario. Prescribed (reference) values of the flow rates are marked with lines; computed values of the flow rates are marked with circles. R _ICA = 100, R _ECA(t) = R _ICA(Q _ICA(t)/Q _ECA (t)), C _ICA = 0.002, C _ECA = 0.005. mean(Re _inlet) = 350, min(Re _inlet) = 46, max(Re _inlet) = 1345, Ws = 4.375