Differential impact of local stiffening and narrowing on hemodynamics in repaired aortic coarctation: an FSI study

Abstract

Even after successful treatment of aortic coarctation, a high risk of cardiovascular morbidity and mortality remains. Uncertainty exists on the factors contributing to this increased risk among which are the presence of (1) a residual narrowing leading to an additional resistance and (2) a less distensible zone disturbing the buffer function of the aorta. As the many interfering factors and adaptive physiological mechanisms present in vivo prohibit the study of the isolated impact of these individual factors, a numerical fluid–structure interaction model is developed to predict central hemodynamics in coarctation treatment. The overall impact of a stiffening on the hemodynamics is limited, with a small increase in systolic pressure (up to 8 mmHg) proximal to the stiffening which is amplified with increasing stiffening and length. A residual narrowing, on the other hand, affects the hemodynamics significantly. For a short segment (10 mm), the combination of a stiffening and narrowing (coarctation index 0.5) causes an increase in systolic pressure of 58 mmHg, with 31 mmHg due to narrowing and an additional 27 mmHg due to stiffening. For a longer segment (25 mm), an increase in systolic pressure of 50 mmHg is found, of which only 9 mmHg is due to stiffening.

Appendix: Mesh and time step sensitivity study

Given the focus of the study on pressure, the criterion for the grid convergence and time step dependency study was the accuracy of the predicted pressures along the aortic arch. The case used for the analysis was the case with the shortest stenosis length (1 cm) and highest degree of stenosis (coarctation index 0.5).

1.1 Mesh sensitivity analysis

Four different, full hexahedral uniform meshes (R1, R2, R3 and R4) with an increasing number of elements in the boundary layer (ranging from 4 to 8 layers), the transversal and axial direction were constructed, with R4 considered as reference. A conforming mesh was applied in the fluid and solid domain. The number of cells is depicted in Table 1, together with the calculation time required to compute one cardiac cycle (on two 10-core Intel Xeon E5-2680v2 processors).

Table 1 Grid refinement study of the pressure in an FSI model of aortic coarctation

As the flow distal to the stenosis is complex and highly disturbed (Reynolds numbers up to 11,139), a high mesh density is required in this region to resolve the flow field in space. This is realized by locally adapting the fluid mesh. The resulting fluid mesh (R5) has, compared to the finest mesh (R4), a higher mesh density in the coarctation zone, but a coarser grid proximal to the stenosis and in the lower part of the descending aorta (see Fig. 8).

Fig. 8

Left indication of sections where pressure was calculated (p proximal; d distal). Right a mesh for the fluid domain (blue) and the arterial wall (red) of an aortic arch with aortic coarctation. Note the axial coarsening toward the descending aorta (R5). b, c The cross-sectional grids of the fluid mesh at the coarctation (coa) and the descending aorta (desc), which result from multiblock structures R4 (uniform grid refinement) and R5 (local grid refinement) (color figure online)

Figure 9 depicts the pressure along the centerline of the aorta at peak systole, with Table 1 tabulating the mean error of the pressure in different cross sections (indicated by the dashed lines in Fig. 9). These errors are defined with respect to the reference grid R4 and relative to the pressure amplitude in the corresponding cross section. The mean error thereby denotes the error averaged over one cardiac cycle and over the respective cross section. From the results in Table 1, it can be seen that even for meshes with a low cell density, the mean errors proximal, halfway and distal to the coarctation zone remain low (<2 %). When comparing the locally refined grid R5 with the uniformly refined grids R3 and R4, an important reduction in computation time is gained (23 h 38 min per cardiac cycle vs 30 h 33 min and 40 h 38 min) without a loss in accuracy. The mean error obtained with the mesh R5 stays below 1.15 %, and comparable errors are found as for the mesh R3.

Fig. 9

Pressure along the centerline at peak systole for increasing mesh densities (R1–R4) and a grid with a local refinement at the coarctation region and a gradual coarsening toward the descending aorta (R5)

1.2 Time step sensitivity analysis

Figure 10 illustrates the impact of the time step size on the pressure evolution at different cross sections along the aorta. It can be observed that the results in the proximal part and at the coarctation zone are more or less time step independent, whereas the small pressure oscillations in the distal part are not captured with a large time step size (of 4 or 5 ms). Moreover, the oscillations developing in d3 are not even resolved properly with a time step size of 1 ms. As such, the time step size was further decreased to 0.5 and 0.25 ms and the results are shown in Fig. 11. An unstable behavior was found if a small time step size was applied. The observed oscillations responsible for this behavior were indeed not resolved for the simulations using larger time step sizes. Because the oscillations itself are resolved by multiple time steps and the frequency of the oscillations is more or less time step independent, it is presumed that these oscillations do not arise from a numerical instability but have a physical origin, triggered by the disturbed blood flow. In a physiological setting, this oscillation would, however, be cushioned by the damping nature of the surrounding tissue. We believe that the lack of physical damping in our model resulted in the observed oscillations that eventually got unstable when using a time step size smaller than 1 ms.

Fig. 10

Influence of the time step size on the pressure evolution at proximal cross section p1, halfway the coarctation zone (coa) and at two distal cross sections (d1 and d3). See Fig. 8 for an indication of these plane locations

Fig. 11

Left detail of the pressure evolution, illustrating the temporal resolution of the oscillations. Right influence of Rayleigh damping on the pressure evolution at cross section d2. Inclusion of Rayleigh damping prevents the simulation from unstable behavior

To test this hypothesis, Rayleigh damping was added to the structural model and the simulation using a time step size of 0.5 ms was repeated. The Rayleigh damping coefficients α and β were selected such that 1 % damping of the waves with a 1-Hz frequency (close to the frequency of the cardiac cycle) was obtained and 20 % damping for the 250-Hz waves (i.e., the frequency of the observed oscillations). These constrictions resulted in a value of 0.116 for the mass proportional damping parameter α and 0.000255 for the stiffness proportional damping β. It is demonstrated in Fig. 11 that the unstable behavior indeed disappears with the use of Rayleigh damping.

The larger pressure oscillations at the start of the simulation (t < 0.2 s) for the case with Rayleigh damping is explained by the temporal discretization schemes used at the start. The simulation without damping is started with a first-order scheme, to facilitate the startup. After 0.2 s, the accuracy is improved by switching to a second-order scheme. For the case with damping, a second-order scheme can be used from the start on. In this study, a time step size of 2 ms has been used as a compromise between accuracy and computation time. The error obtained with this time step size is sufficiently smaller than the mutual differences in results.