- Email: [email protected]
Predicting flow in aortic dissection: Comparison of computational model with PC-MRI velocity measurements
Medical Engineering & Physics 36 (2014) 1176–1184
Contents lists available at ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Predicting flow in aortic dissection: Comparison of computational model with PC-MRI velocity measurements Z. Cheng a , C. Juli c , N.B. Wood a , R.G.J. Gibbs b , X.Y. Xu a,∗ a b c
Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, UK Imperial Vascular Unit, St Mary’s Hospital, Imperial College Healthcare NHS Trust, London W2 1NY, UK Department of Radiology, St Mary’s Hospital, Imperial College Healthcare NHS Trust, London W2 1NY, UK
a r t i c l e
i n f o
Article history: Received 6 November 2013 Received in revised form 9 June 2014 Accepted 2 July 2014 Keywords: Aortic dissection Computational fluid dynamics Magnetic resonance velocity imaging
a b s t r a c t Aortic dissection is a life-threatening process in which the weakened wall develops a tear, causing separation of wall layers. The dissected layers separate the original true aortic lumen and a newly created false lumen. If untreated, the condition can be fatal. Flow rate in the false lumen is a key feature for false lumen patency, which has been regarded as one of the most important predictors of adverse early and later outcomes. Detailed flow analysis in the dissected aorta may assist vascular surgeons in making treatment decisions, but computational models to simulate flow in aortic dissections often involve several assumptions. The purpose of this study is to assess the computational models adopted in previous studies by comparison with in vivo velocity data obtained by means of phase-contrast magnetic resonance imaging (PC-MRI). Aortic dissection geometry was reconstructed from computed tomography (CT) images, while PC-MRI velocity data were used to define inflow conditions and to provide distal velocity components for comparison with the simulation results. The computational fluid dynamics (CFD) simulation incorporated a laminar–turbulent transition model, which is necessary for adequate flow simulation in aortic conditions. Velocity contours from PC-MRI and CFD in the two lumens at the distal plane were compared at four representative time points in the pulse cycle. The computational model successfully captured the complex regions of flow reversal and recirculation qualitatively, although quantitative differences exist. With a rigid wall assumption and exclusion of arch branches, the CFD model over-predicted the false lumen flow rate by 25% at peak systole. Nevertheless, an overall good agreement was achieved, confirming the physiological relevance and validity of the computational model for type B aortic dissection with a relatively stiff dissection flap. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction Aortic dissection is a life-threatening process in which a tear in the aortic wall allows blood to track within the wall, separating its component layers. The dissected layers (dissection lamellae or flap) separate the original true aortic lumen and a newly created false lumen. The anatomical position of the entry tear is used to define Type A dissection (starting in the ascending aorta) and Type B dissection (starting in the distal aortic arch). Type A dissection normally requires immediate surgery, whereas for Type B the decision for medical treatment or surgical intervention depends on the clinical presentation. In these cases, detailed analysis of flow patterns in the dissection aorta can provide useful information for clinical decision-making.
∗ Corresponding author. E-mail address: [email protected] (X.Y. Xu). http://dx.doi.org/10.1016/j.medengphy.2014.07.006 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.
The development of high-resolution medical imaging and robust image processing techniques has promoted image-based computational fluid dynamics (CFD) modelling of in vivo haemodynamics [1]. A number of studies have been carried out to investigate the flow patterns and haemodynamic environment in dissected aorta. Karmonik et al. [2] examined changes in anatomy as well as pressure and wall shear stress in an aortic dissection case with aneurysmal dilatation. In a longitudinal study, Tse et al. [3] investigated variations of certain haemodynamic parameters in the development of a dissecting aneurismal aorta. In our previous studies, an advanced transitional turbulence model was employed to investigate the complex flow patterns and detailed haemodynamic environment in multiple cases of type B dissection [4,5]. Results were analysed in each case and compared between different groups. One of the most important findings was that flow rate into the false lumen was strongly correlated with the primary entry tear size. Since false lumen patency is a key predictor of adverse early and later outcomes [6], false lumen flow rate (or entry
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1177
Fig. 1. Reconstructed geometries of the aortic dissection used in validation. (a) The entire aorta including the ascending aorta lumen (AL), the aortic arch branches, the descending aorta lumen (DL), and the iliac arteries (IA), but excluding the branches off the descending thoracic and abdominal aorta; (b) the locations of three imaging planes where phase-contrast MR velocity mapping was acquired; and (c) the computational model used in validation, with arrows indicating the locations of inlet, outlet, tear, false, and true lumen.
tear size) may be an important parameter for risk stratification of dissection patients. However, computational models for the analysis of flow in aortic dissections (e.g. [2–5]) often involve several assumptions, including compensated exclusion of the aortic arch branches, flat velocity profile at the inlet, and rigid walls accompanied by fixed, constant distal outflow pressures. The purpose of this study is to assess the computational models adopted in our previous studies by comparison with in vivo velocity data obtained by means of phase-contrast magnetic resonance imaging (PC-MRI). Computed tomography (CT) is the standard imaging modality for diagnostic imaging of aortic dissection, but for flow simulations, 4D ciné PC-MRI would be desirable to provide patient-specific flow boundary conditions [7]. PC-MRI flow imaging is considered the only imaging technique that can provide measurement of blood velocities in vivo with good resolution both spatially and temporally. In this study, a typical type B aortic dissection case was analysed by using the same computational model design and numerical analysis methods as adopted in our previous studies [4,5]. This case was chosen for its comprehensive set of PC-MRI data which provided not only the necessary boundary conditions, but also additional data required for validation. 2. Methods 2.1. Medical image processing and geometry reconstruction Anatomical information was extracted from CT scan of a patient (male, age 56) diagnosed with a typical type B aortic dissection. Subsequently, PC-MRI data were acquired for further examination. The study complied with the Declaration of Helsinki and was approved by the local Research Ethics Committee. The patient gave written, informed consent. The model geometry was based on the initial CT images covering the entire aorta (approximately 300 slices, 2 mm thick, spacing 1 mm; pixel resolution 0.86 mm), as shown in Fig. 1a. The proximal entry tear can be seen clearly in the transparent image (Fig. 1b). Anatomical MR images had relatively low axial resolution and were only used for registration and localization of imaging planes
between the CT and PC-MRI data. The CT images were imported into an image-processing package Mimics (Materialise HQ, Louvain, Belgium) and the geometry of the dissected aorta was reconstructed from the contours of the segmented lumen area from the serial CT slices and was smoothed accordingly before meshing. PC-MRI was performed at three planes in the aortic dissection (Fig. 1b), providing sufficient data for inlet boundary condition and validation (TE 2.6 ms, TR 4.7 ms, flip angle 15◦ , pixel spacing 1.3281 mm, slice thickness 8 mm). At the imaging planes 1 and 3, encoding velocity (Venc ) was 200 mm/s for IS (inferior–superior) direction and 100 mm/s for AP (anterior–posterior) and RL (right–left) directions. At imaging plane 2, Venc was 200 mm/s for AP direction, and 100 mm/s for the other two directions. Plane 1 was located in the proximal ascending aorta in order to obtain velocity profiles distal to the aortic valve. Plane 2 was in the aortic arch between the root of the left subclavian artery and the site of proximal entry tear. Since the tear occurred very close to the left subclavian artery, there was limited space to fit plane 2 between the artery and the entry tear. Flow imaging at plane 2 was required to determine the amount of flow distal to the arch branches, and velocity profiles proximal to the tear. Plane 3 was in the mid-thoracic descending aorta distal to the proximal entry tear, at the level of pulmonary bifurcation, covering both the true and false lumen. The PC-MRI data from this plane were used to compare with those from the computational simulation for model validation. Images acquired with PC-MRI were processed by using an in-house programme. The region of interest was segmented from magnitude images and then mapped on to the corresponding phase images; the grayscale of the pixels inside the region of interest was converted into velocity information according to the applied encoding sequence of the scanner (Philips Medical Systems). For imaging planes 1 and 2, the flow rate through each plane was calculated based on the converted velocity and pixel size; this was utilized as boundary conditions for the flow model described in the next section. The PC-MRI data acquired at plane 3 were processed in the same way, and 2D velocity contours of the true and false lumen were plotted for the purpose of comparison with the simulation results.
1178
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
The computational model (Fig. 1c) shows the inlet at imaging plane 1. To keep the approach consistent with previous patientspecific studies, the model was simplified by excluding the three branches in the aortic arch and other small branches in the descending aorta, whilst matching the flow at planes 1 and 2 to provide the correct amount of flow through planes 2 and 3. The false lumen propagated down the lateral side of the thoracic descending aorta, before spiralling around the true lumen toward its anterior side at the mid-descending aorta. It continued to twist to the rightposterior side of the aorta, then to the lateral side in the distal segment of the abdominal aorta. As it extended to the left iliac artery below the abdominal bifurcation, two adjacent outlets (one for each lumen) were included in the model, both located in the abdominal aorta far distal to the imaging plane 3 (Fig. 1c). The length of the descending aorta was approximately 381 mm. The proximal entry tear was about 22 mm in width and 25 mm in height, and only 0.4 mm below the top of the aortic arch. The false lumen was prevented from propagating proximally by the left subclavian artery. No distal tears could be identified. 2.2. Flow simulation and boundary conditions Simulations were carried out with ANSYS CFX 12, and the – correlation-based laminar–turbulent transition model (‘SST Tran’) was adopted [8], as described previously [4]. Blood was assumed to be Newtonian with viscosity 4.0 mPa s and density 1060 kg/m3 . The aortic walls were assumed to be rigid with no-slip. Flow boundary conditions were extracted from the PC-MRI velocity images. The flow data at each imaging plane were recorded at 29 time points along a cardiac cycle of 0.98 s duration. The flow rate waveforms at imaging plane 1 (ascending aorta root) and imaging plane 2 (aortic arch after the branches) are presented in Fig. 2. The peak systolic flow rate was 446 cm3 /s in the ascending aorta, and 370 cm3 /s in the distal aortic arch. Hence, the flow exiting through the three main branches on the arch at peak systole was 76 cm3 /s, which was 17% of the peak flow rate from the heart. This ratio is close to the value (15%) assumed previously [4,5]. The flow waveform derived from imaging plane 2 was specified at the model inlet in the ascending aorta with a flat through-plane profile, as adopted in previous patient-specific studies. Based on the inlet area, the mean and maximum Reynolds numbers were 850 and 3512, respectively, and the Womersley parameter was 10.4. A low inlet turbulence level of 1.5% was chosen to represent initial disturbances in the flow and allow true transition to occur realistically. At the model outlets, zero relative static pressure was applied across both surfaces. A uniform time step of 0.001 s was used. All simulations were carried out for three cardiac cycles to achieve a
Fig. 2. The flow rate waveforms at imaging plane 1 and plane 2 as defined in Fig. 1b. The flow rate was calculated based on the converted velocity at each pixel and the pixel size of the PC-MRI data.
periodic solution, and the results of the last cycle were used for comparison. 2.3. Mesh sensitivity test The reconstructed model was imported into ANSYS ICEM CFX for mesh generation. Unstructured mesh consisting of 3D tetrahedral cells combined with prismatic cells near the wall was created. A very fine near wall resolution in the tear region was ensured by using a minimum of 10 layers of prismatic cells at the wall, height of wall cells in wall units (y+) < 2, together with local refinement. Mesh independence tests were carried out with an initial mesh consisting of 2,600,000 elements and finer meshes with increased number of prismatic wall elements. The results in terms of maximum timeaveraged wall shear stress (WSS) and turbulence kinetic energy were compared between different meshes. The differences in these parameters were less than 4% between the initial 2,600,000 elements mesh and the 3,200,000 elements mesh, and less than 3.5% between the 3,200,000 elements mesh and the 4,700,000 elements mesh. Hence, the mesh of 3,200,000 elements was deemed to be sufficient. 2.4. Analysis of simulation results The simulation results were analysed using ANSYS CFX-Post 12. For the purpose of comparison, velocities are displayed by plotting colour-coded contours at the cross-sectional plane corresponding to imaging plane 3 where PC-MRI data were acquired. Instantaneous streamlines are presented at different time points for better visualisation of time-varying flow patterns. In order to examine the turbulence level in the model, turbulence intensity (Tu) was evaluated, which is defined as:
Tu =
2 k/V 3
(1)
where k is the turbulence kinetic energy and V is the instantaneous local velocity. The iso-surface of Tu is presented to indicate the corresponding turbulence zone with specified value at the systolic mid-deceleration phase, when Tu is usually at its highest value. 3. Results 3.1. Flow distribution in the true and false lumen As reported by Cheng et al. [5], the flow division between the two lumens was strongly dependent on the size and location of the proximal entry tear. Fig. 3 shows the comparison of timevarying flow rate at the location of imaging plane 3 obtained from
Fig. 3. The flow rate waveforms of the true and false lumen from both CFD results and PC-MRI velocity data at imaging plane 3.
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1179
Fig. 4. Comparison of velocity profiles in FH, AP, and RL directions at mid-systolic acceleration phase of the cardiac cycle between PC-MRI velocity data (in m/s) and simulation results. In Figs. 4–7 the red colour represents positive velocities from superior to inferior, anterior to posterior (IS contours), right to left (AP and RL contours); the blue colour represents negative velocities. Note different colour scales for different directions and time points according to the range of velocity values. The secondary velocity vectors are also presented as project on the plane.
the CFD simulation and the corresponding PC-MRI data. Based on the simulation results, the peak flow rate in the false lumen was 293 cm3 /s, and 72 cm3 /s in the true lumen, whereas the cycleaveraged values were 101.6 and 27.1 cm3 /s, respectively. Hence, about 80% of the peak flow was diverted to the false lumen. As shown in Fig. 3, the measured and computed flow rate waveforms are similar in shape but the CFD model over-predicted flow in the false lumen and under-predicted flow in the true lumen. The differences between the computed and measured flow rates were 8.5% for the true lumen and 25% for the false lumen at peak flow rate (the cycle-averaged values were 35% for the true lumen and 28% for the false lumen). A number of factors could have contributed to the differences; delineation of the lumen contours on PC-MRI images; motion of the intra-arterial septum; limited spatial resolution, noise, and apparent artefacts associated with the PC-MRI velocity data. 3.2. Comparison of flow structure and velocity profiles Comparisons were made between the time-varying velocity contours from PC-MRI measurements and flow simulation results at plane 3 (defined in Fig. 1b) for all three velocity components. Figs. 4–7 present the velocity contours in the main inferior–superior (labelled FH, foot-head), secondary AP (anterior–posterior), and RL (right–left) directions extracted from
the PC-MRI data and simulations at four time points along the cardiac cycle. Both sets of data are displayed in the conventional MRI direction, looking towards the head. The corresponding magnitude images from the PC-MRI are included, indicating the shape and size of the false (FL) and true lumen (TL). As can be seen, the FL is approximately four times larger than the TL at this location. For better understanding of the complex flow patterns, secondary velocity vectors have been superimposed on the FH velocity contours. Complementary instantaneous streamlines for the four time points are shown in Fig. 8. At mid-systolic acceleration, the velocity contours (Fig. 4) show that the main direction of the flow in the TL was from superior to inferior (distal) with a maximum velocity of 0.45 m/s based on PCMRI. The high velocity flow was skewed towards the left-anterior side of the TL owing to the curvature of the arch; a similar pattern can be seen in the simulation result. The flow in the FL showed more spatial variation than in the TL. The simulation results captured the relatively high velocity at the left anterior side near the septum and in the middle of the FL. The corresponding secondary velocity vectors show that in the TL the direction of secondary flow was predominantly anterior–posterior, pointing to the right. Similar patterns can also be observed in the AP and RL contours. The proportion of the flow remaining in the TL was accelerated by the narrowing of the lumen, so stretching the helical pattern induced by the non-planar arch. In the FL the pattern was more complex,
1180
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
Fig. 5. Comparison of velocity profiles (in m/s) between PC-MRI velocity data and simulation results in FH, AP, and RL directions at peak systole.
with a weak vortex centred near the right (R) wall, but extending across most of the area. The vortex originated from the flow entering at the tear proximally to this plane, curving along the opposite wall, setting up a large proximal vortex (Fig. 8a), which expanded to fill the lumen by half-way along the descending thoracic aorta. At peak systole (Fig. 5), the PC-MRI flow patterns in the TL showed fairly uniform high velocities at this location, which were well captured by the CFD simulation. The corresponding streamlines and secondary velocity vectors again reflected the deviation of the TL axis from the FH direction. In the FL the pattern had changed in response to the higher velocities through the tear, shown clearly in Fig. 8b. The vortex previously centred near the right wall had moved anteriorly and the effect of the stronger jet deflected along the left wall could be seen (recall flow is shown from the inferior aspect, as with the MRI convention). The high velocities along the left wall resulted from the jet through the tear as it impacted the opposite wall and followed the curvature of the arch; at plane 3 (Fig. 8b) the jet had spread to fill the cross-section, with lower velocities towards the right-anterior wall resulting from the greater expansion of the jet distal to the separated, recirculating flow immediately below the tear. At mid-systolic deceleration (Fig. 6), flow in the TL was still dominated by the superior–inferior component, and the velocity distribution was similar to that at peak systole, with a good agreement between measured and simulated velocities. There was a subtle change in flow near the anterior wall as a result of the helical flow in the aortic arch induced by flow deceleration, as seen in
the instantaneous streamlines (Fig. 8c). The effect of the retardation was clearer in the FL. The flow distal to the weaker jet through the tear was more disorganised, as shown in Fig. 8c, whilst in the region of Plane 3 the somewhat disorganised flow became more evident. In Fig. 6, the simulation showed the high velocity zone having moved farther along the right wall than the measurements, while the secondary vectors revealed that the main vortex had moved towards the centre of the area and a weak secondary vortex had formed at the left-posterior wall. The simulation predicted significant turbulence in the FL distal to the tear, particularly in late systole (Fig. 9). Comparisons of velocities in early diastole, when the maximum retrograde flow occurred in the proximal ascending aorta, are given in Fig. 7. All velocity contours from MR and CFD showed excellent agreement in both lumens. The TL was dominated by a helical flow in the superior–inferior direction, despite the area mean flow in the ascending aorta being retrograde. The corresponding streamlines and secondary vectors revealed a complex flow structure in the arch (Fig. 8d). Flow in the FL at Plane 3 showed retrograde flow on the left-anterior area, occupying two-thirds of the crosssectional area, whilst it was antegrade in the remaining third, at the right-posterior aspect (Fig. 7). The streamlines reflected this feature, where it would appear from Fig. 8d that part of the flow changed direction approximately one diameter distal to Plane 3. Quantitative comparisons were made for the maximum velocities in the TL and FL between the PC-MRI data and CFD results. The maximum difference of instantaneous velocities in both TL and
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1181
Fig. 6. Comparison of velocity profiles (in m/s) between PC-MRI velocity data and simulation results in FH, AP, and RL directions at systolic deceleration phase.
FL was found at the systolic mid-deceleration phase, when the TL velocity was 29% higher with the PC-MRI measurement than with the CFD prediction, whilst the FL velocity was 30% higher with the CFD than with the PC-MRI. For the overall time-averaged comparison, the equivalent differences were 20 and 16%. Reasons for the over-prediction of FL velocity have been explained in the preceding section. 4. Discussion There have been a few limited examples of such studies in the literature; Tse et al. [3] examined a dissected aorta pre- and postaneurysm, whilst Karmonik et al. [2] reported on a chronic type III dissection both at initial examination and at 10-month followup. Neither of the studies incorporated a transitional or turbulence model in their CFD simulations. Based on the stability diagram adopted by Kousera et al. [9], the combination of peak Reynolds number and Womersley parameter considered here is likely in the region of disturbed flow, justifying the need for a transitional flow model. With regard to boundary conditions. Karmonik et al. [2] used patient-specific flow waveform measured with PC-MRI as the inlet boundary condition, while Tse et al. [3] applied a typical aortic flow waveform at the inlet and a corresponding pressure waveform at the exits in the iliac arteries. Although Karmonik et al. [2] showed comparison of PC-MRI and CFD predicted velocities at a transverse plane in the descending aorta, it was for the axial through-plane velocity magnitudes only. Detailed assessment of CFD predictions using MR velocity imaging in three orthogonal directions is lacking. False lumen flow rate is a key predictor of FL patency, being the most important predictor of adverse early and later outcomes
[6]. An interesting finding of Cheng et al. [5] was the correlation between the FL flow rate and the size and location of the primary entry tear, based on the computational results of four patient-specific models with the same assumptions and simplified boundary conditions as employed here. Comparison of the FL flow rate obtained from the CFD simulation and that from the PC-MRI data (Fig. 3) showed that the computational model gave comparable results to the MR measurements with similar flow waveforms for each lumen, and the percentage of flow in the FL (80% vs. 76%). The size and location of the primary entry tear and the percentage of FL flow rate extracted from the MR data were found to fit well to the regression lines presented by Cheng et al. [5]. Further careful comparison was made between the 4D PC-MRI velocity data and the simulation results at the proximal middescending aorta. The velocity data were compared in each of the three components. It was found that (Figs. 4–7) velocity contours derived from the CFD simulation were in good agreement with those from PC-MRI. No retrograde flow was found in the TL, as the dominant velocity component (in the IS or FH direction) was directed from head to feet throughout the cardiac cycle. The other two velocity components also showed consistent patterns in pointing towards the posterior and right directions, respectively, caused mainly by the inclination of the TL to the FH direction. The organized flow pattern in the TL, which was also demonstrated in previous studies [2–5], was related to its compressed size and relatively straight geometry. Flow patterns in the FL were much more complicated, seen both in the contours (Figs. 4–7) and related instantaneous streamlines (Fig. 8). The main flow was generally in the superior–inferior direction and the flow velocity was lower in the other two directions.
1182
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
Fig. 7. Comparison of velocity profiles (in m/s) between PC-MRI velocity data and simulation results in FH, AP, and RL directions at maximum retrograde flow in early diastole.
At mid-systolic acceleration, no clear velocity pattern can be observed in the FH direction but CFD results successfully captured the patterns in the AP and RL directions. By peak systole, the FH velocity was higher near the left wall where CFD and PC-MRI data showed the same pattern. However, the other two directions showed weaker agreement. One likely reason for this was the patchy nature of the PC-MRI data which was especially amplified when velocity was low. Another key reason might be the rigid wall assumption in the CFD simulations. The movement of the dissection flap is likely to be most significant at peak systole when the pressure difference between the TL and FL is at the highest [10]. Although the intimal flap was found to have been calcified and quite stiff in the particular case concerned here, small motions could still happen, hence affect the flow velocity near the flap especially in the AP and RL directions. During retardation (Fig. 6), the agreement between simulations and PC-MRI was good for all velocity components. The TL flow remained antegrade, whilst there was a substantial retrograde zone in the FL at plane 3. It could be discerned from the streamlines in Fig. 8c that flow was largely disturbed and recirculating. There was a good agreement between PC-MRI and simulations at the maximal retrograde flow in early diastole as seen in the velocity contours of Fig. 7. Considerably disorganised flow in the FL was also evident in Fig. 8d. In the region of the tear, there was major recirculation but, because the septum does not lie perpendicular to the oblique sagittal plane, it remains difficult to discern the boundary between FL and TL.
In summary, the computational model successfully captured the complex regions of flow reversal and recirculation qualitatively, although quantitative differences exist. The overall good degree of agreement confirmed the validity of the computational model utilized in our previous case studies, providing reliable data on flow patterns and haemodynamic parameters.
4.1. Limitations The CFD model presented here involved several assumptions. First, the aortic wall and dissection flap were assumed to be rigid. While this assumption was acceptable in the case discussed here owing to its stiff flap, this might not be true in cases when the flap is more compliant. The dissection flap is initially very mobile immediately after the dissection has occurred, but becomes more rigid over time. Ganten et al. [11] examined the intimal flap and wall motion in 32 type B aortic dissection patients by using dynamic, ECG-gated CT imaging technique. It was found that distensibility of the dissected descending aorta was significantly lower than that in healthy subjects of the same age, and the average true lumen narrowing by intimal flap motion was 4.4% (maximum 29%). However, in 14 of the 32 patients there was no detectable flap motion. Another study based on a smaller number of cases demonstrated the highly variable flap motion in type B dissections [10], with up to 30% change in true lumen area. Using 2D PC-MRI in a longitudinal
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1183
Fig. 8. Instantaneous streamlines of blood flow from left posterior aspect at (a) systolic mid-acceleration phase, (b) peak systole, (c) systolic mid-deceleration phase, and (d) maximum retrograde flow in diastole.
Fig. 9. Turbulence intensity isosurfaces at mid-systolic deceleration from left-posterior (left) and right-anterior (right) aspects.
1184
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
study, Karmonik et al. [12] found that flap motion in type B aortic dissections reduced over time. The exclusion of the arch branches is another assumption in this study. To compensate for the flow loss through the three branches, the inlet flow rate at the aortic root was scaled down to maintain the correct amount of flow into the descending aorta. This appeared to give satisfactory results, although the flow structure in the entry tear region would have been compromised. However, owing to the nature of the geometry of dissected aorta, the flow pattern is significantly changed after passing through the tear into the false lumen, as well as in the highly compressed true lumen. In conclusion, we have demonstrated good agreement in flow structure and velocity contours in the dissected aorta, proving the validity of the computational model and the simulation approach for type B aortic dissection with a relatively stiff dissection flap. Future studies are needed to examine the effect of flap motion on flow in the true and false lumen. Declarations Funding: Zhuo Cheng was supported by a PhD scholarship from the Institute of Biomedical Engineering and the National Heart and Lung Institute, Imperial College London. Competing Interests: There is no conflict of interest that could inappropriately influence this research work. Ethical Approval: Ethical approval was given by the local Research Ethics Committee at the St Mary’s Hospital. References [1] Steinman DA. Image-based computational fluid dynamics modeling in realistic arterial geometries. Ann Biomed Eng 2002;30(4):483–97.
[2] Karmonik C, Partovi S, Müller-Eschner M, Bismuth J, Davies MG, Shah DJ, Loebe M, Böckler D, Lumsden AB, von Tengg-Kobligk H. Longitudinal computational fluid dynamics study of aneurysmal dilatation in a chronic DeBakey type III aortic dissection. J Vasc Surg 2012;56:260–3. [3] Tse KM, Chiu P, Lee HP, Ho P. Investigation of hemodynamics in the development of dissecting aneurysm within patient-specific dissecting aneurismal aortas using computational fluid dynamics (CFD) simulations. J Biomech 2011;44:827–36. [4] Cheng Z, Tan FPP, Riga CV, Bicknell CD, Hamady MS, Gibbs RGJ, Wood NB, Xu XY. Analysis of flow patterns in a patient-specific aortic dissection model. J Biomech Eng 2010, 132051007-1-9. [5] Cheng Z, Riga CV, Chan J, Hamady MS, Wood NB, Cheshire NJW, Xu XY, Gibbs RGJ. Computational simulation of the aorta in acute type B dissection: initial findings and potential applicability. J Vasc Surg 2013;57: 35S–43S. [6] Eggebrecht H, Herold U, Kuhnt O, Schmermund A, Bartel T, Martini S, Lind A, Naber CK, Kienbaum P, Kuhl H, Peters J, Jakob H, Erbel R, Baumgart D. Endovascular stent-graft treatment of aortic dissection: determinants of postinterventional outcome. Eur Heart J 2005;26:489–97. [7] Clough RE, Waltham M, Giese D, Taylor PR, Schaeffter T. A new imaging method for assessment of aortic dissection using four-dimensional phase contrast magnetic resonance imaging. J Vasc Surg 2012;55:914–23. [8] Menter F, Langtry R, Volker S. Transition modelling for general purpose CFD codes. Flow Turbulence Combust 2006;77:277–303. [9] Kousera CA, Wood NB, Seed WA, Torii R, O’Regan D, Xu XY. A numerical study of aortic flow stability and comparison with in vivo flow measurements. J Biomech Eng 2013;135, 011003-1-9. [10] Murayyama T, Funabashi N, Uehara M, Takaoka H, Komuro I. New classification of aortic dissection during the cardiac cycle as pulsating type and static type evaluated by electrocardiogram-gated multislice CT. Int J Cardiol 2010;142:177–86. [11] Ganten MK, Weber TF, von Tengg-Kobligk H, Bockler D, Stiller W, Geisbusch P, Kauffmann GW, Delorme S, Bock M, Kauczor HU. Motion characterization of aortic wall and intimal flap by ECG-gated CT in patients with chronic Bdissection. Eur J Radiol 2009;72:146–53. [12] Karmonik C, Duran C, Shah DJ, Anaya-Ayala JE, Davies MG, Lumsden AB, Bismuth J. Preliminary findings in quantification of changes in septal motion during follow-up of type B aortic dissections. J Vasc Surg 2012;55: 1419–26.
Contents lists available at ScienceDirect
Medical Engineering & Physics journal homepage: www.elsevier.com/locate/medengphy
Predicting flow in aortic dissection: Comparison of computational model with PC-MRI velocity measurements Z. Cheng a , C. Juli c , N.B. Wood a , R.G.J. Gibbs b , X.Y. Xu a,∗ a b c
Department of Chemical Engineering, Imperial College London, South Kensington Campus, London SW7 2AZ, UK Imperial Vascular Unit, St Mary’s Hospital, Imperial College Healthcare NHS Trust, London W2 1NY, UK Department of Radiology, St Mary’s Hospital, Imperial College Healthcare NHS Trust, London W2 1NY, UK
a r t i c l e
i n f o
Article history: Received 6 November 2013 Received in revised form 9 June 2014 Accepted 2 July 2014 Keywords: Aortic dissection Computational fluid dynamics Magnetic resonance velocity imaging
a b s t r a c t Aortic dissection is a life-threatening process in which the weakened wall develops a tear, causing separation of wall layers. The dissected layers separate the original true aortic lumen and a newly created false lumen. If untreated, the condition can be fatal. Flow rate in the false lumen is a key feature for false lumen patency, which has been regarded as one of the most important predictors of adverse early and later outcomes. Detailed flow analysis in the dissected aorta may assist vascular surgeons in making treatment decisions, but computational models to simulate flow in aortic dissections often involve several assumptions. The purpose of this study is to assess the computational models adopted in previous studies by comparison with in vivo velocity data obtained by means of phase-contrast magnetic resonance imaging (PC-MRI). Aortic dissection geometry was reconstructed from computed tomography (CT) images, while PC-MRI velocity data were used to define inflow conditions and to provide distal velocity components for comparison with the simulation results. The computational fluid dynamics (CFD) simulation incorporated a laminar–turbulent transition model, which is necessary for adequate flow simulation in aortic conditions. Velocity contours from PC-MRI and CFD in the two lumens at the distal plane were compared at four representative time points in the pulse cycle. The computational model successfully captured the complex regions of flow reversal and recirculation qualitatively, although quantitative differences exist. With a rigid wall assumption and exclusion of arch branches, the CFD model over-predicted the false lumen flow rate by 25% at peak systole. Nevertheless, an overall good agreement was achieved, confirming the physiological relevance and validity of the computational model for type B aortic dissection with a relatively stiff dissection flap. © 2014 IPEM. Published by Elsevier Ltd. All rights reserved.
1. Introduction Aortic dissection is a life-threatening process in which a tear in the aortic wall allows blood to track within the wall, separating its component layers. The dissected layers (dissection lamellae or flap) separate the original true aortic lumen and a newly created false lumen. The anatomical position of the entry tear is used to define Type A dissection (starting in the ascending aorta) and Type B dissection (starting in the distal aortic arch). Type A dissection normally requires immediate surgery, whereas for Type B the decision for medical treatment or surgical intervention depends on the clinical presentation. In these cases, detailed analysis of flow patterns in the dissection aorta can provide useful information for clinical decision-making.
∗ Corresponding author. E-mail address: [email protected] (X.Y. Xu). http://dx.doi.org/10.1016/j.medengphy.2014.07.006 1350-4533/© 2014 IPEM. Published by Elsevier Ltd. All rights reserved.
The development of high-resolution medical imaging and robust image processing techniques has promoted image-based computational fluid dynamics (CFD) modelling of in vivo haemodynamics [1]. A number of studies have been carried out to investigate the flow patterns and haemodynamic environment in dissected aorta. Karmonik et al. [2] examined changes in anatomy as well as pressure and wall shear stress in an aortic dissection case with aneurysmal dilatation. In a longitudinal study, Tse et al. [3] investigated variations of certain haemodynamic parameters in the development of a dissecting aneurismal aorta. In our previous studies, an advanced transitional turbulence model was employed to investigate the complex flow patterns and detailed haemodynamic environment in multiple cases of type B dissection [4,5]. Results were analysed in each case and compared between different groups. One of the most important findings was that flow rate into the false lumen was strongly correlated with the primary entry tear size. Since false lumen patency is a key predictor of adverse early and later outcomes [6], false lumen flow rate (or entry
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1177
Fig. 1. Reconstructed geometries of the aortic dissection used in validation. (a) The entire aorta including the ascending aorta lumen (AL), the aortic arch branches, the descending aorta lumen (DL), and the iliac arteries (IA), but excluding the branches off the descending thoracic and abdominal aorta; (b) the locations of three imaging planes where phase-contrast MR velocity mapping was acquired; and (c) the computational model used in validation, with arrows indicating the locations of inlet, outlet, tear, false, and true lumen.
tear size) may be an important parameter for risk stratification of dissection patients. However, computational models for the analysis of flow in aortic dissections (e.g. [2–5]) often involve several assumptions, including compensated exclusion of the aortic arch branches, flat velocity profile at the inlet, and rigid walls accompanied by fixed, constant distal outflow pressures. The purpose of this study is to assess the computational models adopted in our previous studies by comparison with in vivo velocity data obtained by means of phase-contrast magnetic resonance imaging (PC-MRI). Computed tomography (CT) is the standard imaging modality for diagnostic imaging of aortic dissection, but for flow simulations, 4D ciné PC-MRI would be desirable to provide patient-specific flow boundary conditions [7]. PC-MRI flow imaging is considered the only imaging technique that can provide measurement of blood velocities in vivo with good resolution both spatially and temporally. In this study, a typical type B aortic dissection case was analysed by using the same computational model design and numerical analysis methods as adopted in our previous studies [4,5]. This case was chosen for its comprehensive set of PC-MRI data which provided not only the necessary boundary conditions, but also additional data required for validation. 2. Methods 2.1. Medical image processing and geometry reconstruction Anatomical information was extracted from CT scan of a patient (male, age 56) diagnosed with a typical type B aortic dissection. Subsequently, PC-MRI data were acquired for further examination. The study complied with the Declaration of Helsinki and was approved by the local Research Ethics Committee. The patient gave written, informed consent. The model geometry was based on the initial CT images covering the entire aorta (approximately 300 slices, 2 mm thick, spacing 1 mm; pixel resolution 0.86 mm), as shown in Fig. 1a. The proximal entry tear can be seen clearly in the transparent image (Fig. 1b). Anatomical MR images had relatively low axial resolution and were only used for registration and localization of imaging planes
between the CT and PC-MRI data. The CT images were imported into an image-processing package Mimics (Materialise HQ, Louvain, Belgium) and the geometry of the dissected aorta was reconstructed from the contours of the segmented lumen area from the serial CT slices and was smoothed accordingly before meshing. PC-MRI was performed at three planes in the aortic dissection (Fig. 1b), providing sufficient data for inlet boundary condition and validation (TE 2.6 ms, TR 4.7 ms, flip angle 15◦ , pixel spacing 1.3281 mm, slice thickness 8 mm). At the imaging planes 1 and 3, encoding velocity (Venc ) was 200 mm/s for IS (inferior–superior) direction and 100 mm/s for AP (anterior–posterior) and RL (right–left) directions. At imaging plane 2, Venc was 200 mm/s for AP direction, and 100 mm/s for the other two directions. Plane 1 was located in the proximal ascending aorta in order to obtain velocity profiles distal to the aortic valve. Plane 2 was in the aortic arch between the root of the left subclavian artery and the site of proximal entry tear. Since the tear occurred very close to the left subclavian artery, there was limited space to fit plane 2 between the artery and the entry tear. Flow imaging at plane 2 was required to determine the amount of flow distal to the arch branches, and velocity profiles proximal to the tear. Plane 3 was in the mid-thoracic descending aorta distal to the proximal entry tear, at the level of pulmonary bifurcation, covering both the true and false lumen. The PC-MRI data from this plane were used to compare with those from the computational simulation for model validation. Images acquired with PC-MRI were processed by using an in-house programme. The region of interest was segmented from magnitude images and then mapped on to the corresponding phase images; the grayscale of the pixels inside the region of interest was converted into velocity information according to the applied encoding sequence of the scanner (Philips Medical Systems). For imaging planes 1 and 2, the flow rate through each plane was calculated based on the converted velocity and pixel size; this was utilized as boundary conditions for the flow model described in the next section. The PC-MRI data acquired at plane 3 were processed in the same way, and 2D velocity contours of the true and false lumen were plotted for the purpose of comparison with the simulation results.
1178
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
The computational model (Fig. 1c) shows the inlet at imaging plane 1. To keep the approach consistent with previous patientspecific studies, the model was simplified by excluding the three branches in the aortic arch and other small branches in the descending aorta, whilst matching the flow at planes 1 and 2 to provide the correct amount of flow through planes 2 and 3. The false lumen propagated down the lateral side of the thoracic descending aorta, before spiralling around the true lumen toward its anterior side at the mid-descending aorta. It continued to twist to the rightposterior side of the aorta, then to the lateral side in the distal segment of the abdominal aorta. As it extended to the left iliac artery below the abdominal bifurcation, two adjacent outlets (one for each lumen) were included in the model, both located in the abdominal aorta far distal to the imaging plane 3 (Fig. 1c). The length of the descending aorta was approximately 381 mm. The proximal entry tear was about 22 mm in width and 25 mm in height, and only 0.4 mm below the top of the aortic arch. The false lumen was prevented from propagating proximally by the left subclavian artery. No distal tears could be identified. 2.2. Flow simulation and boundary conditions Simulations were carried out with ANSYS CFX 12, and the – correlation-based laminar–turbulent transition model (‘SST Tran’) was adopted [8], as described previously [4]. Blood was assumed to be Newtonian with viscosity 4.0 mPa s and density 1060 kg/m3 . The aortic walls were assumed to be rigid with no-slip. Flow boundary conditions were extracted from the PC-MRI velocity images. The flow data at each imaging plane were recorded at 29 time points along a cardiac cycle of 0.98 s duration. The flow rate waveforms at imaging plane 1 (ascending aorta root) and imaging plane 2 (aortic arch after the branches) are presented in Fig. 2. The peak systolic flow rate was 446 cm3 /s in the ascending aorta, and 370 cm3 /s in the distal aortic arch. Hence, the flow exiting through the three main branches on the arch at peak systole was 76 cm3 /s, which was 17% of the peak flow rate from the heart. This ratio is close to the value (15%) assumed previously [4,5]. The flow waveform derived from imaging plane 2 was specified at the model inlet in the ascending aorta with a flat through-plane profile, as adopted in previous patient-specific studies. Based on the inlet area, the mean and maximum Reynolds numbers were 850 and 3512, respectively, and the Womersley parameter was 10.4. A low inlet turbulence level of 1.5% was chosen to represent initial disturbances in the flow and allow true transition to occur realistically. At the model outlets, zero relative static pressure was applied across both surfaces. A uniform time step of 0.001 s was used. All simulations were carried out for three cardiac cycles to achieve a
Fig. 2. The flow rate waveforms at imaging plane 1 and plane 2 as defined in Fig. 1b. The flow rate was calculated based on the converted velocity at each pixel and the pixel size of the PC-MRI data.
periodic solution, and the results of the last cycle were used for comparison. 2.3. Mesh sensitivity test The reconstructed model was imported into ANSYS ICEM CFX for mesh generation. Unstructured mesh consisting of 3D tetrahedral cells combined with prismatic cells near the wall was created. A very fine near wall resolution in the tear region was ensured by using a minimum of 10 layers of prismatic cells at the wall, height of wall cells in wall units (y+) < 2, together with local refinement. Mesh independence tests were carried out with an initial mesh consisting of 2,600,000 elements and finer meshes with increased number of prismatic wall elements. The results in terms of maximum timeaveraged wall shear stress (WSS) and turbulence kinetic energy were compared between different meshes. The differences in these parameters were less than 4% between the initial 2,600,000 elements mesh and the 3,200,000 elements mesh, and less than 3.5% between the 3,200,000 elements mesh and the 4,700,000 elements mesh. Hence, the mesh of 3,200,000 elements was deemed to be sufficient. 2.4. Analysis of simulation results The simulation results were analysed using ANSYS CFX-Post 12. For the purpose of comparison, velocities are displayed by plotting colour-coded contours at the cross-sectional plane corresponding to imaging plane 3 where PC-MRI data were acquired. Instantaneous streamlines are presented at different time points for better visualisation of time-varying flow patterns. In order to examine the turbulence level in the model, turbulence intensity (Tu) was evaluated, which is defined as:
Tu =
2 k/V 3
(1)
where k is the turbulence kinetic energy and V is the instantaneous local velocity. The iso-surface of Tu is presented to indicate the corresponding turbulence zone with specified value at the systolic mid-deceleration phase, when Tu is usually at its highest value. 3. Results 3.1. Flow distribution in the true and false lumen As reported by Cheng et al. [5], the flow division between the two lumens was strongly dependent on the size and location of the proximal entry tear. Fig. 3 shows the comparison of timevarying flow rate at the location of imaging plane 3 obtained from
Fig. 3. The flow rate waveforms of the true and false lumen from both CFD results and PC-MRI velocity data at imaging plane 3.
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1179
Fig. 4. Comparison of velocity profiles in FH, AP, and RL directions at mid-systolic acceleration phase of the cardiac cycle between PC-MRI velocity data (in m/s) and simulation results. In Figs. 4–7 the red colour represents positive velocities from superior to inferior, anterior to posterior (IS contours), right to left (AP and RL contours); the blue colour represents negative velocities. Note different colour scales for different directions and time points according to the range of velocity values. The secondary velocity vectors are also presented as project on the plane.
the CFD simulation and the corresponding PC-MRI data. Based on the simulation results, the peak flow rate in the false lumen was 293 cm3 /s, and 72 cm3 /s in the true lumen, whereas the cycleaveraged values were 101.6 and 27.1 cm3 /s, respectively. Hence, about 80% of the peak flow was diverted to the false lumen. As shown in Fig. 3, the measured and computed flow rate waveforms are similar in shape but the CFD model over-predicted flow in the false lumen and under-predicted flow in the true lumen. The differences between the computed and measured flow rates were 8.5% for the true lumen and 25% for the false lumen at peak flow rate (the cycle-averaged values were 35% for the true lumen and 28% for the false lumen). A number of factors could have contributed to the differences; delineation of the lumen contours on PC-MRI images; motion of the intra-arterial septum; limited spatial resolution, noise, and apparent artefacts associated with the PC-MRI velocity data. 3.2. Comparison of flow structure and velocity profiles Comparisons were made between the time-varying velocity contours from PC-MRI measurements and flow simulation results at plane 3 (defined in Fig. 1b) for all three velocity components. Figs. 4–7 present the velocity contours in the main inferior–superior (labelled FH, foot-head), secondary AP (anterior–posterior), and RL (right–left) directions extracted from
the PC-MRI data and simulations at four time points along the cardiac cycle. Both sets of data are displayed in the conventional MRI direction, looking towards the head. The corresponding magnitude images from the PC-MRI are included, indicating the shape and size of the false (FL) and true lumen (TL). As can be seen, the FL is approximately four times larger than the TL at this location. For better understanding of the complex flow patterns, secondary velocity vectors have been superimposed on the FH velocity contours. Complementary instantaneous streamlines for the four time points are shown in Fig. 8. At mid-systolic acceleration, the velocity contours (Fig. 4) show that the main direction of the flow in the TL was from superior to inferior (distal) with a maximum velocity of 0.45 m/s based on PCMRI. The high velocity flow was skewed towards the left-anterior side of the TL owing to the curvature of the arch; a similar pattern can be seen in the simulation result. The flow in the FL showed more spatial variation than in the TL. The simulation results captured the relatively high velocity at the left anterior side near the septum and in the middle of the FL. The corresponding secondary velocity vectors show that in the TL the direction of secondary flow was predominantly anterior–posterior, pointing to the right. Similar patterns can also be observed in the AP and RL contours. The proportion of the flow remaining in the TL was accelerated by the narrowing of the lumen, so stretching the helical pattern induced by the non-planar arch. In the FL the pattern was more complex,
1180
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
Fig. 5. Comparison of velocity profiles (in m/s) between PC-MRI velocity data and simulation results in FH, AP, and RL directions at peak systole.
with a weak vortex centred near the right (R) wall, but extending across most of the area. The vortex originated from the flow entering at the tear proximally to this plane, curving along the opposite wall, setting up a large proximal vortex (Fig. 8a), which expanded to fill the lumen by half-way along the descending thoracic aorta. At peak systole (Fig. 5), the PC-MRI flow patterns in the TL showed fairly uniform high velocities at this location, which were well captured by the CFD simulation. The corresponding streamlines and secondary velocity vectors again reflected the deviation of the TL axis from the FH direction. In the FL the pattern had changed in response to the higher velocities through the tear, shown clearly in Fig. 8b. The vortex previously centred near the right wall had moved anteriorly and the effect of the stronger jet deflected along the left wall could be seen (recall flow is shown from the inferior aspect, as with the MRI convention). The high velocities along the left wall resulted from the jet through the tear as it impacted the opposite wall and followed the curvature of the arch; at plane 3 (Fig. 8b) the jet had spread to fill the cross-section, with lower velocities towards the right-anterior wall resulting from the greater expansion of the jet distal to the separated, recirculating flow immediately below the tear. At mid-systolic deceleration (Fig. 6), flow in the TL was still dominated by the superior–inferior component, and the velocity distribution was similar to that at peak systole, with a good agreement between measured and simulated velocities. There was a subtle change in flow near the anterior wall as a result of the helical flow in the aortic arch induced by flow deceleration, as seen in
the instantaneous streamlines (Fig. 8c). The effect of the retardation was clearer in the FL. The flow distal to the weaker jet through the tear was more disorganised, as shown in Fig. 8c, whilst in the region of Plane 3 the somewhat disorganised flow became more evident. In Fig. 6, the simulation showed the high velocity zone having moved farther along the right wall than the measurements, while the secondary vectors revealed that the main vortex had moved towards the centre of the area and a weak secondary vortex had formed at the left-posterior wall. The simulation predicted significant turbulence in the FL distal to the tear, particularly in late systole (Fig. 9). Comparisons of velocities in early diastole, when the maximum retrograde flow occurred in the proximal ascending aorta, are given in Fig. 7. All velocity contours from MR and CFD showed excellent agreement in both lumens. The TL was dominated by a helical flow in the superior–inferior direction, despite the area mean flow in the ascending aorta being retrograde. The corresponding streamlines and secondary vectors revealed a complex flow structure in the arch (Fig. 8d). Flow in the FL at Plane 3 showed retrograde flow on the left-anterior area, occupying two-thirds of the crosssectional area, whilst it was antegrade in the remaining third, at the right-posterior aspect (Fig. 7). The streamlines reflected this feature, where it would appear from Fig. 8d that part of the flow changed direction approximately one diameter distal to Plane 3. Quantitative comparisons were made for the maximum velocities in the TL and FL between the PC-MRI data and CFD results. The maximum difference of instantaneous velocities in both TL and
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1181
Fig. 6. Comparison of velocity profiles (in m/s) between PC-MRI velocity data and simulation results in FH, AP, and RL directions at systolic deceleration phase.
FL was found at the systolic mid-deceleration phase, when the TL velocity was 29% higher with the PC-MRI measurement than with the CFD prediction, whilst the FL velocity was 30% higher with the CFD than with the PC-MRI. For the overall time-averaged comparison, the equivalent differences were 20 and 16%. Reasons for the over-prediction of FL velocity have been explained in the preceding section. 4. Discussion There have been a few limited examples of such studies in the literature; Tse et al. [3] examined a dissected aorta pre- and postaneurysm, whilst Karmonik et al. [2] reported on a chronic type III dissection both at initial examination and at 10-month followup. Neither of the studies incorporated a transitional or turbulence model in their CFD simulations. Based on the stability diagram adopted by Kousera et al. [9], the combination of peak Reynolds number and Womersley parameter considered here is likely in the region of disturbed flow, justifying the need for a transitional flow model. With regard to boundary conditions. Karmonik et al. [2] used patient-specific flow waveform measured with PC-MRI as the inlet boundary condition, while Tse et al. [3] applied a typical aortic flow waveform at the inlet and a corresponding pressure waveform at the exits in the iliac arteries. Although Karmonik et al. [2] showed comparison of PC-MRI and CFD predicted velocities at a transverse plane in the descending aorta, it was for the axial through-plane velocity magnitudes only. Detailed assessment of CFD predictions using MR velocity imaging in three orthogonal directions is lacking. False lumen flow rate is a key predictor of FL patency, being the most important predictor of adverse early and later outcomes
[6]. An interesting finding of Cheng et al. [5] was the correlation between the FL flow rate and the size and location of the primary entry tear, based on the computational results of four patient-specific models with the same assumptions and simplified boundary conditions as employed here. Comparison of the FL flow rate obtained from the CFD simulation and that from the PC-MRI data (Fig. 3) showed that the computational model gave comparable results to the MR measurements with similar flow waveforms for each lumen, and the percentage of flow in the FL (80% vs. 76%). The size and location of the primary entry tear and the percentage of FL flow rate extracted from the MR data were found to fit well to the regression lines presented by Cheng et al. [5]. Further careful comparison was made between the 4D PC-MRI velocity data and the simulation results at the proximal middescending aorta. The velocity data were compared in each of the three components. It was found that (Figs. 4–7) velocity contours derived from the CFD simulation were in good agreement with those from PC-MRI. No retrograde flow was found in the TL, as the dominant velocity component (in the IS or FH direction) was directed from head to feet throughout the cardiac cycle. The other two velocity components also showed consistent patterns in pointing towards the posterior and right directions, respectively, caused mainly by the inclination of the TL to the FH direction. The organized flow pattern in the TL, which was also demonstrated in previous studies [2–5], was related to its compressed size and relatively straight geometry. Flow patterns in the FL were much more complicated, seen both in the contours (Figs. 4–7) and related instantaneous streamlines (Fig. 8). The main flow was generally in the superior–inferior direction and the flow velocity was lower in the other two directions.
1182
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
Fig. 7. Comparison of velocity profiles (in m/s) between PC-MRI velocity data and simulation results in FH, AP, and RL directions at maximum retrograde flow in early diastole.
At mid-systolic acceleration, no clear velocity pattern can be observed in the FH direction but CFD results successfully captured the patterns in the AP and RL directions. By peak systole, the FH velocity was higher near the left wall where CFD and PC-MRI data showed the same pattern. However, the other two directions showed weaker agreement. One likely reason for this was the patchy nature of the PC-MRI data which was especially amplified when velocity was low. Another key reason might be the rigid wall assumption in the CFD simulations. The movement of the dissection flap is likely to be most significant at peak systole when the pressure difference between the TL and FL is at the highest [10]. Although the intimal flap was found to have been calcified and quite stiff in the particular case concerned here, small motions could still happen, hence affect the flow velocity near the flap especially in the AP and RL directions. During retardation (Fig. 6), the agreement between simulations and PC-MRI was good for all velocity components. The TL flow remained antegrade, whilst there was a substantial retrograde zone in the FL at plane 3. It could be discerned from the streamlines in Fig. 8c that flow was largely disturbed and recirculating. There was a good agreement between PC-MRI and simulations at the maximal retrograde flow in early diastole as seen in the velocity contours of Fig. 7. Considerably disorganised flow in the FL was also evident in Fig. 8d. In the region of the tear, there was major recirculation but, because the septum does not lie perpendicular to the oblique sagittal plane, it remains difficult to discern the boundary between FL and TL.
In summary, the computational model successfully captured the complex regions of flow reversal and recirculation qualitatively, although quantitative differences exist. The overall good degree of agreement confirmed the validity of the computational model utilized in our previous case studies, providing reliable data on flow patterns and haemodynamic parameters.
4.1. Limitations The CFD model presented here involved several assumptions. First, the aortic wall and dissection flap were assumed to be rigid. While this assumption was acceptable in the case discussed here owing to its stiff flap, this might not be true in cases when the flap is more compliant. The dissection flap is initially very mobile immediately after the dissection has occurred, but becomes more rigid over time. Ganten et al. [11] examined the intimal flap and wall motion in 32 type B aortic dissection patients by using dynamic, ECG-gated CT imaging technique. It was found that distensibility of the dissected descending aorta was significantly lower than that in healthy subjects of the same age, and the average true lumen narrowing by intimal flap motion was 4.4% (maximum 29%). However, in 14 of the 32 patients there was no detectable flap motion. Another study based on a smaller number of cases demonstrated the highly variable flap motion in type B dissections [10], with up to 30% change in true lumen area. Using 2D PC-MRI in a longitudinal
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
1183
Fig. 8. Instantaneous streamlines of blood flow from left posterior aspect at (a) systolic mid-acceleration phase, (b) peak systole, (c) systolic mid-deceleration phase, and (d) maximum retrograde flow in diastole.
Fig. 9. Turbulence intensity isosurfaces at mid-systolic deceleration from left-posterior (left) and right-anterior (right) aspects.
1184
Z. Cheng et al. / Medical Engineering & Physics 36 (2014) 1176–1184
study, Karmonik et al. [12] found that flap motion in type B aortic dissections reduced over time. The exclusion of the arch branches is another assumption in this study. To compensate for the flow loss through the three branches, the inlet flow rate at the aortic root was scaled down to maintain the correct amount of flow into the descending aorta. This appeared to give satisfactory results, although the flow structure in the entry tear region would have been compromised. However, owing to the nature of the geometry of dissected aorta, the flow pattern is significantly changed after passing through the tear into the false lumen, as well as in the highly compressed true lumen. In conclusion, we have demonstrated good agreement in flow structure and velocity contours in the dissected aorta, proving the validity of the computational model and the simulation approach for type B aortic dissection with a relatively stiff dissection flap. Future studies are needed to examine the effect of flap motion on flow in the true and false lumen. Declarations Funding: Zhuo Cheng was supported by a PhD scholarship from the Institute of Biomedical Engineering and the National Heart and Lung Institute, Imperial College London. Competing Interests: There is no conflict of interest that could inappropriately influence this research work. Ethical Approval: Ethical approval was given by the local Research Ethics Committee at the St Mary’s Hospital. References [1] Steinman DA. Image-based computational fluid dynamics modeling in realistic arterial geometries. Ann Biomed Eng 2002;30(4):483–97.
[2] Karmonik C, Partovi S, Müller-Eschner M, Bismuth J, Davies MG, Shah DJ, Loebe M, Böckler D, Lumsden AB, von Tengg-Kobligk H. Longitudinal computational fluid dynamics study of aneurysmal dilatation in a chronic DeBakey type III aortic dissection. J Vasc Surg 2012;56:260–3. [3] Tse KM, Chiu P, Lee HP, Ho P. Investigation of hemodynamics in the development of dissecting aneurysm within patient-specific dissecting aneurismal aortas using computational fluid dynamics (CFD) simulations. J Biomech 2011;44:827–36. [4] Cheng Z, Tan FPP, Riga CV, Bicknell CD, Hamady MS, Gibbs RGJ, Wood NB, Xu XY. Analysis of flow patterns in a patient-specific aortic dissection model. J Biomech Eng 2010, 132051007-1-9. [5] Cheng Z, Riga CV, Chan J, Hamady MS, Wood NB, Cheshire NJW, Xu XY, Gibbs RGJ. Computational simulation of the aorta in acute type B dissection: initial findings and potential applicability. J Vasc Surg 2013;57: 35S–43S. [6] Eggebrecht H, Herold U, Kuhnt O, Schmermund A, Bartel T, Martini S, Lind A, Naber CK, Kienbaum P, Kuhl H, Peters J, Jakob H, Erbel R, Baumgart D. Endovascular stent-graft treatment of aortic dissection: determinants of postinterventional outcome. Eur Heart J 2005;26:489–97. [7] Clough RE, Waltham M, Giese D, Taylor PR, Schaeffter T. A new imaging method for assessment of aortic dissection using four-dimensional phase contrast magnetic resonance imaging. J Vasc Surg 2012;55:914–23. [8] Menter F, Langtry R, Volker S. Transition modelling for general purpose CFD codes. Flow Turbulence Combust 2006;77:277–303. [9] Kousera CA, Wood NB, Seed WA, Torii R, O’Regan D, Xu XY. A numerical study of aortic flow stability and comparison with in vivo flow measurements. J Biomech Eng 2013;135, 011003-1-9. [10] Murayyama T, Funabashi N, Uehara M, Takaoka H, Komuro I. New classification of aortic dissection during the cardiac cycle as pulsating type and static type evaluated by electrocardiogram-gated multislice CT. Int J Cardiol 2010;142:177–86. [11] Ganten MK, Weber TF, von Tengg-Kobligk H, Bockler D, Stiller W, Geisbusch P, Kauffmann GW, Delorme S, Bock M, Kauczor HU. Motion characterization of aortic wall and intimal flap by ECG-gated CT in patients with chronic Bdissection. Eur J Radiol 2009;72:146–53. [12] Karmonik C, Duran C, Shah DJ, Anaya-Ayala JE, Davies MG, Lumsden AB, Bismuth J. Preliminary findings in quantification of changes in septal motion during follow-up of type B aortic dissections. J Vasc Surg 2012;55: 1419–26.