Turbulent Systolic Flow Downstream of a Bioprosthetic Aortic Valve: Velocity Spectra, Wall Shear Stresses, and Turbulent Dissipation Rates

Every year, a quarter million patients receive prosthetic heart valves in aortic valve replacement therapy. Prosthetic heart valves are known to lead to turbulent blood flow. This turbulent flow field may have adverse effects on blood itself, on the aortic wall and on the valve performance. A detailed characterization of the turbulent flow downstream of a valve could yield better understanding of its effect on shear-induced thrombocyte activation, unphysiological wall shear stresses and hemodynamic valve performance. Therefore, computational simulations of the flow past a bioprosthetic heart valve were performed. The computational results were validated against experimental measurements of the turbulent flow field with tomographic particle image velocimetry. The turbulent flow was analyzed for disturbance amplitudes, dissipation rates and shear stress distributions. It was found that approximately 26% of the hydrodynamic resistance of the valve was due to turbulent dissipation and that this dissipation mainly took place in a region about one valve diameter downstream of the valve orifice. Farther downstream, the turbulent fluctuations became weaker which was also reflected in the turbulent velocity spectra of the flow field. Viscous shear stresses were found to be in the range of the critical level for blood platelet activation. The turbulent flow led to elevated shear stress levels along the wall of the ascending aorta with strongly fluctuating and chaotic wall shear stress patterns. Further, we identified leaflet fluttering at 40 Hz which was connected to repeated shedding of vortex rings that appeared to feed the turbulent flow downstream of the valve.


INTRODUCTION
Aortic valve replacement (AVR) is a common therapy for moderate to severe aortic stenosis (Kheradvar et al., 2015). In AVR, the diseased native aortic valve is replaced by mechanical (MHV) or biological heart valve prostheses (BHV) in an estimated 250'000 annual interventions worldwide (Yoganathan et al., 2004). Heart valve prostheses are known to create turbulent flow in the aortic root (AR) and the ascending aorta (AAo) (Sotiropoulos et al., 2016). This turbulent flow may be connected to blood platelet activation triggering thrombus formation (Quinlan and Dooley, 2007) and to unphysiological wall shear stresses in the AAo possibly leading to endothelial lesions (Davies et al., 1986). Moreover, turbulent dissipation contributes to the pressure drop across the valve (clinically known as trans-valvular pressure gradient, TVPG). Flow instabilities at the level of the valve orifice could be related to leaflet fluttering in BHVs (Peacock, 1990) which may contribute to structural valve deterioration (SVD) in BHVs due to wear and fatigue of the leaflet tissue. Therefore, hydrodynamic instabilities and turbulent flow past BHVs are important phenomena which may be directly linked to clinically adverse events, and it may be advisable to design prostheses which lead to less turbulent flow.
Despite the relevance of turbulent flow past heart valve prostheses, our current understanding of this flow phenomenon remains incomplete. Early investigations of turbulence past heart valve prostheses include a study by Yoganathan et al. (1986) who based their findings on the turbulent shear stresses obtained from a laser Doppler anemometer for several MHV and BHV designs. They found that BHVs exhibit regions of high turbulent shear stresses at the edge of the central aortic jet issuing from the valve orifice during systole. Lim et al. (1998) studied three types of MHVs and one BHV made from porcine tissue with particle image velocimetry (PIV) measurements of steady systolic flow. Their investigation was mainly based on the analysis of Reynolds shear stresses. The authors detected the highest Reynolds shear stresses downstream of the studied BHV at constant flow rate. This result was attributed to the small orifice area of this valve which led to a strongly accelerated aortic jet with sharp shear layers.
The review on heart valve fluid mechanics by Yoganathan et al. (2004) emphasizes the role of the central aortic jet in BHVs and its regions of high shear. It provides a range of values for turbulent stresses that are thought to be relevant for blood platelet activation and thrombus formation. The review by Sotiropoulos et al. (2016) critically assesses the role of the Reynolds stresses as they do not contain information on instantaneous stresses acting on blood cells or platelets. They refer to the study by Ge et al. (2008) who investigated Reynolds and viscous stresses in the wake of a bileaflet MHV and who argued that Reynolds stresses should be used with caution when assessing blood cell damage since they only provide a statistical description on inertial phenomena and generally overestimate the stresses experienced by blood cells. They further emphasized that instantaneous viscous stresses are of greater interest for quantifying regions of excessive stresses on blood particles. The viscous stress levels in the study of Ge et al. (2008) were found to be large enough for platelet activation, while they did not have the magnitude for causing red blood cell damage or hemolysis. A similar finding was reported by Quinlan and Dooley (2007) who based their study on measurements of the flow field in the wake of three different bileaflet MHVs by Liu et al. (2000). Fluctuation values were obtained from phase-averaged quantities. It was shown that for different phases of the heart pulse, the inertial subrange for turbulent momentum transport agreed roughly with Kolmogorov's −5/3-power law (Pope, 2000, p. 230). Quinlan and Dooley (2007) further developed a mathematical model for the estimation of shear stress induced on blood cells based on a spherical cell model. They showed that in laminar flow, the induced blood cell stress is approximately equal to the bulk shear stress of the flow. In turbulent flow, they argue, that the stress experienced by blood particles is approximately one order of magnitude less than the Reynolds shear stresses. They concluded that the Reynolds stresses themselves cannot be used as a sole indicator for cell loading. Rather, the spectral energy distribution of the flow field should be studied and compared to the size of blood particles to identify relevant scales and their fluctuation magnitude for flow-induced blood trauma.
Among others, Morbiducci et al. (2009), Alemu et al. (2010, and Min Yun et al. (2014) used computational models of bileaflet MHV to study shear-induced platelet activation. The computational study by Hedayat et al. (2017) compared platelet activation in the turbulent flow fields past bileaflet MHV with a model of a BHV and found that the activation potential is higher for MHV. Hedayat and Borazjani (2019) showed that small vortical structures in systolic flow of a bileaflet MHV contribute significantly to their thrombogenic potential. Experimental investigations by Hasler et al. (2016) and Hasler and Obrist (2018) presented three-dimensional PIV data of the flow past a BHV. The measured fields showed turbulent flow with instantaneous shear strain rates beyond 2, 000 s −1 (Figure 1). Given that the limited resolution of these PIV measurements probably led to an underestimation of the strain rates, these results indicate that also the turbulent flow behind BHVs may lead to platelet activation.
Despite the intriguing results from Figure 1, the experimental data by Hasler and Obrist (2018) does not allow for a full characterization of the turbulent flow field in the systolic flow downstream of the BHV. The limited spatial resolution of the measurement effectively acts as a low-pass filter which prevents an accurate calculation of turbulent viscous dissipation and wall shear stresses, and the phase-averaged nature of the data renders a temporal analysis of flow phenomena difficult. Therefore, the present study aims to provide more detailed data on the turbulent systolic flow using a computational model. The results include the statistical characterization of the turbulent flow field during peak systole and a quantification of viscous shear stresses in the bulk flow and of wall shear stresses along the AAo wall. It is shown that approximately 26% of the transvalvular pressure gradient are due to turbulent dissipation. Viscous shear stresses in the turbulent flow can be as high as 40 Pa and turbulent wall shear stresses along the AAo wall are found reach values beyond 12 Pa. Velocity spectra, that were extracted at different locations in the turbulent flow field, exhibit an inertial subrange which adheres to Kolmogorov's −5/3-power law. Furthermore, a vortex shedding phenomenon leading to leaflet fluttering is described and quantified. It is not the aim of this study to directly assess the thrombogenic potential of a specific valve type or to predict other adverse events.
The results of the present study were obtained with a numerical solver for fluid-structure interaction (Nestola et al., 2019) which comprises a high-order flow solver for the simulation of transitional and turbulent flow. The computational model included a BHV and an anatomically correct model of an aortic root including three sinus portions (without coronary FIGURE 1 | Shear rates in the systolic flow behind a BHV. Adapted from experimental results by Hasler and Obrist (2018). ostia). The flow rate was smoothly ramped up until a quasi-steady systolic flow was obtained. This flow state was maintained for 0.1 s to allow for a temporal statistical assessment of the turbulent flow field.

Governing Equations
For this study we considered blood as a Newtonian fluid with density ρ f and dynamic viscosity µ f whose velocity field v f = (v f ,x , v f ,y , v f ,z ) and pressure p f was governed by the Navier-Stokes equations for incompressible flow, The solid structures (BHV, aortic root) satisfied the elastodynamics equations where u s is the displacement field, ρ s the density of the solid and P the Piola-Kirchhoff stress tensor. Boundary conditions were applied at selected boundary points to hold the structures in place within the fluid domain. The symbol · indicates values in the reference configuration.
The interaction between solid and fluid was described by interface conditions which were satisfied on the surface of the structures with unit normal vector n, where F = ∇ u s + I denotes the deformation gradient plus the identity matrix I, J = det F is its associated determinant and σ f is the Cauchy stress of the fluid. More details on the governing equations of the fluid-structure interaction problem are given in Nestola et al. (2019). Different constitutive laws were employed to simulate the structural response. The valve ring and aortic root were modeled as linearly elastic. Following the approach of Auricchio et al. (2014), the valve leaflets were modeled using the anisotropic fiber-reinforced Holzapfel-Gasser-Ogden material with two fiber families (Holzapfel et al., 2000). Its Piola-Kirchhoff stress tensor P( u s ) = ∂ /∂ F was derived from the scalar strain energy function with constitutive parameters µ s for the bulk material mechanics and k ij for the two fiber families j = 1, 2.Ī 1 ,Ī 4,1 , andĪ 4,2 denote the modified invariants I 4,1 = J −2/3 g 0,1 · C g 0,1 I 4,2 = J −2/3 g 0,2 · C g 0,2 The fiber orientations were defined by the unit vectors g 0,1 and g 0,2 , while C denoted the right Cauchy-Green strain tensor. Incompressibility was controlled by adding a volumetric energy term with the penalty coefficient κ to the strain energy function.

Model Configuration
We studied the systolic flow through a BHV of nominal size 21 mm which was inserted in a model of an aortic root (AR). The fluid density and viscosity were ρ f = 1, 060 kg/m 3 and µ f = 0.006 Pa s which is higher than the typical viscosity of blood (0.003 . . . 0.004 Pa s). This viscosity was chosen to reduce the computational cost. The effect of this choice on the results will be critically discussed below. Apart from an overestimation of the stress levels for blood, it will be shown that the overall structure of the flow field is not affected by this choice. The model was configured to yield an average TVPG (difference between aortic and left ventricular pressure) of 8 mmHg. The flow rate was approximately 14.5 l/ min which is a typical systolic flow rate for a mean cardiac output of 4 to 5 l/ min (Betts, 2013;Jahren et al., 2015).

Valve Leaflets
The geometry of the valve model was based on a commercial BHV (Edwards Intuity Elite 21 by Edwards Lifescience, Irvine, CA, USA) valve which was also used in the reference experiments (Figure 2A). Figure 2B shows the corresponding CAD drawing obtained from manually measuring the dimensions of the valve. The dimensions used for the CAD model are presented in Figure 2D. It was assumed that this configuration corresponds to the stress-free reference configuration of the leaflets (Kamensky et al., 2015). The material parameters µ s = 20.1 kPa, k 11 = k 12 = 54.62 kPa, and k 21 = k 22 = 30.86 for the Holzapfel-Gasser-Ogden material in Equation (2) were chosen as in Auricchio et al. (2014). Accordingly, the fiber orientation with respect to the horizontal component of the orthogonal basis was set to β = 30 • which results in a relative orientation between the two fiber families of 2β = 60 • (Figure 2C). Incompressibility of this part of the geometry was achieved through the penalization in Equation (6) with κ = 3 · 10 4 . The density of the material was set to a value of ρ s = 1, 100 kg/m 3 .

Valve Ring
The valve ring (dark gray structure in Figure 2B) comprised all structural elements of the BHV except for the valve leaflets. In the real BHV, the valve ring is made from a thin nitinol wireframe, which shapes the valve posts. The remaining voids in the geometry are then mostly composed from polymer material and covered with textile material. This complex composite structure was modeled as a cylindrical solid ring with three valve posts (dimensions were obtained manually with a caliper and are given in Figure 2D) using a single linearly elastic material. As the mechanical properties for the valve ring were not known, they were estimated and adjusted ad hoc to a shear-modulus of µ s = 0.6 MPa, a bulk modulus of κ s = 3 MPa and a density of ρ s = 1, 500 kg/m 3 .

Aortic Root (AR)
The geometry for the AR, including the sinus portions and part of the AAo, was a parametrized geometry which was assembled from data by Reul et al. (1990) and Haj-Ali et al. (2012). A silicone phantom based on the same parametrized model has been previously used by Hasler and Obrist (2018) in tomographic PIV measurements (geometry M in that study). Figure 3 shows the AR geometry in a cross-section through the xy-plane and a cross-section through the xz-plane together with the values of the individual radii and heights of the sections. The complete setup used for the simulations is shown in Figure 4, which also includes the dimensions of the computational domain for the fluid.
The coordinate system (x, y, z) originates at a point on the central axis of the AR at a height corresponding to the maximum sinus radius r s . The sinotubular junction (STJ) is an important landmark located at z stj = 0.0088 m and marks the point of transition from the bulbous sinus portion to the straight AAo.
The AR was modeled as a linearly elastic material with a shear-modulus µ s = 0.3 MPa, a bulk modulus of κ s = 3 MPa and a density of ρ s = 1, 100 kg/m 3 . We chose a linear elastic material description because (a) the resulting strains in the AR remained small and (b) the numerical results were compared to experimental results with an aortic root phantom manufactured from silicone (Hasler and Obrist, 2018). The chosen values of µ s and κ s therefore yield a material that corresponds to a volume ratio of 1 : 5 between curing agent and polymer according to Armani et al. (1999) if the Poisson ratio is chosen to be ν s = 0.45.

Numerical Method
The governing Equation (1) were solved with a dedicated solver for fluid-structure interaction problems Nestola et al. (2019) which is based on the immersed boundary method (Peskin, 2002). To this end, the valve and the AR were immersed in a computational fluid domain with dimensions 45 × 45 × 97.5 mm (Figure 4). The flow was discretized with a 6 th -order finite difference scheme on a rectilinear grid with 129×129×193 grid points. Grid stretching yielded increased resolution in the center of the fluid domain with a mesh width below 0.25 mm.
The structure was discretized with a finite-element method using approximately 98'000 P1 tetrahedral elements (73'163 for the AR; 18'204 for the valve ring; 6'523 for the valve leaflets). Data was transferred between fluid and structural meshes using a variational formulation of transfer operators for arbitrary meshes (Krause and Zulian, 2016). The reader is referred to Nestola et al. (2019) for details on the numerical implementation, convergence and code validation.
The simulations were run with a time step size of t = 5 · 10 −6 s on 256 Intel Xeon E5-2690 v3 CPU cores (Piz Daint, Cray XC40/50 supercomputer at the Swiss National Supercomputing Center CSCS) resulting in a wall time on the order of 5 days for the simulation of 0.3 s of physical time.
Instead of directly imposing pressure or velocity boundary conditions, we used periodic boundary conditions for the fluid domain. The flow through the AR was driven by a forcing term f TVPG which was added to the right-hand side of (1a), where λ(x) is equal to unity within a cylindrical domain of height 6 mm (magenta region in Figure 4 extending from z = −0.0277 m to z = −0.0217 m with radius r a = 0.011 m) and zero otherwise. The z-component of this forcing term yielded a local pressure elevation of 8 mmHg which led to a flow through the valve with a mean TVPG of 8 mmHg. The x-and y-components of the forcing term (7) have the character of a fringe forcing which attenuated the inflow by penalizing lateral velocities. In summary, the forcing f TVPG fulfilled two roles: (a) creating a flow through the valve with a mean TVPG of 8 mmHg, (b) eliminating velocity fluctuations at the inflow (which may otherwise re-enter the AR from the outflow and through the periodic boundary conditions).  However, we will see in the following, that the fringe forcing does not suppress all incoming velocity fluctuations such that residual fluctuations remain present in the inflow to the BHV.

Analysis of the Flow Field
The resulting flow field was analyzed for various instantaneous and statistical physical quantities. Table 1 gives an overview on the basic physical quantities used to characterize the flow through the BHV. The computational model was configured to maintain a statistically steady flow through the valve model after an initial transient. This quasi-steady state was attained at t ≈ 0.2 s. This

Quantity Symbol Formula
Aortic flow rate For the present configuration, these values are taken at z 0 =0.01 m which is slightly above the STJ. The integrals were taken over the AR lumen and H(·) is the Heaviside step function which is zero for negative arguments and unity for positive arguments.
allowed us to compute the mean of a quantity q(t) as For the present configuration, we used t 1 = 0.2 s and t 2 = 0.3 s. For the analysis of the turbulent flow downstream of the BHV, the Reynolds decomposition of the flow field v f yielded the mean It will be shown in the following that the flow field comprised a periodic component related to vortex shedding at 40 Hz. To separate the coherent structures of this periodic component from the turbulent fluctuations, we performed a triple decomposition of the flow field according to Reynolds and Hussain (1972) which further decomposed the Reynolds fluctuation v ′ f into a periodic componentṽ f and a turbulent component The periodic componentṽ f was computed as which includes a phase averaging the flow field over N periods of length T. The magnitude of the fluctuations was characterized by the root-mean-square (mrs) of the fields By normalizing the rms values with the mean flow magnitude V f = |V f |, we obtained turbulence intensities I and I ′′ and the intensityĨ of the periodic fluctuation, The shear strain rate and shear stress in the flow field were quantified via the strain rate tensor S, The maximum instantaneous viscous shear stress was then computed as where s max and s min are the maximum and minimum eigenvalues of S. The Reynolds stress tensor was defined as from which we derived a maximum Reynolds shear stress as where σ max and σ min are the maximum and minimum eigenvalues of { ij }.
The turbulent dissipation rate ǫ was computed according to The wall shear stress τ w acting on the aortic wall was computed as where n was the unit normal vector of the wall pointing into the fluid domain. The magnitude of the wall shear stress was then given as from which we computed the time averaged wall shear stress as The oscillatory and/or turbulent character of the wall shear stress is commonly quantified by the oscillatory shear index which was computed as A value of OSI = 0 indicates unidirectional wall shear stress and 0.5 indicates oscillatory wall shear stress (Lee et al., 2009). Furthermore, we computed the time averaged magnitude of the WSS fluctuations as OSI and TAWSSF were used to characterize turbulent character of the wall shear stress. Figure 5 shows the evolution of the flow field from t = 0 to 0.3 s along with snapshots of vortical structures in the flow field at three distinct stages (see also the supplementary flow field video for an animated sequence of the flow field). At t = 0.046 s, immediately after valve opening, a starting vortex (V0 in Figure 5) was shed from the leaflet tips. This structured vortex ring was accompanied by a first peak of the average jet velocity v jet,avg . A similar three-lobed starting vortex was described in Sotiropoulos et al. (2016). Continued acceleration of the flow (increasing Q and v jet,avg , cf. Table 1) led to more complex vortex patterns. The second snapshot at t = 0.125 s highlights this development. At this point in time, the starting vortex V0 had already been advected out of the computational domain and three additional vortex rings (V1, V2, V3) had been shed from the valve orifice at regular intervals with a frequency close to 40 Hz (associated with peaks of v jet,avg ). While the most recent vortex ring (V3 in Figure 5) was still an intact closed ring, the vortex ring V2 had been advected approximately one valve diameter downstream and was in the process of breaking down. Vortex V1 had already broken down completely to small scale vortical structures. After t = 0.2 s, the quantities in Figure 5, which were taken close to the valve orifice at z = 0.01 m, continued to show oscillations suggesting that there was continued shedding of vortex rings at 40 Hz. The flow field (third snapshot in Figure 5 taken at t = 0.24 s) exhibited a fully developed turbulent flow indicating rapid laminar-turbulent transition of the aortic jet and the vortex rings. This led to a Reynolds number for the flow through the AAo of Re AAo = 2ρ f Q/(πµ f r AAo ) ≈ 2000. The Reynolds number for the aortic jet was estimated as Re jet = ρ f v jet,avg r jet /µ f ≈ 2100.

Leaflet Fluttering
The vortex shedding at 40 Hz was directly connected to repeated transversal deflections of the valve leaflets which is also known as leaflet fluttering (Peacock, 1990). This fluttering showed amplitudes of several millimeters and had the same frequency f flutter ≈ 40 Hz as the vortex shedding. Figure 6 depicts the leaflet motion throughout one fluttering period T flutter = 1/f flutter ≈ 0.025 s in a cross-section along its symmetry plane. The motion pattern had the form of a wave traveling from the base to the tip of the leaflet. The wave peak traveled at a mean velocity of approximately 0.4 m/s. The wave speed increased toward the end of the flutter period leading to a whiplash motion of the leaflet tip.

Mean Flow Field
The fully developed turbulent aortic jet in the center of the AR showed mean velocities beyond 2 m/s. The axial component V f ,z of the mean flow is depicted in Figure 7 for various cross-sections cutting the domain centrally in the xz-plane and yz-plane and at different xy-planes. Besides the central aortic jet, there were retrograde flow regions between the AAo wall and the aortic jet.
These regions started about one valve diameter downstream of  The cross-sections in Figure 7C show that the aortic jet was not axisymmetric. At z = 0.01 m, it featured a star-shaped region with higher velocities. The three points of the star were aligned with the valve posts and commissures. At z = 0.03 m, the cross-section resembled a six-pointed star which is the reason why the attachment points of the jet in th xz-plane ( Figure 7A) were farther downstream than the attachment points in the yz-plane (Figure 7B). At z = 0.05 m, these features had mostly disappeared and the character of the mean flow field had transitioned from a free jet to a pipe flow.

Turbulent Flow Field
The breakdown of the aortic jet into small eddies is illustrated in Figure 8 and in the supplementary flow field video. The vortex rings, which were issued at 40 Hz, broke down rapidly into turbulent structures. Visual inspection of Figure 8 suggests that the flow field was most chaotic approximately one valve diameter downstream of the orifice. Farther downstream, the turbulent flow field appeared to calm down. These qualitative observations are in line with experimental observations shown in Figure 1 and we will demonstrate in the following that they can be substantiated quantitatively.
Figures 9A-C shows the root-mean-square (rms) fluctuation field v ′ rms in the same cross-sectional planes as in Figure 7. The highest rms values were found above the valve posts which coincided with the orientation of the tree-pointed star within the central jet. Therefore, the highest fluctuation amplitudes were found in regions where the mean flow featured the sharpest shear layers. The center of the aortic jet showed very little velocity fluctuations as it is common for potential cores of jets. With increasing downstream distance, the fluctuation field became more homogeneous and had a reduced amplitude suggesting decaying turbulent intensity and increasingly homogeneous turbulent flow.
The triple decomposition of the flow field, Equation (10), was computed to better discriminate between the coherent vortex structures shed from the BHV at a frequency of 40 Hz and the turbulent fluctuations. To this end, the flow field between t = 0.2 s and 0.3 s was phase averaged over four fluttering periods of T flutter = 0.025 s according to Equation (11 Figure 9F shows the turbulence intensities I and I ′′ and the intensityĨ of the periodic fluctuations along the centerline and along an axial line (−0.007 m, 0 m, z) which passed through the shear layer. This data confirms the observations from above and shows that the intensityĨ of the periodic fluctuations had a distinct peak at z ≈ 0.015 m. This peak also marks the point where the turbulence intensity I ′′ in the shear layer began to increase. Further, we find that the turbulence intensity at the inflow (after the fringe region) was approximately 10% which indicates that the fringe region did not suppress all fluctuations. However, we also observe that the turbulence intensity I ′′ on the centerline decreased to 5% within the valve suggesting that the residual turbulent fluctuations were attenuated in the accelerated flow through the valve. The turbulence intensity I ′′ increased again beyond 10% only after the peak ofĨ at z ≈ 0.015 m.

Viscous Shear Stress and Reynolds Shear Stress
Maximum viscous shear stresses τ max (15) were computed for the flow field at t = 0.3 s (Figures 10A-C). The resulting stress field was consistent with the above observations in that it showed peak values at a height of z ≈ 0.03 m. The colorbar is set to an upper limit of 16 Pa but peak values within the whole domain reached up to 40 Pa. Figures 10D-F shows the maximum Reynolds shear stress RSS max (17). The regions of high RSS max at the location FIGURE 8 | Vortical structures at t = 0.24 s visualized by isosurfaces of the λ 2 -criterion by Jeong and Hussain (1995).
of the leaflets were artifacts due to the immersed boundary implementation of the fluid-structure interaction and should not be considered for the analysis of the turbulent flow field. Otherwise, the Reynolds stress field showed a similar structure as the rms velocity fluctuations (Figure 9) and exhibited values up to 300 Pa within the shear layers and above the valve posts.

Turbulent Dissipation Rate
The turbulent dissipation rate ǫ (18) is illustrated in Figure 11. Consistent with the distribution of the rms fluctuations (Figure 9), the highest dissipation rates were in the shear layers above the valve posts. Beyond z = 0.03 m, the dissipation rates decreased in amplitude. There was an almost dissipation-free region in the jet core which closed at z ≈ 0.05 m.
Similar to the plots for RSS max , Figures 11A,B include artifacts due to the immersed leaflet structures. Although the dissipation rates at these locations are physically not meaningful, they visualize the extent of the leaflet excursions during multiple fluttering periods. Figure 12 shows the evolution of the energy dissipation rate ǫ along the centerline of the AR (dotted red curve) and along a parallel line (dotted blue curve) which passed through the shear layer at the edge of the aortic jet at {x 0 , y 0 , z} = {−0.007 m, 0 m, z}. Whereas, the dissipation rate along the centerline was very small at the valve orifice and increased in the downstream direction, the dissipation rate along the shear layer had a peak between z = 0.02 and 0.03 m after which it decayed toward the value of the centerline dissipation rate.
The average dissipation rate ǫ avg (average over the crosssection at the given height) is shown in Figure 12 by solid blue and red lines (the two curves are identical). Evidently, the average dissipation rate reached a distinct peak at z ≈ 0.027 m after which it decreased monotonically. This suggests that the flow experienced the most intense turbulent dissipation approximately one valve diameter past the valve orifice which is in line with our previous observations.
The total turbulent dissipation rate in the AAo was where the integral was evaluated over the AAo volume from z = 0.01 m to z = 0.05 m. This turbulent dissipation rate can be compared to the rate of mechanical energy loss of the flow through the BHV which amounted to P mech = p mean · Q mean = 0.245 W, where p mean is the mean TVPG. Comparison of the two values shows that 26% of the total mechanical losses can be attributed to turbulence. Or in other words, approximately 2.5 mmHg of the TVPG of 8 mmHg were due to turbulent flow.

Velocity Spectra
Spectra of the velocity fluctuations were studied at different locations. To this end, the power spectral density of time series of velocity fluctuations were computed with MATLAB's pwelch function. Figure 12 shows the power spectral density of the axial velocity fluctuation v ′ f ,z sampled at different heights z on the centerline and on a parallel line passing through the shear layer of the jet.
In both data sets, the leaflet fluttering frequency f flutter ≈ 40 Hz dominated the spectra close to the valve orifice (z = 0.01 . . . 0.02 m). These peaks can be associated with the vortex shedding from the valve orifice. Farther downstream, the peaks vanished while broad-banded spectra were established. This process was completed at z ≈ 0.03 m where the average turbulent dissipation rate reached its peak value. Figure 13 shows two velocity spectra E zz of the turbulent flow field. They were obtained from autocorrelation functions R zz (x ′ ) for two points on the centerline at z = 0.03 and z = 0.05 m. Both spectra comprised a wavenumber range with −5/3-slope which is reminiscent of a turbulent inertial subrange (Pope, 2000, p. 230). At z = 0.05 m, the energy in this range was somewhat lower than at z = 0.03 m indicating a decay in turbulent intensity in downstream direction. Figure 14A illustrates the WSS distribution for the fully developed flow at t = 0.24 s. We found irregular WSS patterns mainly downstream of the STJ with peak values beyond The mean WSS ( Figure 14B) exhibited values mostly below 3 Pa. This value is roughly four times higher than the WSS of a We conclude from these observations that the turbulent systolic flow led to high WSS levels. Highest WSS fluctuations were noted beyond z = 0.03 m where the free aortic jet attached to the AAo wall (cf. Figure 7). Finally, it should be pointed out that (somewhat weaker) turbulent WSS patterns were also identified between z = 0.01 and 0.03 m where the flow field close to the wall was dominated by retrograde flow.

DISCUSSION
It was the primary aim of this study to characterize the turbulent systolic flow of a BHV. To this end, we designed a computational model of a BHV in an anatomically correct model of an AR. The governing equations were solved with a FSI solver (Nestola et al., 2019) which comprises a finite-element structural solver for soft tissue and a high-order Navier-Stokes solver that has been designed for the study of laminar-turbulent transition (e.g., Obrist et al., 2012;John et al., 2016). The FSI solver was verified and validated for canonical benchmarks in Nestola et al. (2019). For validation of the present model, we will show in the following that the computational model yielded results which were consistent with experimental observations. In particular, we will compare our results to experimental results by Hasler and Obrist (2018) and Vennemann et al. (2018) who studied the same bovine BHV as in the present work.
The topology of the mean flow with a central aortic jet with maximum velocities of approximately 2.5 m/s and retrograde flow at approximately 0.5 m/s along the AAo wall (Figure 7) agrees very well with the experimental observations of Hasler and Obrist (2018). These observations are also consistent with the instantaneous flow field at peak systole reported by Lee et al. (2020) which was obtained with a computational model for a similar BHV configuration and which exhibited the same flow topology with a non-axisymmetric central jet and retrograde flow along the AAo wall.
For further validation of the computational results against experimental data, selected results from the study of Hasler and Obrist (2018) are shown in the Supplementary Material (PIV data) where the experimental data were processed and plotted in the same way as the results of the present study. Inspection of this data shows that there is also good agreement between experimental and computational data for the turbulent flow fields: The rms fluctuation velocities (Figure 9 and   jet at z ≈ 0.03 m. Likewise, RSS max (Figures 10D-F  and Supplementary Figure 2) peak in the same region at approximately 300 Pa. The concentration of rms fluctuations and RSS above the valve posts could not be observed in the experiment. We suspect that this difference was due to different inflow conditions in experiment and simulation. It is possible that higher inflow disturbance levels in the experiment led to a less organized flow field which made it more difficult to identify local structures above the valve posts. Finally, Supplementary Figure 3 shows velocity spectra E zz (k x ) for the same locations as the spectra in Figure 13. Although these spectra exhibit the same overall structure as the computational data, the experimental velocity spectrum at z = 0.03 m (Supplementary Figure 3A) shows a significantly lower energy level than the computational data. This may be attributed to the laminar core which had not yet fully closed at z = 0.03 m in the experiment. At z = 0.05 m, the laminar core had closed in the experiment and in the computational model and the turbulent flow field attained a more homogeneous character. Accordingly, the experiential velocity spectrum at z = 0.05 m (Supplementary Figure 3B) is in very good agreement with the computational data.
The leaflet fluttering at 40 Hz agrees very well with observations of Vennemann et al. (2018) who described leaflet fluttering at 36 Hz in an experiment with similar hemodynamic configuration and the same valve design. The present results are also close to the experimental observations of Peacock (1990) who reported fluttering at 15 to 30 Hz for a different bovine BHV in water. Recent in vivo electrocardiographic measurements in stentless bovine BHV found fluttering frequencies of 15 Hz (Aljalloud et al., 2018). Lee et al. (2020) reported fluttering of a bovine BHV in a pulse duplicator in saline at approximately 30 Hz while their numerical model exhibited fluttering at approximately 60 Hz (estimated from Figure 4 in Lee et al., 2020). These comparisons suggest that the present model was appropriate to reproduce the phenomenon of leaflet fluttering, although the fluttering frequency was somewhat higher than in corresponding experimental settings. Quantitative differences to the experimental results of Moore and Dasi (2014), who observed fluttering at 50 − 100 Hz for a porcine BHV in saline and no fluttering when the a water/glycerin mixture was used, may be due to the porcine tissue with results in thinner and more supple leaflet structures than in bovine BHV. The amplitude of the fluttering appears rather high when compared to experimental data. Whereas, the leaflet tips in the present study moved approximately 5 mm during a fluttering period, Peacock (1990) reports amplitudes of only 2 mm and Moore and Dasi (2014) observed fluttering amplitudes of 4 mm. It is not clear at this point, what causes these large amplitudes and further studies are necessary to understand which mechanical or geometrical parameters determine the fluttering amplitude.
We showed that leaflet fluttering was directly connected to repeated shedding of vortex rings at 40Hz. These vortex rings are an inherent feature of the fully developed aortic jet and must not be confused with the starting vortex which is only issued once after valve opening (Sotiropoulos et al., 2016). The triple decomposition of the flow field (Figures 9D-F) suggests that the fluctuations due to the vortex rings had a peak at z ≈ 0.015 m after which the turbulence intensity I ′′ started to increase, marking the start of the vortex ring breakdown. This process was also reflected in the temporal spectra of Figure 12 which were dominated close to the valve orifice (at z = 0.01 m) by the shedding frequency of 40 Hz and then evolved to turbulent spectra featuring the classical −5/3-slope of the inertial subrange (Figure 13). The peak in turbulent dissipation at z ≈ 0.03 m suggested that this region may have been the place of turbulent breakdown. This region was also associated with the highest viscous shear stresses and Reynolds shear stresses (Figure 10).
While RSS levels are of limited relevance for the prediction of blood trauma (Quinlan and Dooley, 2007;Ge et al., 2008), the viscous shear stress levels (up to 40 Pa) indicate that shearinduced thrombocyte activation may occur in the shear layers of the aortic jet. As already pointed out by Hedayat et al. (2017) and Hedayat and Borazjani (2019) for the case of bileaflet MHV, this activation in the aortic jet may become clinically relevant, if activated thrombocytes are transported by the retrograde flow along the wall toward the sinus portions where thrombi may form (Chakravarty et al., 2017;Jahren et al., 2018;Hatoum et al., 2019).
The turbulent flow behind the BHV left its footprint also on the AAo wall. We found elevated and fluctuating WSS up to 14 Pa mainly downstream of the point where the aortic jet attached to the wall (Figure 14). These turbulent WSS may play a role in lesions of the endothelial layer (Davies et al., 1986). The present results also indicate that the distribution and intensity of turbulent WSS on the AAo wall are related to the point of attachment of the central aortic jet. This suggests that the ratio between aortic jet diameter and AAo diameter could be relevant for the location and magnitude of peak WSS: the larger the AAo diameter, the later the attachment and the lower and farther downstream the peak WSS. Of course, these conjectures must be confirmed quantitatively and it should also be pointed out that the present model with a straight AAo did not reflect the full geometrical complexity with the bending of the AAo which may lead to an early impingement (rather than an attachment) of the aortic jet.
Further limitations of the present study include the inflow, the modeling of the BHV and the high fluid viscosity used in our model. As we have seen in Figure 9F, the inflow to the BHV had a residual turbulence intensity I of approximately 10% due to incomplete suppression of fluctuations in the fringe region. It is unclear whether this level of fluctuations is a good representation of the inflow to a BHV coming from the left ventricle or the inflow in an experimental setup. It can be expected that the complex and unsteady flow within the left ventricle will lead to significant inflow disturbances. However, their magnitude is not known and further studies are necessary to assess to what extent the inflow perturbations affect the transition process behind the valve.
The BHV was modeled using manual measurements of the valve dimensions and literature data on the mechanical properties of bovine pericard tissue used for the valve leaflets. Although the valve kinematics matches experimental observations generally quite well, the high fluttering amplitude suggests that the present BHV model should be further refined by optimizing mechanical and geometrical parameters of the model.
To assess the effect of the fluid viscosity on the results, we performed ad hoc simulations with a simplified model with viscosities of 0.004, 0.006, and 0.008 Pa s. These tests indicated that the viscosity had a negligible effect on the mean flow field and only a very small effect on the location of the turbulent breakdown. For the lower viscosity (0.004 Pa s), we found an increase of 10% of the turbulent kinetic energy with respect to the nominal viscosity (0.006 Pa s) and an decrease of 5% for the highest viscosity (0.008 Pa s). The average viscous stresses changed nearly proportionally to the viscosity which also indicates that the shear rates remained nearly unchanged. Therefore, it should be expected that the viscous shear stress levels and the wall shear stress levels for blood are reduced by 30 to 40%.
Finally, we would like to discuss differences between the pulsatile flow through a real heart valve and the quasi-steady systolic flow configuration in the present study. Certainly, the present configuration eliminated some transient effects due to flow pulsatility and we found that fluctuations of r jet , v jet,avg , and v jet,max slightly decayed toward the end of our simulation (cf. Figure 5). This could indicate that the fluttering is not sustained. However, careful inspection of the leaflet kinematics showed that one of the leaflets started to flutter slightly out-of-phase from the other leaflets. Therefore, the reduced fluctuations toward the end of the simulation were due to asynchronous leaflet motion and not due to reduced individual leaflet fluttering amplitudes. Further, we believe that there is sufficient time scale separation between central flow phenomena (e.g., fluttering at 40 Hz and the typical duration of the systolic phase approximately 0.3 s) and we found that the flow had enough time to establish the turbulent flow field during the early systolic phase. Therefore, we believe that the studied flow field was representative for the flow during peak systole. Additionally, the quasi-steady configuration had the advantage that turbulence statistics could use time-averaged quantities and that there were no artifacts due to cycle-to-cycle variations which are known to occur in pulsatile flow past heart valves (Sotiropoulos et al., 2016).

CONCLUSION
The present computational study characterized the turbulent systolic flow downstream of a bovine BHV. Similar to the study by Lee et al. (2020), we validated our numerical model against experimental data obtained for the same valve design.
We identified a peak of turbulent dissipation approximately one diameter downstream of the valve orifice. This was also the region where the vortex rings broke down, that were shed at 40 Hz from the fluttering leaflet tips. The total turbulent dissipation was found to be accountable for 26% of the total pressure loss over the valve (TVPG) indicating that turbulence is a significant and detrimental factor for hemodynamic valve performance. Further, we found elevated turbulent viscous shear stresses up to 14 Pa which may be connected to shearinduced thrombocyte activation. This could indicate that BHV thrombogenicity (Chakravarty et al., 2017) is connected to the turbulent flow past the valve.
To our best knowledge, this computational study is the first to present turbulent spectra behind BHV. As pointed out by Quinlan and Dooley (2007), whose study was limited to data for mechanical heart valves by Yoganathan et al. (1986) and Liu et al. (2000), such velocity spectra are an important basis to study shear-induced thrombocyte activation. Finally, the present computational model allowed us to study for the first time the wall shear stress patterns along the AAo wall (Figure 14 and WSS video in the Supplementary Material). We found elevated levels of turbulent WSS which suggest a possible connection between BHV turbulence, endothelial lesions and long-term AAo health.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
BB contributed to the development of computational model, running the computational model, post-processing, analysis and interpretation of data, writing and critical revision of the manuscript, and study design. LP contributed to the postprocessing and analysis of experimental data and critical revision of the manuscript. DO contributed to the post-processing, analysis and interpretation of data, writing and critical revision of the manuscript, and study design and supervision. All authors contributed to the article and approved the submitted version.