3.1. Phase-averaged flow analysis

Phase-averaged axial velocity fields are reported in figure 5 at different cross-sections ($z/D=1.5$, $z/D=2.5$, $z/D=3.5$) of successive cycle stages, as they provide meaningful insights about the different flow features of the three cases investigated. In the PF stage the typical triple jet pattern is found to dominate the mean flow field, in agreement with previous investigations on SJMVs, both with straight (Guivier-Curien, Deplano & Bertrand Reference Guivier-Curien, Deplano and Bertrand2009; Hatoum, Maureira & Dasi Reference Hatoum, Maureira and Dasi2020) and anatomic (Ha et al. Reference Ha, Kim, Kweon, Lee, Kim, Kim and Yang2016) aortic root configurations. The streamwise mean velocity profile at PF, over the valve symmetry plane, replicates the trend observed in Dasi et al. (Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007) by PIV measurements and in Yun et al. (Reference Yun, Dasi, Aidun and Yoganathan2014a) by means of a lattice-Boltzmann computational approach. Actually, the mean flow field issuing from our simulations matches the maximum flow velocity and the width of the shear layer, whereas the backflow part of the profile is not represented. This is due to the significant differences in the geometric modelling of the aortic sinuses and in the valve orientation: it is clearly evident from figure 5(a) that the $y$-slice (centred on the valve axis) does not include recirculation regions. Analogous considerations can be drawn for the velocity root-mean-square (RMS) profile, bearing also in mind that different peak Reynolds numbers are involved in the comparison.

Figure 5.

Dimensionless axial component $\bar {w}/U$ of the mean flow field at subsequent cross-sections of the aortic root: (a) PF, (b) ED phase, (c) LD phase.

Figure 6.

Streamwise distribution of phase-averaged pressure drop along the valve axis. (a) Peak flow rate. (b) Early deceleration.

The LTFV model instead provides four jets with a wider central jet. At this stage the recirculation regions reach the maximum intensity in all configurations, with axial velocity peaks of $\bar {w}/U \approx -0.6$ in the SJMV and SJMVVG configurations, and $\bar {w}/U \approx -0.45$ in the LTFV. The reverse flow regions are observed to be localised in correspondence of the Valsalva sinuses in the bileaflet configurations, with the greatest intensity in the non-bisected sinus, in accordance with the prediction by Akutsu, Matsumoto & Takahashi (Reference Akutsu, Matsumoto and Takahashi2011). Conversely, the LTFV model provides a significantly lower amount of recirculation, being uniformly distributed in the near-wall regions of the aortic sinuses. Such a feature is a consequence of the installation of the LTFV leaflets within the housing (the reader should refer to figure 1), which puts the leaflets in a further forward position with respect to the aortic sinuses, when compared with the SJMV configuration. Thus, in the fully open position the LTFV leaflets shield the sinus cavity from the flow to create extended recirculation regions with large negative bulk flow velocity. This difference can be also perceived from the instantaneous vorticity fields over the valve longitudinal sections plotted in figures 7 and 8. At PF the LTFV jets propagate up to the streamwise position $z/D=2.2$ while holding their coherency. At ED the periodic vortex shedding triggered by the inflow rate reduction and by the development of the adverse pressure gradient in all models can be clearly identified in the instantaneous vorticity contours depicted in figure 8. In the ED stage the bileaflet valve jets extend up to the coordinate $z/D=3.5$ with reduced bulk velocity with respect to PF, while the characteristic triple jet pattern and the recirculation regions are maintained in space up to the coordinate $z/D=2.8$. In the bileaflet design the side jets impinge the aortic root generating a thin wall shear layer, whereas the trileaflet arrangement generates narrower lateral jets, symmetrically distributed in the azimuthal direction, which in turn provide a thicker wall shear layer. In section $z/D=2.5$ the central LTFV jet achieves a slightly larger bulk velocity ($\bar {w}/U \approx 0.90$) than the bileaflet lateral jets ($\bar {w}/U \approx 0.86$). It is worth pointing out that at ED the LTFV configuration shows significant changes in the topology of the mean flow between cross-sections $z/D=1.5$ and $z/D=2.5$. In the space of just one diameter the lateral jets disperse their momentum between the central jet and three circumferential jets generated from the interaction of the former with the structures shed from the housing. Again, a better understanding of the flow topology can be achieved comparing the mean flow (figure 5) with the instantaneous vorticity field (figure 8). In the LD phase the lateral jets of the SJMV model are still well separated in the aortic root, whereas the central jet starts to disintegrate due to the transition to turbulence. Instead, the circumferential jets observed at ED in the LTFV mix with the core jet to generate a three-lobed jet topology with $120^\circ$ symmetry. In section $z/D=2.5$ the bulk velocity of such a structure is almost unchanged from previous cycle instants, whereas the bulk velocity of bileaflet lateral jets has reduced to $\bar {w}/U \approx 0.79$. The same structure quickly breaks up as the jet advances downstream ($z/D=3.5$), introducing increased stochasticity in the flow behaviour. Immediately downstream of the leaflet both bileaflet cases produce strong reverse flow regions within the sinus cavity. On the contrary, the LTFV presents negative flow rate further downstream of the leaflets ($z/D=2.5$), which are not present in the other configurations.

Figure 7.

Contours of instantaneous out-of-plane vorticity ($\omega _x D/U$ and $\omega _y D/U$) over the longitudinal cross-sections of the aortic root in the PF stage.

Figure 8.

Contours of instantaneous out-of-plane vorticity ($\omega _x D/U$ and $\omega _y D/U$) over the longitudinal cross-sections of the aortic root in the ED stage.

The axial distribution of the mean pressure drop $(\bar {p}-\overline {p_0})/(\rho U^2)$ along the valve axis is plotted in figure 6 at PF and ED, where $\overline {p_0}$ indicates the mean pressure at the inlet section. At PF, the transvalvular pressure drop (TPD) is assessed, corresponding to 0.629 $\rho U^2$ (5.04 mmHg) for the SJMV model and 0.667 $\rho U^2$ (5.35 mmHg) for the SJMVVG model. Conversely, the trileaflet valve provides a much smaller pressure drop, 0.390 $\rho U^2$ (3.12 mmHg), owing to the presence of a wide central orifice. The TPD value found for the SJMV at PF matches the systolic gradient measured one year post-surgery, equal to 5.7 mmHg (St Jude Medical, Inc. 2010). No benchmark value for the TPD within a trileaflet valve is found. We found that LTFV is characterized by a smoother pressure recovery crossing the valve, while VGs bring a slight delay in the pressure recovery and negligible changes in the mean pressure drop. The minimal pressure increment with VGs is due to the further reduction of the orifice area.

It is worth noting that in our simulations we prescribe the same time-varying flow rate at the inflow boundary, and a radiative outflow condition at the outlet boundary for all valve models, whereas in the experimental set-up of Hatoum & Dasi (Reference Hatoum and Dasi2019), the flow rate emerges from the coupled problem, as a consequence of the pressure settling within the mock loop. Pressure recovery also depends critically on the root geometry that in our case includes the three sinuses of Valsalva, while in the experimental case it was straight (cylindrical). Although the VGs entail some differences in the mean flow field we emphasise that the cross-section contours in figure 5, as well as any other cuts of the mean flow field, do not provide a fully clear picture of the effect of the VGs on the flow field. Thus, we can speculate that they provide barely discernible differences on the phase-averaged mean flow, and we rely on instantaneous fields and second-order statistics to figure out the influence of these geometric appendices.

Ultimately, the observed mean flow fields contain some undamped vortical structures due to the relative brevity of the sampling time window, which is not wide enough to encompass the full characteristic period of the largest Kármán-like dipoles shed from the obstructions.

3.2. Fluctuating stress analysis

We illustrate stress distributions through the different valve configurations. Once acknowledged the ubiquitous nature of blood damaging mechanisms in turbulent flows (briefly mentioned in the introductory part of the manuscript), we limit our speculations to an Eulerian perspective of turbulent stresses. Reynolds stresses arise from the chaotic fluctuations primarily induced by the turbulent transition, and represent a fictitious stress contribution correlated with an enhanced blood damaging (Faghih & Sharp Reference Faghih and Sharp2019). A widely accepted conjecture postulates that abnormal flow conditions such as regurgitant jets and turbulent flows provide an instantaneous increase in viscous shear stress, exposure time and phase collision probability (Leverett et al. Reference Leverett, Hellums, Alfrey and Lynch1972; Antiga & Steinman Reference Antiga and Steinman2009; de Tullio et al. Reference de Tullio, Nam, Pascazio, Balaras and Verzicco2012). Hence, the evaluation of phase-averaged fluctuating stresses certainly provides a heuristic estimation of the stress levels experienced by blood elements, as well as the topology of the high-damaging regions. Although the fluctuating viscous stresses have been regarded as a far more appropriate metric than Reynolds stresses for quantifying blood damage (Ge et al. Reference Ge, Dasi, Sotiropoulos and Yoganathan2008), we argue that in a phase-averaging environment free from cycle-to-cycle variations, the latter provide a picture of the stochastic convective acceleration causing instantaneous enhancement in the viscous peak stresses and residence time experienced by blood constituents.

In the biomechanical literature the levels of both viscous and Reynolds stresses are often compared with haemolysis and platelet activation thresholds which have been obtained from a simple shear or extensional stress condition in laminar flows. The extension of these thresholds to turbulent stresses and complex stress states is fairly questionable, therefore, we interpret these thresholds as merely indicative.

A significant advantage of a numerical approach with respect to experimental investigations lies on the fact that three-dimensional evaluation of turbulent stress components can be performed without recurring to two-dimensional simplifying assumptions. In the experimental practice, the haemodynamic compliance of prosthetic heart valves is usually assessed from velocity information extracted over significant planes. However, the flow through a mechanical heart valve with realistic aortic root geometry was extensively proved to be fully three dimensional (Sotiropoulos & Borazjani Reference Sotiropoulos and Borazjani2009), as well as the turbulent stress distribution. Two-dimensional analysis has been proved to underestimate the normal and shear stresses, depending on the valve model and aortic root geometry (de Tullio et al. Reference de Tullio, Cristallo, Balaras and Verzicco2009). Since the peak Reynolds stress is not invariant under coordinate system rotation, a coordinate-independent metric is needed to quantify mechanical loads. Thus, we evaluate the full Reynolds stress tensor for every Eulerian node within the aortic root and compute the maximum turbulent shear stress (TSS) acting on a surface element (Malvern Reference Malvern1969) as the principal shear stress component

(3.1)\begin{equation} {\rm TSS}_{max}=\tfrac{1}{2} \left(\sigma_{max}-\sigma_{min} \right), \end{equation}

where $\sigma _{max}$ and $\sigma _{min}$ represent the maximum and minimum principal normal stress components. It is recognised that flow through mechanical heart valves exhibits principal Reynolds normal stresses twice larger than the principal shear Reynolds stresses (Liu et al. Reference Liu, Lu and Chu2000; Nyboe et al. Reference Nyboe, Funder, Smerup, Nygaard and Hasenkam2006). However, the normal stresses are not regarded in this section since they are believed less harmful to RBC membrane (Faghih & Sharp Reference Faghih and Sharp2019). Here $\sigma _{max}$ and $\sigma _{min}$ can be obtained by solving the characteristic polynomial equation associated with the Reynolds stress tensor $\overline {u'_i u'_j}$,

(3.2)\begin{equation} \sigma^3-I_1 \sigma^2 + I_2 \sigma -I_3 =0, \end{equation}

being $I_1$, $I_2$ and $I_3$ the second-order tensor invariants

(3.3)\begin{equation} \left.\begin{gathered} I_1 = \overline{u'_i u'_i}, \\ I_2 = \tfrac{1}{2} \left[\left(\overline{u'_i u'_i} \right)^2 - \overline{u'_i u'_i}^2 \right] , \\ I_3 = \tfrac{1}{6} \left(\overline{u'_i u'_i} \right)^3 - \tfrac{1}{2} \overline{u'_i u'_i}^2 \left(\overline{u'_j u'_j} \right)^2 +\tfrac{1}{3} \overline{u'_i u'_i}^3 . \end{gathered}\right\} \end{equation}

In order to provide a fair comparison among different valve models and different cycle instants we restrict the evaluation of the $\textrm {TSS}_{max}$ to the interior node cloud from $z/D=0.8$ to $z/D=4.0$. Since the leaflet kinematics is prescribed, the velocity fluctuations at the solid–fluid interface are negligible. Before discussing the $\textrm {TSS}_{max}$ distributions, we highlight that the usefulness of the principal stress transformation to characterise fluid stresses for blood damage is uncertain and may depend on the local flow configuration, nevertheless, the complexity of the flow though prosthetic devices enables this descriptor as the most representative item for macro-scale evaluations.

The comparison carried out in figure 9 confirms that the level of accuracy achieved in our simulation guarantees a satisfactory agreement with experimental data both on first-order and second-order statistics. Figures 11–13 display two cutaway planes (corresponding to the spanwise planes $XZ$ and $YZ$) of the three cases investigated with $\textrm {TSS}_{max}$ contours, along with probability density distributions of $\textrm {TSS}_{max}$, computed from three-dimensional fields. The interpretation of these figures can be facilitated by inspecting the streamwise distribution of turbulent kinetic energy (TKE), defined as $(\overline {u'_k u'_k})/U^2$, which is collected computing a RMS value at each cross-section (see figure 10). Previous experimental investigations on the SJMV model (Liu et al. Reference Liu, Lu and Chu2000; Hatoum et al. Reference Hatoum, Maureira and Dasi2020) found the peak systole to be characterised by the highest magnitude of $\textrm {TSS}_{max}$ within the cardiac cycle, whereas we identify a more complicated scenario. Certainly, the largest TKE RMS values are encountered at ED, whereas the phase location of the largest TSS$_{max}$ depends on the valve configuration. Such a discrepancy might be motivated by the difference in the boundary conditions governing the flow field, which in the experimental procedures consists of the pressure gradient applied by the hydrodynamic loop. At ED the TKE RMS peaks are found one diameter downstream of the leaflets for both SJMV and SJMVVG (see figure 10), whereas the LTFV provided similar peak values at PF and ED. In the PF SJMV configuration, the peak $\textrm {TSS}_{max}$ value, corresponding to a $10^{-5}$ probability density value, is found to be $(\textrm {TSS}_{max} /\rho U^2)_{P=1.0\times 10^{-5}}=0.142$, i.e. 151.4 Pa. The presence of VGs halves the same quantity to 0.0761 (81.15 Pa), whereas the LTFV generates similar values, namely 0.134 (142.9 Pa). Such values refer to a very limited event probability since they are confined to narrow spots in the flow field, located immediately downstream of the leaflets. A greater relevance in blood damage sourcing can be attributed to the volumetric median, whose values take 1.038, 0.479 and 0.294 Pa for the SJMV, SJMVVG and LTFV cases, respectively. This matches the experimental findings of Hatoum & Dasi (Reference Hatoum and Dasi2019) who correlated the reduction of the turbulent stress level in the co-rotating VGs arrangement with the reduction of vorticity intensity induced by the corruption of the most energetic coherent structures shed from the leaflet edges. This effect is captured in the $\textrm {TSS}_{max}$ contours in figure 12, where the VGs clearly mitigate the intensity of the stress spots connected to the large-scale vortex shedding, and it is reflected in a lower and spatially delayed TKE peak. In light of these findings, the additional instabilities brought by VGs seem to reduce the intensity of both normal and shear second-order fluctuations. The mechanism beyond such stress distribution and the spatial scales of the eddies that contribute to their production will be reinforced in the section concerning the spectral analysis. The LTFV configuration induces a similar peak value to SJMV, but a much lower median value, probably because the additional leaflet increases the probability of large stress levels. On the other hand, the limited mixing occurring in the LTFV entails narrower low-stress regions, mostly localised in the near-wall shear layer.

Figure 9.

Profile of mean- and RMS value of the phase-averaged axial velocity component over the SJMV symmetry line at distance $z/D=0.22$ from the leaflet tip, compared with data from literature, at PF.

Figure 10.

Streamwise distribution of the cross-sectional RMS value of turbulent kinetic energy.

Figure 11.

Contours of $\textrm {TSS}_{max}$ over two spanwise cross-sections of the aortic root at the PF stage and probability density distribution computed from the corresponding three-dimensional field.

Figure 12.

Contours of $\textrm {TSS}_{max}$ over two spanwise cross-sections of the aortic root at the ED stage and probability density distribution computed from the corresponding three-dimensional field.

At ED, we find $\textrm {TSS}_{max}$ peak values equal $0.131 \rho U^2$ (139.7 Pa) for the SJMV, and $0.119 \rho U^2$ (126.9 Pa) for the SJMVVG, mostly distributed in the leaflet wake and in the near-wall region at the sino-tubular junction. At this stage, where the vortex shedding from leaflet tips is fully developed, the stress level up to $z/D=1.5$ is much lower in the SJMVVG case than the SJMV case, due to the influence of VGs. Conversely, the LTFV provides a much lower peak stress level, with $(\textrm {TSS}_{max} / \rho U^2)_{P=1.0\times 10^{-5}}=0.0907$ (96.72 Pa), distributed on the edges of the core jet, which still preserve full coherency. This value is significantly larger than what was found by Hatoum et al. (Reference Hatoum, Maureira and Dasi2020), who denoted a $\textrm {TSS}_{max}$ peak of 72 Pa for the SJMV model, and approximately close to the values reported by Ge et al. (Reference Ge, Dasi, Sotiropoulos and Yoganathan2008) from three-dimensional numerical investigations. We stress that experimental evaluations based on two-dimensional velocity measurements provide systematically lower stress estimates. Differently, the computational investigation by Ge et al. (Reference Ge, Dasi, Sotiropoulos and Yoganathan2008) led to peak $\textrm {TSS}_{max}$ values close to 100 Pa. When comparing these values with the incipient haemolysis threshold brought to light from experimental correlations by Jhun et al. (Reference Jhun, Stauffer, Reibson, Yeager, Newswanger, Taylor, Manning, Weiss and Rosenberg2018), equal to 3000 Pa at $Re \ge 6 \times 10^4$, one might argue that the level of turbulent stress exerted on blood constituents across a BMAV are far below the critical value. However, they are recognised to be sufficient to induce significant sublethal damage to RBCs and platelet lysis for a given exposure time. In a representative study by Sutera (Reference Sutera1977), platelet clustering and lysis under shear load was found by simple shear tests for stress equal to 10 Pa. The transition to turbulence allows the bedding of the median stress level to 2.026 Pa for SJMV and SJMVVG, and to 0.906 Pa for the LTFV. The raise of the median values from PF to LD in all configurations bears witness to the transition to a distributed turbulent state across the computational domain. As the inflow rate decreases and a more uniform level of turbulent production occupies the downstream domain region, a larger probability of getting near-median stress values is observed, whereas the peak values are significantly reduced. At LD the peak and median stress values among the investigated cases conform to $0.0496 \rho U^2$ (52.89 Pa) and $0.00545 \rho U^2$ (5.81 Pa), respectively. It is worth noting that the LTFV provide a displaced TKE peak, at $z/D=2.7$, with respect to the other configurations ($z/D=1.9$), possibly caused by the fluctuations occurring in the formation of the tri-lobed jet observed in the mean flow at $z/D=2.5$ (the reader should refer to figure 5).

A statistical hypothesis testing is performed on the probability density distributions in figures 11–13, comparing SJMVVG and LTFV against the SJMV. The chi-square test provided $p<0.001$, confirming that the null hypothesis is incorrect, therefore, the differences among the distributions can be considered statistically significant. The same outcome was obtained for all other distributions provided in the present work.

The picture of turbulence-induced flow environment experienced by RBCs can be completed by evaluating the fluctuating viscous stresses. These can be interpreted as the actual stress exerted on the volume element due to phase-averaged strain fluctuations. Following Jones (Reference Jones1995) the viscous stress distribution is evaluated by means of an equivalent scalar field $\tau _{eq}$, named equivalent viscous stress (EVS) computed by double contraction of the fluctuating strain-rate tensor $s_{ij}$,

(3.4)\begin{equation} \tau_{eq}= \mu \sqrt{2 \overline{s_{ij} s_{ij}}}, \end{equation}

with $s_{ij}=1/2 (\partial u'_i/\partial x_j + \partial u'_j/\partial x_i)$. Despite the fact that this information is deprived of the directionality of the stress, we rely on it because it comprises all the strain components, and because the majority of the stress-based models for haemolysis are built out of a scalar measure (Faghih & Sharp Reference Faghih and Sharp2019). Stress-based models are acknowledged to overestimate the damaging level, thus providing conservative estimates. Furthermore, the damaging level of blood constituents in high shear flows is affected by their degree of alignment with the flow direction (Goldsmith, Marlow & MacIntosh Reference Goldsmith, Marlow and MacIntosh1972) rather than the direction of the stress principal axis. Thus, we consider $\tau _{eq}$ a relevant and sufficiently general descriptor. It is worth noting that the EVS corresponds to the TKE dissipation rate $\varepsilon =2 \nu \overline {s_{ij} s_{ij}}$, by means of the relation $\tau _{eq}=\rho \sqrt {\nu \varepsilon }$. Thus, the EVS also depicts the rate at which the turbulence energy is dissipated by viscous forces.

At PF fluctuating viscous stresses are localised in the near-leaflet wake, extending $1.2D$ downstream of the trailing edge, as previously acknowledged by Yun et al. (Reference Yun, Dasi, Aidun and Yoganathan2014b) by inspecting the TKE dissipation rate contours. Likewise in simplified aortic root geometries (Ge et al. Reference Ge, Dasi, Sotiropoulos and Yoganathan2008; Yun et al. Reference Yun, Dasi, Aidun and Yoganathan2014b), a second viscous stress peak is found in correspondence of the housing edge, where the diffusion of the shear layer generates recirculation eddies constrained within the sinus cavity by the lateral jet. We denote from figure 14 the presence of these spots in all the configurations considered. Peak stress levels are identified with the same criteria used for $\textrm {TSS}_{max}$. However, we speculate that the intensity of the dissipation in the Valsalva sinuses strongly depends on the type of valve implantation. The LTFV presents the lowest EVS peak, corresponding to $0.00683 \rho U^2$ (7.28 Pa), whereas in the other configurations the EVS peak takes the approximate value $0.00797 \rho U^2$ (8.50 Pa). Since the EVS is negligible within most of the aortic root, the median values are lower than $10^{-4} \rho U^2$ (0.107 Pa). The flow mechanisms characterising the ED make the peak stress value raise in all configurations, although larger increments are recorded for bileaflet cases. Specifically, $0.0106 \rho U^2$ (11.3 Pa), $0.0105 \rho U^2$ (11.19 Pa) and $0.00763 \rho U^2$ (8.13 Pa) for SJMV, SJMVVG and LTFV, respectively. In this phase the bileaflet models provide identical probability of the $\tau _{eq}$ event, although the spanwise contours in figure 15 show a larger shear stress in the early wake region in the SJMV. Once again this difference can be correlated with the reduced intensity of turbulent fluctuations brought by the VG arrays. The LTFV configuration generates definitively less viscous dissipation both in terms of peak and median value, mostly connected with the mixing of side jets and the core jet also noticed in the mean flow field (see figure 5).

Figure 13.

Contours of $\textrm {TSS}_{max}$ over two spanwise cross-sections of the aortic root at the LD stage and probability density distribution computed from the corresponding three-dimensional field.

Figure 14.

Contours of $\tau _{eq}$ over two spanwise cross-sections of the aortic root at the PF stage and probability density distribution computed from the corresponding three-dimensional field.

Figure 15.

Contours of $\tau _{eq}$ over two spanwise cross-sections of the aortic root at the ED stage and probability density distribution computed from the corresponding three-dimensional field.

Figure 16.

Contours of $\tau _{eq}$ over two spanwise cross-sections of the aortic root at the LD stage and probability density distribution computed from the corresponding three-dimensional field.

Throughout the cycle phases characterised by the largest TKE, in the leaflet wake, the maximum shear stress level falls well below the haemolysis threshold previously discussed in the turbulent shear stress assessment. On the other hand, the $\tau _{eq}$ values recorded for the bileaflet cases are beyond the platelet activation threshold measured in Hung et al. (Reference Hung, Hochmuth, Joist and Sutera1976), corresponding to $10.0 \div 16.5$ Pa with an exposure time of 102 s, and in Ramstack, Zuckerman & Mockros (Reference Ramstack, Zuckerman and Mockros1979) corresponding to $30 \div 100$ Pa with an averaged residence time of 10 s. If relying on these thresholds, the LTFV can be regarded as less harmful for platelet activation. Eventually, the peak viscous stresses drop in the LD phase to the values 4.52 Pa for bileaflet cases and to 4.06 Pa for the LTFV (see figure 16).

An evaluation of the mean viscous shear stress could be useful to provide a broader picture of the flow stresses, and their correlation with the blood damage. However, these are not investigated since this study is mainly focused on the turbulence-induced stresses.

In view of the complex distribution of turbulent stresses, an integral interpretation would be impractical without considering more than one scalar quantity (Banerjee et al. Reference Banerjee, Krahl, Durst and Zenger2007), despite the fact they mostly act at a larger scale than the characteristic size of RBC (Antiga & Steinman Reference Antiga and Steinman2009). This is especially true in such a confined flow, at nearly transitional Reynolds number, where the uncertainty about the collocation of Reynolds stresses in the wavenumber space is much higher than simple isotropic homogeneous flows. To this extent, we quantify the level of anisotropy of turbulent stresses at ED and LD phases, where the transition to turbulence is definitely achieved.

The turbulent stress anisotropy is assessed by a colour mapping of the traceless anisotropy tensor $B_{ij}$, defined by

(3.5)\begin{equation} B_{ij}=\frac{\overline{u'_i u'_j}}{\overline{u'_k u'_k}} - \frac{\delta_{ij}}{3}, \end{equation}

with $\delta _{ij}$ being the Kroneker delta. The magnitude and orientation of turbulent fluctuations do conform, in a time-averaged sense, to specific states that can be characterised by a suitable mapping of the anisotropy tensor eigenvalues. Following Lumley & Newman (Reference Lumley and Newman1977), we explore the different turbulent states by building, at some relevant cross-sections, the anisotropy invariant map (AIM), whose coordinates are obtained as

(3.6a–c)\begin{equation} C_{1C}=\lambda_1-\lambda_2, \quad C_{2C}= 2 \left( \lambda_2-\lambda_3 \right), \quad C_{3C}= 3 \lambda_3 + 1. \end{equation}

Here $\lambda _i$ are the eigenvalues of $B_{ij}$, sorted in descending order, and $C_{iC} \ge 0$ are the eigenvalue-related weights of $B_{ij}$ in the appropriate tensor basis (Emory & Iaccarino Reference Emory and Iaccarino2014), making up the three turbulence corner states. To make broad-scale physical interpretations easier in our framework, we represent the AIM distribution at subsequent cross-sections of the aortic root by means of a contour representation whose scalar values are obtained assigning the RGB (Red-Green-Blue) colour triplet $[C_{1C},C_{2C},C_{3C}]$. This representation was proposed by Emory & Iaccarino (Reference Emory and Iaccarino2014), and named anisotropy componentality contours. Each component of the colour triplet represents the contribution to a AIM corner state (see figure 17), with corresponding stress shape: one-component turbulence (1C), two-component axisymmetric turbulence (2C) and isotropic turbulence (3C). The three borders joining the AIM corners provide intermediate limiting stages, described in figure 17. The global level of anisotropy of turbulent fluctuations at a certain cross-section is simply denoted by the anisotropy index, defined as

(3.7)\begin{equation} AI=\sqrt{\frac{1}{n} \sum_{i=1}^{n} (C_{1C}+C_{2C})_i^2} \end{equation}

over the $n$ nodal locations within the observation domain.

At ED the turbulent stresses show a high anisotropy level, with $AI>0.75$ in section $z/D=2.6$, where the stresses mostly approach the one-component limit, both in the SJMV and LTFV configurations (see figure 18). At ED the cross-sections corresponding to $z/D=3.2$ and $z/D=3.8$ cannot be considered meaningful for this purpose because they do not lay in a region where significant turbulent fluctuations occur, so the related contours refer to the residual fluctuation being washed off the domain. In section $z/D=2.0$, which is pervaded by intense vorticity fluctuations at ED (see figure 8), the anisotropy level settles to $AI=0.61$ for SJMV and to $AI=0.74$ for LTFV. The former configuration presents a wide range of turbulent stress states uniformly distributed over the section, despite the presence of large shear layers in the mean flow field, depicted in figure 5. The latter instead provides a central region approaching the two-component limit, which makes the anisotropy index larger. At LD the anisotropy index is found to decrease and to become uniform across the domain to $AI \approx 0.6$ for both cases. In sections $z/D=2.6$ and $z/D=3.2$ the most elevated anisotropy occurred in correspondence of the shear layer surrounding the central jet, where the stress states are biased towards the one-component limit (1C), and the axisymmetric expansion bound. In the near-wall region, the stress states mostly approach the two-component limit, entailing a combination of red/green regions, as expected in the presence of a wall (Mansour, Kim & Moin Reference Mansour, Kim and Moin1988). The contours regarding the SJMVVG configuration are not reported here for the sake of brevity, since they provide a pattern very similar to that observed for the SJMV case. We speculate that despite the different vortex dynamics induced by the VGs, they do not affect the Reynolds stress directionality.

The considerations listed above on stress anisotropy provide a useful perspective for interpretation of spectral analyses, which are often carried out with the hypothesis of homogeneous isotropic flow.

3.3. Spectral analysis

In order to correlate the reduction in the turbulent stresses observed with the VGs to the large-scale vortex dynamics, we investigate the time spectra of instantaneous kinetic energy $( u_i )^2/U^2$ at some points aligned with the SJMV leaflets in the fully open position. The LTFV configuration is excluded from this comparison. Figure 19 shows the time spectra of kinetic energy at the axial coordinates $z/D=1.5$ and $z/D=3.0$, averaged over the monitor lines 3, 5, 7, 9. Such spectra are computed from the velocity signal within a much larger cycle window than those considered for phase averaging in order to include the extended vortex shedding dynamics from PF $Q/Q_{max}=1.0$ to mid-deceleration $Q/Q_{max}=0.5$. In the leaflet shear layer the flow exhibits strong vortex shedding, with high energy content in the smaller temporal frequencies, whose dynamics is portrayed in figure 8 at ED. As the jets propagate further downstream, and the shear layers breaks up in favour of a mixing dynamics, the energy contained in the smaller frequency reduces by one order of magnitude. We emphasise that within the time window considered for signal sampling, the classical vortex shedding observed downstream of leaflets (Dasi et al. Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007) develops for both SJMV and SJMVG models, but the presence of the VGs leads to a lower energy content at all time scales, confirming that the associated co-rotating vortices locally disrupt the vortex rings periodically shed from the leaflet trailing edge in the high flow rate stage of the systole. Although this observation is provided from local evidences, we postulate that a similar pattern is consistently replicated in the leaflet tangent direction, explaining the significant reduction in the probability of large $\textrm {TSS}_{max}$ events at PF (see figure 11). This trend is reversed at the location $z/D=3.0$, where the SJMVVG provides larger kinetic energy, within a region were full transition occurred for both models. Further analysis on the energetic content of the purely fluctuating components are provided in the following by means of spatial spectra. One of the major motivations of this analysis is to correlate the turbulent stresses at different phases of the cardiac cycle with the geometry of the mechanical prostheses. Once the anisotropy pattern and the dissipative features occurring in the different configurations have been observed, the aforementioned correlations are sought by further investigating the spatial scales characterising the turbulent fluctuations in the frequency domain.

Figure 17.

Barycentric colourmap of turbulent stress anisotropy.

Figure 18.

Anisotropy componentality contours of the SJMV and LTFV configurations at the ED and LD stages, in subsequent aortic root cross-sections.

Figure 19.

Time spectra of the kinetic energy at the axial coordinates (a) $z/D=1.5$ and (b) $z/D=3.0$, averaged over the monitor lines 3, 5, 7 and 9.

Provided the strongly anisotropic nature of turbulent fluctuations we consider a local characterisation of the spatial fluctuations to be of greater impact in this application, therefore, power spectra are computed on the monitor lines (whose configuration is illustrated in figure 4), over segments starting downstream of the leaflets. A Hann window is applied to the velocity signals and the spatial average along the $z$-direction is also subtracted. Line-averaged spectra are plotted for the streamwise TKE (figure 20), as well as for the anisotropic streamwise/leaflet-normal energy (figure 20) at three subsequent cycle instants. It can be immediately grasped that the energetic content of the turbulent fluctuations increases, especially at the largest spatial scales, as the inflow rate decelerates by the progressive shift of all curves towards higher energy levels. As the incoming flow rate decreases, the shear layers surrounding the core jet become more unstable due to the adverse pressure gradient, triggering a larger energy release towards smaller flow structures. None of the plotted spectra provide a clear classic scaling (Pope Reference Pope2001), despite the fact that a good fit with the inertial range scaling ($E_{ww}(k) \propto k^{-5/3}$) was identified already at the PF in previous studies (Yun et al. Reference Yun, Dasi, Aidun and Yoganathan2014b). In the present study the spatial spectral analysis is performed on localised lines within an anisotropic flow, which are not necessarily located where turbulence develops first or it is the strongest. Furthermore, one should bear in mind that at the PF instant most of the aortic root is crossed by residual low-energy structures that are immediately washed out in the acceleration phase. The recent experimental and numerical investigations by Cerbus et al. (Reference Cerbus, Liu, Gioia and Chakraborty2020), providing insights on local energy spectra for transitional pipe flow, reported that the power spectra (on-centreline and off-centreline) conforms to $-5/3$ scaling only for $Re \geq 80\,000$. Since the leaflets introduce a stronger degree of anisotropy in the flow than the simple free pipe, and the Reynolds number conforms to the transitional pipe regime, there is no reason to expect the Kolmogorov universal scaling laws.

Figure 20.

Streamwise isotropic energy spectra at different cycle instants. Solid and dashed lines correspond to spectral values averaged over the inner and outer monitor lines annuli.

When considering the LTFV configuration, the inner line annulus fall within the central jet, which holds spatial coherency longer than the central jet ejected from the SJMV/SJMVVG models. As a matter of fact, the $E_{ww}(k)$ spectrum mediated over the inner annulus has significantly lower energy content at the PF with respect to the outer annulus of the same model (see figure 20). At the ED/LD instants the energy distribution appears to be more uniform among the different configurations due to the enhanced mixing occurring in the deceleration phase. We emphasise that the intrinsic fluctuation associated with the energy-containing scales, corresponding to the largest wavenumber, at the ED/LD instants provides larger energy injection in the LTFV case, for both the inner and outer line annulus, owing to a more intense breakdown of the leaflet shear layer. These considerations apply both to the isotropic spectra and the shear spectra plotted in figure 21. We highlight that, differently from the isotropic spectra, the shear spectra, provide, especially in the LD phase, a better matching with the Lumley's scaling law, by which the small level of anisotropy in the inertial subrange can be assumed to be a small perturbation of the background isotropic state (Pope Reference Pope2001). This confirms that the strongly anisotropic nature of the fluctuations is lost in the central wavenumber range at LD, whereas a significant level of anisotropy characterises the intermediate length scales, and not only the energy-containing motion, at the PF and ED stages. It is worth pointing out the absence of energy pileup at high wavenumbers in any spectral representation, and the smoothness of the energy cascade which covers approximately nine orders of magnitude, meaning that most of the smallest flow scales are well represented without energy accumulation.

Figure 21.

Anisotropic energy spectra at different cycle instants. Solid and dashed lines correspond to spectral values averaged over the inner and outer monitor lines annuli.

Our interpretation of the local energy spectra does not contradict the computational findings of Yun et al. (Reference Yun, Dasi, Aidun and Yoganathan2014b) and Becsek, Pietrasanta & Obrist (Reference Becsek, Pietrasanta and Obrist2020), that identified a clear inertial range from volume-averaged spectra collected over a wider flow rate interval. Since no clear indications on the fundamental micro-scales can be inferred from the energy spectra, we further inspect the nature of the fundamental scales involved in the flow field by computing the longitudinal Taylor micro-scale (LTMS) over the monitor lines. The Taylor micro-scale, denoted by $\lambda$, represents the characteristic length of the largest eddies which are responsible for the dissipation of energy (Pope Reference Pope2001). Since the RBC membrane strain, which is a precursor to membrane failure, has a logical connection to surrounding fluid strain, the estimation of $\lambda$ might be relevant to haemolysis beyond any specific haemolysis mechanism (Faghih & Sharp Reference Faghih and Sharp2019). Longitudinal Taylor micro-scale is associated to the curvature of the velocity autocorrelation function. Thus, the longitudinal two-point correlation is computed for each line, at the generic streamwise location $z_0$ as

(3.8)\begin{equation} R_{ww}(z_0,\Delta z)= \overline{w'\left(z_0-\Delta z/2\right)w' \left(z_0+\Delta z/2\right)}. \end{equation}

This relation, once rescaled with respect to the velocity variance, provides the longitudinal autocorrelation function $f_a(z_0,\Delta z)$. Then, for small separation $\Delta z$, the LTMS scale $\lambda (z_0)$ can be approximated by the roots of the osculating parabola of the autocorrelation function (Pope Reference Pope2001),

(3.9)\begin{equation} \lambda(z_0)=[-\tfrac{1}{2} f_a''(z_0,0)]^{-({1}/{2})}.\end{equation}

The second derivative in (3.9) is computed using a sixth-order accurate centred finite-difference scheme, using a seven-point stencil. This construction is carried out exploiting the symmetric correlation functions at different spots in between the streamwise coordinates $z/D=2.0$ and $z/D=4.0$. Such a construction provides physically meaningful results for small values of $\lambda /\varLambda$ in homogeneous flows, with $\varLambda$ being the integral length scale (Belmabrouk & Michard Reference Belmabrouk and Michard1998). The influence of the additional terms appearing in higher-order expansions for non-homogeneous turbulence can be neglected for small $\Delta z$, namely, for a small deviation of the autocorrelation function from the osculating parabola (Belmabrouk & Michard Reference Belmabrouk and Michard1998).

In the present framework the evaluation of the LTMS does provide a cross-check on the dimension of the largest dissipative length scale with respect to energy spectra. Furthermore, we consider the LTMS to accomplish a more reliable estimation on fundamental scales than the Kolmogorov length scale, since the size of the smallest eddies might be affected by the local grid resolution. In contrast, the smallest LTMS observed, corresponding to $\lambda /D=0.03$, can be overlaid onto approximately nine grid points, meaning that the corresponding eddy can be fully resolved with a second-order scheme. The axial coordinates $z_0$ from which the correlation functions are computed are chosen so that the zero-correlation condition is reached in a region of high vorticity. A polar map of the LTMS is reported in figure 22 for the three valve models analysed at the ED stage. In the SJMV and SJMVVG cases we identified a mean LTMS scale value of $\lambda /D=0.045$, corresponding to 1.0 mm. This value is approximately close to the length scale found by Yun et al. (Reference Yun, Dasi, Aidun and Yoganathan2014b), who identified the start of the energy decay at $l=3$ mm, from visual inspection of the energy spectra. Liu et al. (Reference Liu, Lu and Chu2000) obtained an experimental estimation of the LTMS, falling into the interval $\lambda =0.1\text{--} 0.6$ mm, computed by means of the velocity–time autocorrelation function, which was integrated from laser-Doppler anemometer measurements. They extracted the spatial LTMS approximating the spatial correlations by temporal correlations in view of the frozen field approximation (Taylor Reference Taylor1938). However, in pulsatile flows with a high degree of anisotropy this approximation is generally considered inaccurate (Pope Reference Pope2001). The LTFV configuration presented a larger LTMS, with a mean value of $\lambda /D=0.075$ at LD (see figure 23), corresponding to 1.73 mm. A larger LTMS value reflects the lower TKE RMS profile in figure 10, confirming that the trileaflet model produces a less energetic turbulence decay, with a larger dissipative structure, which can be thought to be less harmful to RBCs.

Figure 22.

Polar map of the LTMS $\lambda /D$ in the ED phase, at two different coordinates.

The mean LTMS value found from our analysis corresponds to the non-dimensional wavelength $k/D \approx 130$, corresponding to a full energy decay region. No clear differences are generated by VGs, meaning that the different mixing mechanism at the energy-containing scales does not affect the dissipation in the smaller flow structures. The LTFV configuration provides highly non-homogeneous values, probably because the breaking of the shear layer carrying the jets mixing takes place further downstream with different energy content.

3.4. Outlooks and clinical implications

The present study can be further enhanced by including geometric features of a realistic aortic root reproducing the characteristic curvature of the ascending aorta. The non-planar shape of the root axis could introduce additional complexity to the flow field, as the swirl illustrated by Morbiducci et al. (Reference Morbiducci, Ponzini, Rizzo, Cadioli, Esposito, Montevecchi and Redaelli2011). The aortic arch can also generate secondary vortical regions which can potentially affect both mean flow field and fluctuations (Zhu, Seo & Mittal Reference Zhu, Seo and Mittal2018). Secondly, the adoption of pressure boundary conditions at the cross-flow domain edges certainly provides a more physiological flow evolution, especially if coupled with a lumped circulation model (Lee et al. Reference Lee, Rygg, Kolahdouz, Rossi, Retta, Duraiswamy, Scotten, Craven and Griffith2020), and consequently, a more realistic prediction of the flow features during valve closing (Dasi et al. Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007).

The FSI dynamics of a realistic valve model should include also the arterial wall elastic effects. A computational investigation by Hsu et al. (Reference Hsu, Kamensky, Bazilevs, Sacks and Hughes2014) proved that the flexibility of the aortic root damps out the fluctuations in the flow rate and attenuates leaflet oscillations in biologic valves. By mitigating these further sources of unsteadiness, we speculate that weaker turbulent fluctuations can be observed. Flamini, DeAnda & Griffith (Reference Flamini, DeAnda and Griffith2016) found that a compliant aortic wall substantially generates smaller peak pressure at systole and lower retrograde flow in the diastolic stage than those attained by rigid models. In the case of a mechanical prosthetic model with prescribed flow rate, limited differences in second-order statistics are expected.

The LTFV and SJMVVG models showed reduced turbulent stress levels if compared with the SJMV model, therefore, they can potentially replace it. However, this trend must be confirmed by clinical experimentation, which will involve further practical issues to account for. Concerning the VGs, their presence might induce a local accumulation of pannus and other fibrovascular tissues (Kim et al. Reference Kim, Choi, Kweon, Yang and Kim2020), which in turn can disrupt the associated haemodynamics. For this reason, the design of VGs must be limited to the simplest geometric elements. One possible investigation would require the application of a proper orthogonal decomposition method (De Cillis et al. Reference De Cillis, Cherubini, Semeraro, Leonardi and De Palma2021) to understand the topology of the modes responsible for the stress reduction and optimise the VGs collocation on the leaflets. We argue that in our study the VGs lay close to the leaflet leading edge, thus, they are subjected to the largest flow velocity. In this configuration the flow might wash the VGs surface and prevent the formation of biologic deposits, but this can be confirmed only by in-vivo experimentation. The VGs can guarantee improved haemodynamics without loosing decades of clinical experience with the SJMV model, since no large-scale design modifications are needed on the valve components. On the other hand, the LTFV certainly provides haemodynamic features closer to those of the natural valve, along with a noiseless closure. Older versions of a trileaflet valve were affected by the risk of leaflet sticking (Ghista Reference Ghista1976) which heavily influenced their development. However, the LTFV configuration is declared to guarantee an unimpeded leaflet opening (Novostia 2019).

Further analyses should address whether the beneficial haemodynamic effects of the VGs and LTFV hold for an intra-annular implantation. The reduced effective orifice area is expected to affect the TPD, the vorticity dynamics downstream of the leaflets as well as the recirculation within the Valvaslva sinuses. However, these arguments remain controversial (Kim et al. Reference Kim, Kim, Kim, Jung, Choo, Chung and Lee2019), and the related speculations should be confirmed by further analyses, especially in the presence of VGs.