 Research
 Open Access
 Published:
Assessment of turbulent blood flow and wall shear stress in aortic coarctation using imagebased simulations
BioMedical Engineering OnLine volume 20, Article number: 84 (2021)
Abstract
In this study, we analyzed turbulent flows through a phantom (a 180\(^{\circ }\) bend with narrowing) at peak systole and a patientspecific coarctation of the aorta (CoA), with a pulsating flow, using magnetic resonance imaging (MRI) and computational fluid dynamics (CFD). For MRI, a 4Dflow MRI is performed using a 3T scanner. For CFD, the standard \(k\epsilon \), shear stress transport \(k\omega \), and Reynolds stress (RSM) models are applied. A good agreement between measured and simulated velocity is obtained for the phantom, especially for CFD with RSM. The wall shear stress (WSS) shows significant differences between CFD and MRI in absolute values, due to the limited nearwall resolution of MRI. However, normalized WSS shows qualitatively very similar distributions of the local values between MRI and CFD. Finally, a direct comparison between in vivo 4Dflow MRI and CFD with the RSM turbulence model is performed in the CoA. MRI can properly identify regions with locally elevated or suppressed WSS. If the exact values of the WSS are necessary, CFD is the preferred method. For future applications, we recommend the use of the combined MRI/CFD method for analysis and evaluation of the local flow patterns and WSS in the aorta.
Background
Coarctation of aorta (CoA) is a congenital condition in which the aorta has a narrowing, usually in the thoracic descending aorta distal to the branching arteries of the aortic arch. The narrowing of the artery causes flow acceleration, where a turbulentlike flow may occur during the systolic phase [1]. It has been shown that the transitional and turbulent flow in CoA leads to aberrant blood flow in the narrowing and a vortexlike recirculation pattern distal to the stenosis [2]. Due to the stenosis and onset of turbulence, the wall shear stress (WSS) is also elevated, and the presence of turbulence may cause oscillations of its values [3]. This type of flow may cause, among others, degradation of the arterial wall, initialization of an aneurysm, and atherosclerosis [4].
Several studies have assessed the hemodynamics of this pathology using magnetic resonance imaging (MRI) [5, 6]. However, due to the relatively low spatial resolution of MRI, the flow velocity in the proximity of the wall and its derived quantities such as WSS may be incorrect [7]. Several recent studies have identified imagebased computational fluid dynamics (CFD) to be a good alternative to study blood flow in CoA [8, 9]. However, these previous studies have not addressed the important effects of locally generated turbulence in CoA, as demonstrated by Gaze et. al. [10]. The turbulent flow in CoA has been simulated by using various turbulence modeling approaches: (i) the Reynoldsaveraged Navier–Stokes (RANS) ([11]; (ii) largeeddy simulation (LES) [12]; and finally, (iii) direct numerical simulation (DNS) [13] methods. The DNS and LES proved to perform very well for transient and turbulent flow regimes for arteries with stenotic regions [14]. However, because of the huge computational costs associated with high temporal and spatial resolution requirements of LES and DNS, these approaches are less suitable for clinical applications [15]. To meet demands on a computationally efficient and sufficiently accurate CFD approach, we propose to employ the unsteady RANS method with an advanced secondorder momentsbased turbulence model (socalled Reynolds stress model, RSM). The advantage of this model lies in its ability to automatically take into account exact production terms of the turbulent stresses (which need to be additionally modeled in the eddyviscosity type of RANS model), as well as to predict turbulence anisotropy (in contrast to the assumption of turbulence isotropy as used in the eddyviscosity turbulence models), which are important features of a flow in turbulent regime.
In the present study, we will first introduce a Ubend phantom that mimics the aorta with coarctation and produces numerous flow features observed in vivo. For the phantom, we will perform a detailed comparison between experiments (performed by 4Dflow MRI) and CFD simulations. In addition to the proposed RSM turbulence model, also two widely used eddyviscositybased turbulence models will be introduced. We will investigate levels of agreement between CFD and phantom experiments by focusing on the mean flow features and local distributions of the wall shear stress. Finally, we will perform a comparative assessment between CFD simulation (based on the turbulence model which performed best for the phantom study) and in vivo 4Dflow MRI for the patientspecific pulsating blood flow in CoA.
Results
Phantom
4Dflow MRI flow rate
The volumetric flow rate was extracted from 11 different locations distributed evenly along the length of the phantom using CAAS MR Solutions v5.0 to test the performance of the MRI acquisition. We calculated the error between the set inlet volumetric flow \(\hbox {Q}_{0}\) and the flow extracted at the different cut planes \(\hbox {Q}_{i}\). The average error in flow rate was \(0.25 \pm 2.11\%\).
Comparison of turbulence models
Next, we moved towards a detailed comparison of the MRI and CFD simulations performed with various turbulence models by comparing the nondimensional velocity (\(v/v_0\), where \(v_0\) is the mean inlet velocity) magnitude profiles at six characteristic locations: inlet (A), the start of the bend (B), middle of bend (C), end of the bend (D), middle of narrowing (E), and distal to narrowing (F), as shown in Fig. 1. The differences between the models emerge in the middle of the bend, i.e., at location C. This location is particularly sensitive due to the generation of the secondary flows (Dean vortices) and flow acceleration along the outer wall curvature. Additionally, at the poststenotic location (location F), numerical simulations captured well the recirculation region, however, the \(k\varepsilon \) model underestimated the velocity magnitude in the center, whereas both SST and RSM models are showing a very good agreement with MRI in the wall vicinity, with a slight overprediction in the center.
The observed changes in the distributions of the velocity profiles for CFD with considered turbulence models are due to different predictions of turbulence levels. To illustrate this, we plot series of the profiles of the nondimensional turbulent kinetic energy (\({\mathrm{k/v}}_0^2\)) at identical locations as previously analyzed for the velocity profiles (locations A–F), Fig. 1g–l. All turbulence models are giving similar profiles at inlet segment location (A), with almost identical values in the center and symmetrical peak values in the proximity of the wall and are in good agreement with the previously reported results [16]. The symmetrical distribution is, with elevated turbulence levels in the proximity of the outer wall, due to the presence of the bent (B–D). It is interesting to note that in the center of stenosis (location E), despite a big overprediction of turbulent kinetic energy by the \(k\varepsilon \) model, resulting velocity magnitudes still agree very well due to the dominance of convective term in the momentum equation (due to a sudden flow acceleration). Finally, it can be seen that the poststenotic region (location F) is characterized by the highest levels of turbulence caused by combined effects of a flow acceleration (in the center) and flow recirculation (in the wall proximity).
Voxeltovoxel velocity and vorticity comparisons
As the next step, we will compare in more detail (voxeltovoxel) results of CFD (with the best performing turbulence model, RSM) against MRI. The contours of the velocity magnitude in the central horizontal crosssection (\(\hbox {y} = 0\)) are shown in Fig. 2. Here, we present the velocity magnitude distribution on original CFD resolution (CFD), downsized CFD resolution (DCFD, where downsizing is done to match original MRI resolution), original measurements (MRI), and the absolute difference between downsized CFD and MRI (DCFDMRI), respectively. It can be seen that an overall good agreement is obtained between DCFD and MRI and that all most salient flow features are well captured with both techniques. The small deviations are located in the proximity of the walls (in the curved part) and central stenotic and poststenotic regions.
To compare secondary flow patterns, we plot contours of the outofplane vorticity component for all cases at characteristic selected crosssections (A–F), Fig. 3. Note that the outofplane vorticity component (defined in here adopted coordinate system as: \(\omega _x=\partial w/\partial y  \partial v/\partial z\)) is a sensitive flow parameter since it captures gradients of both velocity components in the particular plane perpendicular to the flow direction. At the inlet segment location (A), there should not be yet any significant appearance of the secondary motions, as confirmed by CFD and DCFD results. The MRI contours show a more noisy distribution, but its levels are relatively small. By entering the bend curvature (location B), the vorticity starts to be generated. Here, due to limited spatial resolution, the MRI just partially captures some of secondary flow features. The agreement is much better at the most interesting location in the center of the bend (location C) where all cases captured welldetailed structure of the Dean vortices. The traces of Dean vortices are still visible at the end of the curved bend (location D), where a satisfactory agreement between CFD and MRI is obtained. A similar level of agreement is also obtained in the stenosis center (location E). Some visible deviations are observed in the poststenotic region (location F), where the differences in vorticity magnitude are more pronounced.
Wall shear stress
The contours of the wall shear stress along the phantom walls for the original CFD (obtained with RSM turbulence model), downsized CFD results (\(\hbox {DCFD}_{0.2\times 0.2\times 0.2}\) and \(\hbox {DCFD}_{0.7\times 0.7\times 1.5\,\hbox {mm}^3}\)), and original MRI are shown in Fig. 4a. To provide a complete distribution of the WSS along the phantom wall, we generated twodimensional maps of WSS, where the horizontal coordinate represents the nondimensional circumference (expressed in angles, \(\pi \le r \le +\pi \), where \(r=0\) indicates the innercurve and \(r=\pm \pi \) indicate the outer phantom curve), while the vertical coordinate represents the nondimensional enveloped arclength of the phantom (\(0\le l/l_0 \le 1\), where 0 and 1 correspond to the start and the end of the phantom, and \(l_0\) is the centerline length), Fig. 4b. The nondimensional WSS maps (WSS/\(\hbox {WSS}_{\mathrm{mean}}\), where \(\hbox {WSS}_{\mathrm{mean}}\) indicates the spatially averaged WSS over the entire phantom surface, given in Table 1) are shown in Fig. 4c. Finally, profiles of the circumferentially averaged nondimensional WSS are shown in Fig. 5.
The agreement between CFD and MRI is good in the inlet leg of the phantom. The differences emerge in the bent and stenotic regions. Here, with the decreasing resolution, the absolute values of WSS decrease. It can be seen that all cases are predicting high values of WSS in the stenotic region (within the dashed lines), but lowering of the spatial resolution of the CFD results and MRI produced a slight shift of the location where WSS reached its maximum value.
Patientspecific CoA: in vivo MRI and CFD based on RSM turbulence model
After demonstrating that CFD with RSM turbulence model was sufficient to predict the characteristic flow features in the phantom, we next moved to the patientspecific aorta geometry with coarctation, for which in vivo 4Dflow MRI measurements are available, Fig. 7c. The resulting flow pattern, presented in form of streamtraces colored by the velocity magnitude at the peak systole is shown in Fig. 6a. It can be seen that a good agreement between MRI and CFD is obtained in capturing important flow features: a strong helical pattern in the aortic arch, and a sudden flow acceleration in the coarctation. A summary of direct comparison between CFD and MRI in predicting the peak and spatially averaged (mean) values of the WSS is provided in Table 2. It can be seen that in vivo MRI underestimated the local values of WSS, similarly to our previous findings in the phantom geometry. The effect of the surface smoothing revealed relatively small differences between the rough and smoothed geometries. Note that present CFD results agree well with similar numerical studies reported in the literature, e.g., [17,18,19] Instead of focusing on the local differences in WSS from MRI and CFD in their absolute terms, we proceed with qualitative comparisons between simulations and experiments by identifying regions along aorta walls characterized by locally elevated or lowered values of WSS to their spatially averaged mean values (averaged over entire aorta wall). The contours of the nondimensional WSS distribution (WSS/\(\hbox {WSS}_{\mathrm{mean}}\)) for CFD (both rough and smoothed geometry) and MRI are given in Fig. 6b. It can be seen that an overall good agreement is obtained, especially when considering the smoothed CFD and MRI distributions in the coarctation and the descending part of aorta. This is additionally illustrated by showing 2D maps of the local nondimensional WSS, where the entire surface of aorta wall is mapped, Fig. 6c, d. The blank spaces in the mapped surfaces represent the branching arteries that were removed during the mapping procedure. Due to the proximity of the first two branching arteries (i.e., brachiocephalic trunk and left common carotid artery), they are merged on the mapped surface.
Finally, the scatter plots (symbols) and circumferentially averaged nondimensional WSS profiles (lines) are shown in Fig. 6e. Similar to comparisons in the phantom geometry, qualitatively good agreement is obtained with distinct peak values in the coarctation. A shift in the location of the maximal WSS for the MRI is also observed. Note that larger peaks of WSS from CFD at \(l/l_0=0.9\) locations are due to the secondary side branches of the aorta which are not properly resolved in MRI.
Discussion
The present study investigated the flow and WSS for a specific aorta pathology—aortic coarctation—using a combined 4Dflow MRI and CFD techniques for simplified (phantom) and patientspecific geometry. For both studied geometries, the flow is in a fully (phantom) or a partially developed (CoA) turbulent regime. The importance of selected turbulence models for CFD is demonstrated by performing a comparative assessment between two types of eddyviscosity models and RSM for the phantom configuration.
Comparison of turbulence models with MRI
The differences in the crosssection averaged velocity magnitude profiles (at A–F locations) between MRI and CFD with the RSM model did not exceed 15% in the stenotic region and 7% in the rest of the phantom. In contrast, the maximum disagreement between CFD with \(k\varepsilon \) model and MRI was 45% for the poststenotic region, while this disagreement was reaching 18% in the stenotic region for the SST model. The overall best performances of the RSM turbulence models can be explained in terms of the theoretical foundation behind this model. The exact treatment of production terms of the individual turbulent stress components plays a crucial importance in complex threedimensional flows (e.g., curved part of the phantom followed by a stenotic region) as presented here. In comparison with the eddyviscosity models, the RSM predicts well the secondary motions and captures well flow adaptation to sudden changes of the crosssectional area [20]. This was also shown in terms of streamwise velocity and turbulent kinetic energy in a \(60^{\circ }\) bend tube [21], where RSM performed best in comparison to other commonly used eddyviscosity models. Especially for turbulent kinetic energy, while the eddyviscosity models tend to underpredict the DNSbased values, RSM can capture the behavior better. We have shown this also in our results, where RSMbased turbulent kinetic energy showed slightly higher peaks in the curved part of the phantom. Based on this and direct comparison of performances of turbulence models with MRI measurements in the phantom, we conclude that the RSM turbulence model is the most suitable to properly capture the most important flow features. Additionally, in terms of computational efficiency, although a larger number of transport equations needs to be solved by the RSM turbulence model when compared to the eddyviscositybased models, its computational costs are still much smaller when compared to highfidelity LES or DNS methods (i.e., \({{{\mathcal {O}}}} (10^2  10^3)\) faster, respectively), which makes it a good choice for the patientspecific clinical applications.
Wall shear stress based on CFD and MRI
Two main points need to be addressed when comparing simulations (CFD) and experiments (MRI): (i) the absolute values of the WSS, and (ii) the local distributions of the WSS, respectively. Generally, we observed consistent lower values of WSS from MRI in comparison to CFD results. This can be explained in terms of the lower spatial resolution of MRI—especially in the proximity of the wall. The absolute values of the spatially averaged \(\hbox {WSS}_{\mathrm{mean}}\) from CFD simulations for the phantom are \(3.7\times \) (CFD), \(2.5\times \) (\(\hbox {DCFD}_{0.2\times 0.2\times 0.2 \,{\mathrm{mm}}^3}\)), and \(1.4\times \) (\(\hbox {DCFD}_{0.7\times 0.7\times 1.5 \,{\mathrm{mm}}^3}\)) higher than the MRI values, respectively. Similarly, the peak WSS in the stenotic region (\(\hbox {WSS}_{\mathrm{st}}\)) for CFD simulations are \(8.3\times \) (CFD), \(3.9\times \) (\(\hbox {DCFD}_{0.2\times 0.2\times 0.2 \,{\mathrm{mm}}^3}\)), and \(1.83\times \) (\(\hbox {DCFD}_{0.7\times 0.7\times 1.5 \,{\mathrm{mm}}^3}\)) higher than the MRI values. Note that the stenotic region is the most sensitive one due to a sudden flow acceleration, and a reduction of the number of voxels (since the spatial resolution of MRI is fixed). In contrast, the CFD wall resolution in this region is increased since the identical number of control volumes as for the healthy segment is now distributed over the reduced area of stenotic crosssection. It can be seen that reduction of spatial resolution of CFD, lowers values of WSS. A similar systemic undersolving of WSS by MRI was also reported in [7].
In contrast to the absolute values of WSS, the local distributions of WSS calculated from CFD and measured by MRI exhibit more similarities. To illustrate this, we scaled the local WSS with the spatially averaged mean WSS (\(\hbox {WSS}_{\mathrm{mean}}\)) of each modality, as shown in Fig. 4c. This approach enables us to compare variations of WSS associated with locally elevated or suppressed distributions for the reference averaged value. The circumferentially averaged mean WSS profiles, shown in Fig. 5, show a good agreement except at the stenotic region. It can be seen that the reduction of spatial resolution of CFD reduced the peak values, but also introduced a shift of the peak location. Similar behavior was shown also in a related study of [22], but these findings were not addressed. This shift should be taken into account when analyzing cases where the exact location of the peak WSS is of importance.
Effect of assumptions in CFD
The reliability and accuracy of CFD simulations in studying blood flow in the patientspecific conditions are directly connected to realistic representations of vessel geometry, as well as the imposed boundary conditions. The geometry representation can affect the simulations by not including all of the side branches [23] and due to the segmentation variability [24]. For the aortic coarctation studies, the choice of proper inlet and outlet sidebranching boundary conditions was highlighted in [19] and [25], respectively. Finally, the assumption of rigidwall in patientspecific simulations can lead to differences up to 30% [26].
Our approach to perform analysis for a simplified phantom is based on a stepbystep elimination process of specific contributions (e.g., exact wall geometry with all side branches, the exact specification of the inlet and outlet boundary conditions, wall elasticity) which make a fair comparison between CFD and MRI for the patientspecific cases difficult. By eliminating the effects of the abovementioned limitations, we have achieved identical phantom working conditions for simulations and experiments. Despite achieving a good agreement between CFD and MRI for the mean velocity profiles at various crosssections of the phantom, the local distribution of the WSS still exhibited significantly different values. A similar level of disagreement was also reported in the literature [17, 18]. Based on the abovepresented arguments, we conclude that the major contributor to disagreement between MRI and CFD is due to limitations in spatial resolution of MRI in the proximity of the aortic wall.
The effect of preprocessing on the outcomes of simulations is also highlighted in the patientspecific CoA. For this case, a more drastic smoothing of the aortic wall resulted in a decrease in mean WSS of 7.3%, bringing the value closer to MRI. However, this decrease is still relatively small, comparing to almost an orderofmagnitude difference in WSS between MRI and CFD and both simulations resulted in very similar global distributions of WSS.
Our thorough comparison of blood flow and WSS based on MRI and CFD showed, that the blood flow, in general, agreed well for both techniques. However, the derived variables, like WSS, deviate much more. The reasons behind the deviation can be found on both sides. As demonstrated in the phantom, the underestimation of WSS is mostly due to MRI not being able to capture the steep gradients in a region with sudden flow acceleration. However, the boundary condition treatment and preparation of the geometry in CFD has also a big effect. Patientspecific CoA showed that the application of flat parabolic profile leads to discrepancies in the ascending aorta, and using arterial geometry without adequate smoothing slightly overestimates CFDbased WSS. Hence, when evaluating the WSS differences between MRI and CFD attention should be given to the limitation of both techniques.
Clinical applications
While in this work we present the MRICFD coupling using 4Dflow MRI data, the technique can be coupled with different imaging techniques for the acquisition of the anatomy—e.g., MRA imaging. Nevertheless, PCMRI at the inlet is still necessary for accurate definition of boundary conditions [19]. With this integrated CFDMRI approach, where CFD is based on the advanced RSM turbulence model, it is possible to generate simulation results with high spatial resolution and a high level of accuracy within a few hours. This may lead to several interesting clinical applications of the imagebased CFD framework for diseased arteries. Potential examples include pre and postoperative followup for patients suffering from CoA [12]. Furthermore, by studying a wider population of patients, biomarkers for restenosis could be identified similarly as was done for femoropopliteal arteries [27], which could lead to a predictive method for potential complications connected to CoA.
The proposed approach is not only limited to coarctation or aorta. It can be easily applied also in other aortic diseases (e.g., aneurysm and dissection), or to study the blood flow in different parts of the cardiovascular system. Examples of these applications have already been tested, for example in cerebral aneurysms [28], stenosis of major arteries [29], and pulmonary arteries [30], and show very good promise.
Finally, it is important to touch upon the feasibility of using imagebased CFD in the clinical application. As we showed with this study, the methods have a great potential to study the longterm effects of diseases or to model the progression and predict the outcomes of chronic vascular diseases. However, the state of the methods at the moment does not allow for implementation in acute decisionmaking. This is especially due to the fact, that the simulations and their preparation is timeconsuming and requires expert knowledge. Most medical doctors do not have adequate training to perform such simulations and therefore, experts in computational fluid dynamics should be involved. Hence, until improvements in the automatizing the imagebased CFD are made, for example by implementing meshless methods (e.g., Solid Particle Hydrodynamics [31]), the usage of MRI or other imaging methods for diagnostics of acute cases is necessary.
Limitations
For the MRI measurements in the simplified phantom, the 4Dflow MRI with nonsegmented gradientecho and echoplanar imaging (EPI) acceleration has been used. In contrast to that, for in vivo patientspecific CoA, the 4Dflow MRI with segmented gradientecho without EPI has been applied. As result, some minor differences in accuracy in velocity quantitation may be present between these two experiments. Additionally, water was used as a working fluid for the flow in the phantom instead of blood mimicking fluid with nonNewtonian viscosity, e.g., as proposed by Cheng et. al. [32]. Use of such a fluid would represent the blood flow more adequately [32], however, since the shear rate in the aorta was relatively high the nonNewtonian effects should be minimal. The CFD simulations of the patientspecific CoA have been performed with the rigid walls assumption. The dynamic movement of the aorta is present in vivo 4Dflow MRI, which can produce differences in local distributions of WSS along the arterial wall. To circumvent the effects of the dynamic movement of the aorta during the cardiac cycle, the simplified phantom geometry has been considered too. Despite the relatively high Reynolds number of flow in the patientspecific aorta, some local nonNewtonian effects can take place, which is currently not taken into account in CFD simulations. We have considered a single case of the patientspecific CoA and have performed a comparative assessment of 4Dflow MRI and CFDRSM as a first proofofconcept. This study can be easily extended with a larger number of patientspecific aorta conditions.
Conclusions
With this study, we showed that MRIbased CFD simulations are a good alternative tool to use in studying the blood flow in CoA. Using MRI and MRIbased simulations, we assessed the blood flow in a phantom representing a simplified CoA and in a patientspecific aorta with coarctation.
Due to the narrowing and relatively high Re, the flow in the phantom was of turbulent nature. Because of this, a choice of turbulent model had to be made. We have compared \(k\epsilon \), SST, and RSM. The differences between the turbulence model arise after the bend—where the flow gets more complex. The lowerorder CFD models (\(k\epsilon \), SST) cannot accurately model the secondary flow motion that naturally appears in these types of geometries as is shown in the results obtained from MRI. However, the RSM model can predict these motions as was shown for both velocity magnitude as well as the outofplane vorticity. Thanks to using the phantom, where the boundary conditions between MRI and CFD are identical, we were able to accurately study WSS, an important parameter that is often regarded as a biomarker for CoA. We showed that WSS based on MRI is approximately fourtimes lower than WSS based on CFD, however, agrees well in terms of the local distribution.
Finally, we have applied MRICFD coupling on patientspecific CoA to demonstrate usability of the technique for clinical applications. The agreement between 4Dflow MRI and CFD was good in terms of velocity. For WSS, simulations showed again higher values that MRI, however, the local regions of high/low WSS agree well between the different techniques. However, the simulations bring several advantages, like higher resolution and better prediction of the absolute values of WSS. This shows, that imagebased simulations are a good technique to asses the state of this pathology.
Methods
Studied cases
Phantom
The phantom for this study was a 3Dprinted computeraided design object, mimicking the simplified aorta with coarctation. It represents a \(180^{\circ }\) bend structure in which obstruction is present in the distal leg, as shown in Fig. 7a. The threedimensional Utube phantom is fabricated from Duraform Flex material with characteristic stiffness of a shore hardness scale of A75, manufactured by Materialise. The phantom can be considered as a rigid structure. To allow magnetic resonance imaging, the 3D Utube phantom is placed inside a 10 l jerry can, which is filled with gelatin. Gelatin is used instead of water to prevent movement of the liquid due to the sound produced by the MRI system. The exact composition of the gelatin is 9.9 l water, 600 g gelatin, 100 ml paraben, and 1.5 ml Gadovist. The jerry can has two connectors to allow connection with the tubes through which the liquid will be pumped. A constant throughput is provided using a pump connected to the inlet tube with a prescribed flow of \(\dot{Q_0}=4.5\) l/min. Additional morphometric and flow characteristics of the studied phantom can be found in Table 3.
Patientspecific aorta with coarctation
This study protocol was approved by the Medical Ethics Committee of the Leiden University Medical Center (P14.095), and informed consent was signed by both parents/guardians of the subject. The patient with CoA included in this study was female, 14.5 years old, 164 cm, 51,9 kg, and had a tricuspid aortic valve. The final geometry includes both thoracic and abdominal aorta (AbAo) with six branching arteries: brachiocephalic trunk (BT), left common carotid artery (LCCA), left subclavian artery (LSA), celiac trunk (CT), and left/right renal artery (L/RRA) (visualized in Fig. 7c. Additional morphometric and flow characteristics of the studied CoA can be found in Table 3.
Magnetic resonance imaging
For both, the phantom and patientspecific CoA, 4Dflow MRI was performed on a 3T MRI system (Ingenia, Philips Healthcare, Best, The Netherlands). For the patientspecific CoA, aortic 4Dflow MRI was performed using a hemidiaphragm respiratory navigator with retrospective electrocardiogram gating without echoplanar imaging. Additional file 1 about the MRI sequences can be found in Table 4.
The acquired 4Dflow MRI data sets were afterward analyzed using CAAS MR Solutions v5.0 (Pie Medical Imaging BV, Maastricht, The Netherlands). The analysis is similar for both phantom and patientspecific CoA and the most important differences are highlighted. In both, the analysis was initialized by manual placement of starting and ending points. For the phantom, the two points were placed at the same level in the opposite legs, whereas for the patientspecific CoA, the starting point was placed in the aortic root, and the ending points were placed in the abdominal aorta and the six branching arteries. A phasespecific 3D volume was automatically segmented for the specific phase and copied to all phases. The 3D segmentation uses a deformable model algorithm that recursively optimizes the location of the surface towards the vessel luminal boundary based on image gradients, extracted from the appropriate phase within the 4Dflow MRI data, while simultaneously maintaining local smoothness of the 3D segmented surface, [34]. Manual delineation of the vessel lumen boundary was applied with the available adaptation tool from the software in case of segmentation incorrectness for the patientspecific CoA.
CFD model
Numerical mesh
To perform CFD simulations of the patientspecific aorta geometry (obtained from 4Dflow MRI segmentation using CAAS MR Solutions v5.0), we first used Vascular Modeling Toolkit for the geometry reconstruction [35]. To check the sensitivity of CFD results on the imposed levels of geometry smoothing, we have used the rough (\(\hbox {CFD}_{\mathrm{rough}}\)) (i.e., by using a Taubin filter with 100 iterations and passband settings of 0.4), and the smoothed geometry (\(\hbox {CFD}_{\mathrm{smooth}}\)) (i.e., by using a Taubin filter with 100 iterations and passband settings of 0.01). For both geometries, cylindrical extensions were added at all outlets with a length of \(1.5d_0\) (where \(d_0\) is the diameter of the individual blood vessel where the outlet is located). For both phantom and patientspecific aorta simulations, we have employed a hybrid numerical mesh containing prismatic elements in the proximity of the wall (to properly resolve characteristic boundary layers), while the tetrahedrons were used in the central part of the domain. Details about the mesh sizing are shown in Table 3.
We have performed a meshindependency study and the final numerical mesh for the phantom case was approximately 14 million control volumes, while approximately 7 million control volumes were used for the aorta case. Note that the larger number of control volumes for the phantom geometry was due to the addition of a segment with a length of \(10 d_0\) to have a proper capture of the poststenotic flow region (not shown in Fig. 7a).
Governing equations
Since we are dealing with a fully (the phantom case) or a partially (the aorta case) developed turbulent flow regimes, we adopt an unsteady RANS approach to model turbulence. We apply three different classes of turbulence models: two based on the eddyviscosity concept (standard \(k\epsilon \) and shearstress transport (SST) \(k\omega \) model), and an advanced turbulence model based on solving a complete set of the turbulent stress components (the full RSM). Despite its superior theoretical foundation when compared to the eddyviscosity turbulence models [20], applications of the RSM model are very scarce in biomedical flow applications. Here we propose the use of the RSM as an alternative to a highfidelity DNS or LES approaches. The following set of the governing equations is introduced for the abovementioned turbulence models:

the standard eddyviscosity \(k\epsilon \) model with enhanced wall treatment [36]: PDEs for (\(U_i  p  k  \varepsilon \))

the lowReynolds shear stress transport (SST) \(k\omega \) model [37]: PDEs for (\(U_i  p  k  \omega \))

the Reynolds stress model (RSM) with the linear pressure strain term and enhanced wall treatment [38]: PDEs for (\(U_i  p  \overline{u_i u_j}  \varepsilon \))
where \(U_i\) is the velocity vector, p is the pressure, k is turbulent kinetic energy, \(\varepsilon \) is dissipation rate, \(\omega \) is turbulent frequency, and \(\overline{u_i u_j}\) is turbulent stress tensor.
Boundary and initial conditions
For the phantom, the imposed volumetric flow rate of 4.5 l/min is identical to the experimental conditions and corresponds to the inlet Reynolds number of \(Re=4539\) (\(Re=V_0\cdot D_0/\nu \)). For the patientspecific aorta, the timedependent inlet conditions are matched with MRI measurements during the entire cardiac cycle. We have extracted the measured volumetric flow rate at the inlet plane (\(Q_0\)), and converted it to the characteristic mass flow rate (\({\dot{m}}=Q_0 \cdot \rho _{\mathrm{blood}}\)). The mass flow rates were fitted with a smooth spline with piecewise polynomial (with a smoothing parameter \(p=0.99999947\) and \(R^2=0.9995\)), Fig. 7b, which gives the following range of the inlet Reynolds number, \(0 \le Re \le 6276\), and corresponding Womersley number of \(Wo=D_0 \left( \omega _f/\nu \right) ^{1/2} = 29\). In total, we have simulated five cardiac cycles, to obtain results without the influence of initial conditions. Only the last cycle was used for the analysis. For all turbulence parameters, the uniform inlet values were imposed with the following specifications: the intensity of turbulence of 5%, the ratio of turbulent and molecular viscosity (\(\mu _t/\mu \)) of 10, the isotropic assumption of normal turbulent stress components (\(\overline{u_i u_i}=2/3 k\)), and zero values of the turbulent shear stress components (\(\overline{u_i u_j}=0\)). At outlets, a zero diffusion flux was imposed for all transport variables. For the patientspecific aorta, a predefined fixed (MRIbased) percentage of the inlet flow rate was prescribed. The noslip velocity boundary condition was imposed at the walls of blood vessels, and the model was assumed to be rigid.
Physical properties and simulation setup
For the phantom, water was used as a working fluid (\(\rho =998\, \hbox {kg}/\hbox {m}^3\), \(\mu =1.003\,\hbox {mPa}\cdot \hbox {s}\)). For the aorta, the real blood properties were assumed for the simulations (\(\rho =1060\,\hbox {kg}/\hbox {m}^3\), \(\mu =3.5\,\hbox {mPa}\cdot \hbox {s}\)). It was previously demonstrated that the assumption of constant blood viscosity is adequate for aortic blood simulations [39]. The simulations were performed using Ansys Fluent 19.1 (Ansys Inc., Canonsburg, Pennsylvania, USA) with the following simulation settings:

Solver—pressure based

Pressure–velocity coupling—SIMPLE

Spatial discretization

Gradientleastsquares cellbased

Pressure—second order

Momentum—secondorder upwind

Turbulence variables—secondorder upwind


Temporal discretization (CoA case)

Fully implicit secondorder scheme

Time step \(\Delta t = 0.0005\)s


Residuals (all) \(10^{5}\)
Analysis
Downsizing and mapping
For additional analysis and voxeltovoxel comparison, the phantom CFD velocity data were downsampled by applying a bilinear interpolation on two different equidistant meshes (DCFD):

\(\hbox {DCFD}_{0.7\times 0.7 \times 1.5}\): voxel resolution \(0.7\times 0.7\times 1.5\,\hbox {mm}^3\) (identical to MRI)

\(\hbox {DCFD}_{0.2\times 0.2 \times 0.2}\): voxel resolution \(0.2\times 0.2\times 0.2\,\hbox {mm}^3\).
Both downsized CFD velocity fields were analyzed using CAAS MR Solutions v5.0 (Pie Medical Imaging BV, Maastricht, The Netherlands) to obtain WSS. Additionally, the WSS distribution along the vessel walls was mapped on a 2D surface where the horizontal axis indicates the nondimensional radial distance from the vessel centerline (\(\pi \le r \le +\pi \)), and the vertical axis indicates the nondimensional arclength of the vessel (\(0 \le l/l_0 \le 1\)). This mapping was done by an originally developed inhouse tool in Matlab R2019a (MathWorks, Inc., Natick, Massachusetts, U.S.A.). This approach has provided an easy and objective comparison between different results.
Vorticity calculation
Vorticity (\({\varvec{\omega }}\)) was calculated from MRIbased velocity components using Matlab R2019a (MathWorks, Inc., Natick, Massachusetts, U.S.A.) as
where \({\varvec{u}}\) is the velocity vector. In the numerical procedure, the partial derivatives were calculated using a central differencing scheme for the interior points and a singlesided (forward) difference scheme for the edges.
Wall shear stress calculation
The WSS for MRI and the downsized CFD data sets were calculated based on the extracted velocity profile perpendicular to the phasespecific segmented 3D surface using CAAS MR Solutions v5.0. After factorizing the velocity profile into its component parallel to the lumen wall, WSS was computed by the first derivative of a quadratic approximation of that velocity profile at the location of the lumen wall as
where \(U_{}\) is the wallparallel velocity component and n is the wallnormal direction. For the CFD simulations involving turbulence models, the wall shear stress is directly available from calculations of the wallparallel velocity component based on the enhanced wall treatment.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
References
Yamaguchi T, Kikkawa S, Tanishita K, Sugawara M. Spectrum analysis of turbulence in the canine ascending aorta measured with a hotfilm anemometer. J Biomech. 1988;21:489–95. https://doi.org/10.1016/00219290(88)902412.
Krasuski RA, FouadTarazi F. Chapter 74  coarctation of the aorta. In: Lip GY, Hall JE, editors. Comprehensive Hypertension. Philadelphia: Mosby; 2007. p. 923–30. https://doi.org/10.1016/B9780323039611.500775.
Abe H, Kawamura H, Choi H. Very largescale structures and their effects on the wall shearstress fluctuations in a turbulent channel flow up to re\(\tau \)=640. J Fluids Eng. 2004;126:835–43. https://doi.org/10.1115/1.1789528.
Feng Y, Wada S, Tsubota K, Yamaguchi T. A modelbased numerical analysis in the early development of intracranial aneurysms. In: 2005 IEEE engineering in medicine and biology 27th annual conference; 2005. p. 607–10.
Dyverfeldt P, Hope MD, Tseng EE, Saloner D. Magnetic resonance measurement of turbulent kinetic energy for the estimation of irreversible pressure loss in aortic stenosis. JACC Cardiovasc Imaging. 2013;6:64–71. https://doi.org/10.1016/j.jcmg.2012.07.017.
Juffermans JF, Nederend I, van den Boogaard PJ, ten Harkel ADJ, Hazekamp MG, Lamb HJ, Roest AAW, Westenberg JJM. The effects of age at correction of aortic coarctation and recurrent obstruction on adolescent patients: MRI evaluation of wall shear stress and pulse wave velocity. Eur J Radiol. 2019;3:24. https://doi.org/10.1186/s4174701901029.
Zimmermann J, Demedts D, Mirzaee H, Ewert P, Stern H, Meierhofer C, Menze B, Hennemuth A. Wall shear stress estimation in the aorta: impact of wall motion, spatiotemporal resolution, and phase noise. J Magn Reson Imaging. 2018;48:718–28. https://doi.org/10.1002/jmri.26007.
Rinaudo A, D’Ancona G, Baglini R, Amaducci A, Follis F, Pilato M, Pasta S. Computational fluid dynamics simulation to evaluate aortic coarctation gradient with contrastenhanced CT. Comput Methods Biomech Biomed Engin. 2015;18:1066–71. https://doi.org/10.1080/10255842.2013.869321.
Zhu Y, Chen R, Juan YH, Li H, Wang J, Yu Z, Liu H. Clinical validation and assessment of aortic hemodynamics using computational fluid dynamics simulations from computed tomography angiography. BioMed Eng Online. 2018;17:53. https://doi.org/10.1186/s1293801804855.
Gaze DC. Congenital heart disease. In: IntechOpen. 1st ed. 2018. https://doi.org/10.5772/intechopen.74138.
Yang F, Zhai B, Hou LG, Zhang Q, Wang J. Computational fluid dynamics in the numerical simulation analysis of endtoside anastomosis for coarctation of the aorta. J Thorac Dis. 2018;10(12):6578–84.
Andersson M, Lantz J, Ebbers T, Karlsson M. Quantitative assessment of turbulence and flow eccentricity in an aortic coarctation: impact of virtual interventions. Cardiovasc Eng Technol. 2015;6:281–93. https://doi.org/10.1007/s132390150218x.
Arzani A, Dyverfeldt P, Ebbers T, Shadden SC. In vivo validation of numerical prediction for turbulence intensity in an aortic coarctation. Ann Biomed Eng. 2011;40:860–70. https://doi.org/10.1007/s1043901104476.
Gårdhagen R, Lantz J, Carlsson F, Karlsson M. Large eddy simulation of stenotic flow forwall shear stress estimation  validation and application. WSEAS Trans Biol Biomed. 2011;8:86–101.
Wols B. Computational fluid dynamics in drinking water treatment. KWR watercycle research institute series. London: IWA Publishing; 2011.
Malhotra A, Mousa KH. Turbulence modelling in pipe flow. Math Comput Model. 1990;14:755–60. https://doi.org/10.1016/08957177(90)90283S.
Szajer J, HoShon K. A comparison of 4d flow MRIderived wall shear stress with computational fluid dynamics methods for intracranial aneurysms and carotid bifurcations—a review. Magn Reson Imaging. 2018;48:62–9.
Miyazaki S, Itatani K, Furusawa T, Nishino T, Sugiyama M, Takehara Y, Yasukochi S. Validation of numerical simulation methods in aortic arch using 4d flow MRI. Heart Vessels. 2017;32:1032–44. https://doi.org/10.1007/s0038001709792.
Goubergrits L, Mevert R, Yevtushenko P, Schaller J, Kertzscher U, Meier S, Schubert S, Riesenkampff E, Kuehne T. The impact of MRIbased inflow for the hemodynamic evaluation of aortic coarctation. Ann Biomed Eng. 2013;41:2575–87. https://doi.org/10.1007/s1043901308792.
Hanjalić K, Launder B. Modelling turbulence in engineering and the environment: secondmoment routes to closure. Cambridge: Cambridge University Press; 2011. https://doi.org/10.1017/CBO9781139013314.
Yang X, Tucker PG. Assessment of turbulence model performance: Large streamline curvature and integral length scales. Comput Fluids. 2016;126:91–101. https://doi.org/10.1016/j.compfluid.2015.11.010.
Puiseux T, Sewonu A, Meyrignac O, Rousseau H, Nicoud F, Mendez S, Moreno R. Reconciling pcMRI and CFD: an invitro study. NMR Biomed. 2019;32:e4063. https://doi.org/10.1002/nbm.4063.
Yevtushenko P, Hellmeier F, Bruening J, Kuehne T, Goubergrits L. Numerical investigation of the impact of branching vessel boundary conditions on aortic hemodynamics. Curr Dir Biomed Eng. 2017;3:321–4. https://doi.org/10.1515/cdbme20170066.
Perinajová R, Juffermans JF, Westenberg JJ, van der Palen RL, van den Boogaard PJ, Lamb HJ, Kenjereš S. Geometrically induced wall shear stress variability in CFDMRI coupled simulations of blood flow in the thoracic aortas. Comput Biol Med. 2021;133:104385. https://doi.org/10.1016/j.compbiomed.2021.104385.
Madhavan S, Kemmerling EMC. The effect of inlet and outlet boundary conditions in imagebased CFD modeling of aortic flow. BioMed Eng Online. 2018;17:66. https://doi.org/10.1186/s1293801804971.
Lantz J, Renner J, Karlsson M. Wall shear stress in a subject specific human aorta – influence of fluidstructure interaction. Int J Appl Mech. 2011;03:759–78. https://doi.org/10.1142/S1758825111001226.
Gökgöl C, Diehm N, Räber L, Büchler P. Prediction of restenosis based on hemodynamical markers in revascularized femoropopliteal arteries during leg flexion. Biomech Model Mechanobiol. 2019;18:1883–93. https://doi.org/10.1007/s10237019011839.
Castro M, Putman C, Cebral J. Computational fluid dynamics modeling of intracranial aneurysms: effects of parent artery segmentation on intraaneurysmal hemodynamics. Am J Neuroradiol. 2006;27:1703–9.
Zhang D, Xu P, Qiao H, Liu X, Luo L, Huang W, Zhang H, Shi C. Carotid DSA based CFD simulation in assessing the patient with asymptomatic carotid stenosis: a preliminary study. BioMed Eng Online. 2018;17:31. https://doi.org/10.1186/s1293801804659.
Rijnberg FM, van der Woude SFS, van Assen HC, Juffermans JF, Hazekamp MG, Jongbloed MRM, Kenjeres S, Lamb HJ, Westenberg JJM, Wentzel JJ, Roest AAW. Nonuniform mixing of hepatic venous flow and inferior vena cava flow in the fontan conduit. J R Soc Interface. 2021;18:20201027. https://doi.org/10.1098/rsif.2020.1027.
Caballero A, Mao W, Liang L, Oshinski J, Primiano C, McKay R, Kodali S, Sun W. Modeling left ventricular blood flow using smoothed particle hydrodynamics. Cardiovasc Eng Technol. 2017;8:465–79. https://doi.org/10.1007/s132390170324z.
Cheng AL, Wee CP, Pahlevan NM, Wood JC. A 4d flow MRI evaluation of the impact of sheardependent fluid viscosity on in vitro fontan circulation flow. Am J Physiol Heart Circ Physiol. 2019;317:H1243–53. https://doi.org/10.1152/ajpheart.00296.2019.
Fung YC. Biomechanics: circulation. 2nd ed. Berlin: Springer; 1997.
Delingette H. General object reconstruction based on simplex meshes. Int J Comput Vis. 1999;32:111–46. https://doi.org/10.1023/A:1008157432188.
Antiga L, Piccinelli M, Botti L, EneIordache B, Remuzzi A, Steinman DA. An imagebased modeling framework for patientspecific computational hemodynamics. Med Biol Eng Comput. 2008;46:1097. https://doi.org/10.1007/s1151700804201.
Launder BE, Spalding DB. Lectures in mathematical models of turbulence [by] B. E. Launder and D. B. Spalding. New York: Academic Press London; 1972.
Menter FR. Twoequation eddyviscosity turbulence models for engineering applications. AIAA J. 1994;32:1598–605. https://doi.org/10.2514/3.12149.
Gibson MM, Launder BE. Ground effects on pressure fluctuations in the atmospheric boundary layer. J Fluid Mech. 1978;86:491–511. https://doi.org/10.1017/S0022112078001251.
Soares AA, Gonzaga S, Oliveira C, Simões A, Rouboa AI. Computational fluid dynamics in abdominal aorta bifurcation: nonnewtonian versus newtonian blood flow in a real case study. Comput Methods Biomech Biomed Engin. 2017;20:822–31. https://doi.org/10.1080/10255842.2017.1302433.
Acknowledgements
Not applicable.
Funding
This research was financially supported by a grant from the Dutch Heart Foundation, The Netherlands (Grant Number CVON201808RADAR).
Author information
Authors and Affiliations
Contributions
RP prepared the models, ran the CFD simulations, analyzed and interpreted the data, and was a major contributor in writing the manuscript. JJ acquired the MRI data and postprocessed the patientdata. JLM prepared the models, ran the CFD simulations and analyzed the phantom data. JPA developed the phantom model, acquired the MRI data of phantom, analyzed and interpreted the data and was a contributor in writing the manuscript. LL developed the phantom model, analyzed, and interpreted the data. JW acquired the MRI data, analyzed and interpreted the data and was a contributor in writing the manuscript. HL conceived the presented study, and interpreted the data. SK conceived the presented study, analyzed and interpreted the data and was a major contributor in writing of manuscript and final editing. All authors read and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
This study protocol was approved by the Medical Ethics Committee of the Leiden University Medical Center (P14.095).
Consent for publication
Informed consent was signed by both parents/guardians of the subject.
Competing interests
The authors declare that they have no competing interests.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Supplementary Information
Additional file 1.
Mesh Dependency Study. Table S1. Mesh dependency analysis for phantom and patient, with wall shear stress for three different meshes (fine  1, medium  2, coarse  3), refinement ratio r, Richardson extrapolation (fh=0), Grid Convergence Index (GCI1,2  finemedium, GCI2,3  mediumcoarse), and the test whether the studied variables lie in the asymptotic range. Figure S2. Tetrahedral mesh with details of inlet and stenosed region for phantom (a) and patientspecific aorta (b).
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Perinajová, R., Juffermans, J.F., Mercado, J.L. et al. Assessment of turbulent blood flow and wall shear stress in aortic coarctation using imagebased simulations. BioMed Eng OnLine 20, 84 (2021). https://doi.org/10.1186/s12938021009214
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s12938021009214
Keywords
 Magnetic resonance imaging
 Computational fluid dynamics
 Turbulence
 Aorta
 Coarctation
 Phantom