1. Introduction
A mechanical heart valve is a durable prosthetic device used to replace damaged heart valves, ensuring proper blood flow and preventing complications such as heart failure or stroke. Commonly made from materials like pyrolytic carbon, it is desirable for the valve to last for a lifetime, even for younger patients, with an aim to improve quality of life by alleviating symptoms such as fatigue and breathlessness. The St Jude Medical (SJM) valve is among the more prominent ones used for aortic or mitral valve replacement, featuring two semicircular leaflets that open and close to regulate the blood flow optimally. This bileaflet design aims to achieve the desired blood flow control and minimise turbulence, as compared to older mechanical valve designs. Made from durable, biocompatible pyrolytic carbon, this valve appears to be ideally suited for long-term functionality, but requires lifelong anticoagulation therapy to prevent blood clots (Langenaeken et al. Reference Langenaeken, Vanoppen, Janssens, Tanghe, Verbrugghe, Rega and Meuris2023). As a precautionary measure, lifelong anticoagulation therapy is often recommended, increasing bleeding risks and requiring regular international normalised ratio monitoring. It is thus unsuitable for those with poor medication adherence or high bleeding risks. From haemodynamic considerations, turbulent vortices around such mechanical heart valves can cause adverse effects such as haemolysis – red blood cell (RBC) damage – increased risk of thrombosis and platelet activation, and potential valve dysfunction due to wear, compelling design modifications to mitigate these effects (De Tullio et al. Reference de Tullio, Cristallo, Balaras and Verzicco2009; Yun et al. Reference Yun, Aidun and Yoganathan2014b ; Hedayat & Borazjani Reference Hedayat and Borazjani2019; Nitti, De Cillis & de Tullio Reference Nitti, De Cillis and de Tullio2022).
Undesirable turbulence structures may also lead to energy loss, increased workload of the heart, and endothelial damage, impairing the valve efficiency. Additionally, turbulence can elevate the risk of infective endocarditis by promoting endothelial trauma. The underpinnings of platelet activation or RBC damage are attributed to the levels of high mechanical load and exposure times to the blood elements (Bluestein et al. Reference Bluestein, Rambod and Gharib1999, Reference Bluestein, Li and Krukenkamp2002). Therefore, large-scale flow disturbances and local recirculation zones in combination may lead to higher blood damage. Reynolds stress due to turbulent fluctuations may further augment the mechanical stress on the blood cells (Antiga & Steinman Reference Antiga and Steinman2009) if the length scale of the fluctuations (Kolmogorov scale) is closer to the blood cell dimension. These adverse effects emphasise the importance of improved valve designs to minimise the potential complications in patients undergoing prosthetic valve replacement. Refining the leading and trailing edges of a mechanical heart valve could potentially address the above challenges by promoting smoother blood flow, reducing turbulence and wall shear stress, and minimising the risks of haemolysis and thrombosis. However, this approach remains unexplored in the paradigm of mechanical heart valve design. To overcome this deficit, we present here direct numerical simulations (DNS) of a newly designed bileaflet mechanical heart valve (BMHV) with streamlined leading and trailing edges. This design modification aims to reduce endothelial damage and patient discomfort caused by turbulence or audible noise, offering safer and more efficient valve performance by better mimicking physiological blood flow patterns. Our DNS results demonstrate the effectiveness of streamlined edges in minimising abrupt blood flow direction changes, preventing flow separation, and reducing the formation of turbulence-inducing eddies and vortices. Compared to the state-of-the-art SJM valve, the new design facilitates smoother flow transitions across valve surfaces, and minimises flow disturbances. The findings highlight the impact of edge streamlining in reducing obstruction, drag and wake turbulence downstream, ensuring smoother blood flow and balanced pressure gradients. This design significantly reduces shear stress, enhances haemodynamic efficiency, and lowers the risk of complications such as haemolysis and thrombosis, offering a promising advancement in mechanical heart valve technology.

Figure 1. (
$a$
) Representation of blood flow through the left ventricle and aorta – the native aortic valves lie between the left ventricular outflow tract and the aortic root, which consists of three bulges corresponding to the sinuses of Valsalva. (
$b$
) Schematic of the computational domain models the ventricular and aortic chambers as straight cylindrical tubes separated by the valve ring consisting of the valve and the physiological sinuses of Valsalva (Haj-Ali et al. Reference Haj-Ali, Marom, Ben Zekry, Rosenfeld and Raanani2012). (
$c$
) Side view of the domain. The origin is fixed at the leftmost section of the valve ring in the front view (
$y$
and
$z$
directions coincide with the centre of the tube). (
$d$
) Physiological flow rate–time relationship at the inlet of the ventricle (marked by the orange line), where selected phases of interest are highlighted with filled green circles: mid acceleration (MA),
$t=0.1$
s; late acceleration (LA),
$t=0.16$
s; peak flow (PF),
$t=0.21$
s; early deceleration (ED),
$t=0.26$
s; and late deceleration (LD),
$t=0.32$
s. (
$e$
) The base case valve of the study (a 23 mm SJM Regent valve model). (
$f$
) The proposed valve model (STE).
2. Mathematical model
2.1. Geometric configuration
The native aortic valve, positioned between the left ventricular outflow tract and the aorta, is represented in figure 1(
$a$
). We have modelled this configuration as blood flow through a straight, cylindrical chamber of the ventricle and the aorta (figure 1
$b$
). This model also includes the physiological shape of the sinuses of Valsalva corresponding to the healthy aortic valve model given by Haj-Ali et al. (Reference Haj-Ali, Marom, Ben Zekry, Rosenfeld and Raanani2012) (figure 1
$c$
). In this aortic domain, we have considered a 23 mm SJM Regent BMHV model (figure 1
b,e). The SJM valve, positioned between the ventricular and aortic chambers, is free to pivot around its hinges due to the action of flow. The ventricular and aortic chambers have diameter
$d$
and lengths
$2d$
and
$3d$
, respectively. Other researchers have also considered similar geometrical configurations (Yun et al. Reference Yun, Wu, Simon, Arjunon, Sotiropoulos, Aidun and Yoganathan2012; Hedayat, Asgharzadeh & Borazjani Reference Hedayat, Asgharzadeh and Borazjani2017; Zolfaghari et al. Reference Zolfaghari, Kerswell, Obrist and Schmid2022).
The SJM valve exhibits a low angle of incidence (approximately
$5^{\circ}$
as measured from the horizontal plane) when fully open, and closes at a high angle of incidence of approximately
$60^{\circ}$
in each cardiac cycle (Dasi et al. Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007; Borazjani, Ge & Sotiropoulos Reference Borazjani, Ge and Sotiropoulos2008; Yun et al. Reference Yun, Dasi, Aidun and Yoganathan2014a
). In this work, we further propose a valve model obtained by streamlining the leading and trailing edges of the original SJM model (figure 1
$f$
), which is abbreviated as the STE valve model. We imposed identical maximum and minimum leaflet angular positions, the same as the SJM valve, and conducted two-way coupled fluid–structure interaction (FSI) simulations for all cases. Note that we have kept a gap of approximately 0.5 mm between the leaflet’s side faces and the valve ring for computational stability.
2.2. The DNS
The blood flow through the aortic root and valve model is described by the unsteady, three-dimensional, incompressible Navier–Stokes equation, which in dimensional form reads

where,
$\boldsymbol{u}=(u, v, w)$
represents the velocity components in the
$(x, y, z)$
directions, respectively,
$p$
is the pressure (
$P$
) divided by the density (
$\rho$
), and
$\nu ={\unicode{x03BC}} /\rho$
denotes the kinematic viscosity of the fluid.
The motion of the leaflets is governed by Newton’s law of conservation of angular momentum around the hinge axis:

where
$\theta =\theta (t)$
is the leaflet’s angular position,
$I$
is the moment of inertia of the leaflets about their pivots, and
$M$
is the moment exerted by the flow on the leaflets. The moment
$M$
is determined by integrating the moments due to pressure and viscous forces over the respective leaflet surfaces about their hinge location. The value of
$I$
is distinct for different valve models.
2.3. Boundary conditions
In this study, we have considered a physiological temporal inflow profile with period
$T$
, corresponding to a heart rate of
$70 \,\textrm{beats min}^{-1}$
. This profile resembles the one considered in the in vitro pulsatile flow loop for BMHVs by Dasi et al. (Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007) and in the numerical simulations conducted by Borazjani et al. (Reference Borazjani, Ge and Sotiropoulos2008). The volumetric flow rate prescribed at the inlet of the aortic model is depicted in figure 1(
$d$
). We prescribe a sectionally uniform and temporally varying velocity profile at the inflow boundary based on the physiological inflow rate as described in § 2.1. The choice of a uniform profile is made due to the high Womersley number characterising this flow (Zamir Reference Zamir2000).
The aortic walls are considered to be rigid in this work. At the confining walls, we enforce no-slip boundary conditions, while at the outlet of the flow domain, we impose traction-free outflow boundary conditions. The kinematic constraints of the leaflets are imposed based on the design specification of maximum and minimum angle as mentioned in § 2.5 below.
2.4. Details on solution methodology
2.4.1. Fluid flow solver
The Navier–Stokes equation (2.1) is solved numerically using a fractional-step method (Hirt & Cook Reference Hirt and Cook1972) on a staggered Cartesian grid for primitive variables. A second-order upwind scheme is used for the convective terms, while a second-order accurate central difference is used for the diffusive terms. The time integration is done by using the forward Euler scheme. The method results in a Poisson equation for pressure correction, which is the most time-consuming step in the entire algorithm. The overall scheme is explicit, rendering the flow solver to be very fast in predicting solution of time-dependent problems having small length and time scales, such as the current problem.
The basic flow solver is augmented with a sharp interface-based immersed boundary (IB) method to handle non-grid conforming boundaries (Raj et al. Reference Raj, Khan, Alam, Prakash and Roy2023). It is an indirect variant of the discrete forcing IB method. In this method, the solid body is discretised as a set of Lagrangian markers (triangular meshes) that are tracked throughout the domain, segregating the entire domain into fluid, solid and intercepted cells. The Navier–Stokes equation is solved in the fluid cells, while a special reconstruction scheme is used in the intercepted cells, which serves as boundary conditions to the fluid cells. The IB method is particularly suited for moving boundary problems (such as the present case) compared to body-fitted methods since the Eulerian fluid mesh need not be morphed at every time step. Readers are referred to Raj et al. (Reference Raj, Khan, Alam, Prakash and Roy2023) for the IB method in general, and Sarkar et al. (Reference Sarkar, Dawn, Raj, Khan and Roy2024a ) for boundary condition imposition and force calculation in the IB method. The overall flow solver has a second-order spatial accuracy (Raj et al. Reference Raj, Khan, Alam, Prakash and Roy2023) and can handle complex geometries such as aerofoils (Raj & Roy Reference Raj and Roy2023), stenosed straight arteries (Kumar, Roy & Das Reference Kumar, Roy and Das2024), stenosed curved arteries such as S-bend (Khan et al. Reference Khan, Raj, Alam, Chakraborty and Roy2023a ), and stenosed branched arteries such as carotid (Khan et al. Reference Khan, Sharma, Chakraborty and Roy2023b ).
2.4.2. Structural solver and FSI scheme
The structural solver involves the integration of (2.2) for rigid leaflets, and maintaining no-slip and no residual stress conditions at the interface separating the fluid and solid domains. The leaflet surface is discretised using triangular mesh elements with a density such that each fluid cell consists of at least three triangular mesh elements. The forward Euler scheme is used to advance the structural equation.
The FSI scheme (Sarkar et al. Reference Sarkar, Dawn, Raj, Khan and Roy2024a ) is based on a Dirichlet Neumann coupling of the fluid and solid solvers. It calls the fluid and solid solvers iteratively based on a loose coupling approach. First, the fluid solver is advanced in time using the previously known position and velocity of the structure. Then the forces (or moments) applied on the structure are calculated using the new field variables and old position of the structure. Finally, the structural system (2.2) advances to obtain its new position using a suitable time marching scheme. This method of coupling is explicit and is very fast compared to the strong coupling method, and is also found to be computationally stable for simulations involving low-mass-ratio solids (Sarkar et al. Reference Sarkar, Dawn, Raj, Khan and Roy2024a ). The FSI solver has a second-order spatial accuracy, and has been verified for multiple FSI problems involving both rigid and flexible solid bodies: vortex-induced vibrations of an elastically mounted circular cylinder, sedimentation of a circular disk, vortex-pair interaction with thin leaflets, pulsatile flow past a BMHV in an axisymmetric aorta, and flow-induced deformation of a thin elastic plate attached to a circular cylinder (Sarkar et al. Reference Sarkar, Dawn, Raj, Khan and Roy2024a ). The scheme has also been demonstrated in the study of BMHV flows due to Newtonian and non-Newtonian rheology models (Sarkar et al. Reference Sarkar, Sharma, Chakraborty and Roy2024b ).
2.5. Geometrical parameters
The diameter of the aortic and ventricular chamber is
$d=25.4$
mm. Parameters describing the casing, which consists of the ventricle, aorta, valve ring, hinge position and sinuses of Valsalva, are reported in table 1. Also reported are the parameters of the SJM valve and the proposed STE valve.
Table 1. Parameters of the case study.

2.6. Flow parameters and fluid properties
The period of the pulsatile cycle is
$T=0.86$
s. The inflow features a peak flow rate
$Q_{\textit{max}}= 24.2$
$\textrm {min}^{-1}$
and mean flow rate
$Q_{\textit{mean}}= 4.9$
$\textrm {min}^{-1}$
. To ensure optimal valve closure, a maximum reverse flow rate
$Q_{\textit{peak leakage}}=0.2 Q_{\textit{max}}$
is prescribed, corresponding to the adverse pressure gradient generated by ventricular expansion during diastole. The forward flow duration lasts approximately
$T_{\textit{systole}}=0.34$
s. These parameters are also shown in table 1.
The density of blood is assumed to be constant at
$\rho =1000$
kg
$\textrm {m}^{-3}$
, and the kinematic viscosity is
$\nu =3.5 \times 10^{-6} {\textrm{m}}^2 {\textrm{s}}^{-1}$
. The assumption of Newtonian behaviour of blood is supported by the notably high shear rates observed in large arteries such as the aorta (Pedley Reference Pedley1980), along with the insignificant deviations from the linearity of the shear rate–shear stress relationship (Chhabra & Richardson Reference Chhabra and Richardson2008). Despite the evidence, recent numerical studies (De Vita, de Tullio & Verzicco Reference De Vita, de Tullio and Verzicco2016; Sarkar et al. Reference Sarkar, Sharma, Chakraborty and Roy2024b
) have highlighted the impact of blood’s non-Newtonian rheology on estimating output parameters such as leaflet kinematics, flow features, blood damage, and so on, in BMHV flows, in comparison to Newtonian counterparts. However, it is crucial to acknowledge the uncertainty surrounding the input parameters and the structure of the non-Newtonian model, which may influence these resultant output parameters (De Vita et al. Reference De Vita, de Tullio and Verzicco2016; Sarkar et al. Reference Sarkar, Sharma, Chakraborty and Roy2024b
). Additionally, the viscoelastic nature of blood may also contribute to affecting these output parameters, warranting further investigation. Hence in this study, we choose the simplified Newtonian model to facilitate a clear comparison of the performance of different valve types. This approach enables us to elucidate the differences in output parameters resulting from changes in valve designs without the added complexities associated with non-Newtonian models.
2.7. Dimensionless numbers
The Reynolds number (
${\textit{Re}}$
), representing the ratio of inertial to viscous forces, is defined by the diameter of the aorta (
$d$
) and a chosen velocity scale (
$u$
) as
${\textit{Re}}=u d/\nu$
. Various velocity scales may be considered, such as the cycle-averaged bulk inflow velocity and peak bulk inflow velocity, resulting in
${\textit{Re}}_{\textit{mean}}=1170$
and
${\textit{Re}}_{\textit{peak}}=5780$
, respectively.
Another dimensionless number characterising the time-varying inertial forces to viscous force, known as the Womersley number (
$\textit{Wo}$
), can be defined as
$Wo = ({d}/{2})\sqrt {{2 \unicode{x03C0} }/{(\nu T)}}$
. In our study, a high
$Wo=18.34$
is used, which is consistent with the characteristics of large arteries such as the aorta.
2.8. Particulars on simulation parameters – validation and verification
In this subsection, we show the validation of our numerical scheme by comparing velocity profiles downstream of the leaflets with particle image velocimetry (PIV) and lattice Boltzmann method (LBM) results of Yun et al. (Reference Yun, Dasi, Aidun and Yoganathan2014a
) at a constant inflow rate corresponding to different flow regimes persisting in the actual pulsatile cycle. The valves are fixed at the fully open position for the velocity profile comparisons. We employ uniform mesh size 107
${\unicode{x03BC}} \textrm {m}$
in the valves and sinus regions, and a suitable grid stretching is used elsewhere. A time averaging of 400 time steps is used to obtain the mean velocity profiles. Figure 2 shows the axial velocity profile validation at
${\textit{Re}}=5000$
, representing the near peak flow rate phase. Here,
$x^*=0$
represents the downstream tip of the leaflet. Our numerical model predicts the velocity exceptionally well with the experimental PIV data in the central jet region at larger downstream distances. The validations for
${\textit{Re}}=2400$
and
$750$
are also found to match well with the literature (see Sarkar et al. Reference Sarkar, Sharma, Chakraborty and Roy2024b
).

Figure 2. Comparison of mean axial velocity profiles for steady inflow rate at
${\textit{Re}}=5000$
:
$(a)$
$x^*=3.1$
mm,
$(b)$
$x^*=7.5$
mm,
$(c)$
$x^*=13.8$
mm,
$(d)$
$x^*=18.3$
mm, with PIV and LBM results of Yun et al. (Reference Yun, Dasi, Aidun and Yoganathan2014a
). Figure reproduced from Sarkar et al. (Reference Sarkar, Sharma, Chakraborty and Roy2024b
) with permission.
For the grid independence test, a total of three meshes were considered, i.e.
$631 \times 240 \times 240$
(coarse),
$703 \times 301 \times 301$
(medium) and
$879 \times 400 \times 400$
(finest), giving rise to uniform cell sizes
$r_d/160=133.75$
${\unicode{x03BC}} \textrm {m}$
,
$r_d/200=107$
${\unicode{x03BC}} \textrm {m}$
,
$r_d/267=80.15$
${\unicode{x03BC}} \textrm {m}$
, respectively, in the valves and sinus region. The ratios of cell volume for consecutive grids are
$V_{{coarse}}/V_{\textit{medium}}=1.95$
and
$V_{\textit{medium}}/V_{\textit{fine}}=2.38$
.

Figure 3. Comparison of (a–c) mean axial velocity and (d–f) r.m.s axial velocity at steady peak inflow (
${\textit{Re}}=5780$
): (a,d)
$x^*=5.0$
mm, (b,e)
$x^*=10.1$
mm, (c, f)
$x^*=20.1$
mm. Figure reproduced from Sarkar et al. (Reference Sarkar, Sharma, Chakraborty and Roy2024b
) with permission.
Figure 3 shows comparison of the mean and root mean square (r.m.s.) axial velocity values at different axial distances in the perpendicular plane bisecting the valves (
$z=0$
) at steady inflow velocity corresponding to the peak inflow of the pulsation cycle (
${\textit{Re}}=5780$
). A total of 50 000 time-point averaging is utilised to obtain statistical results in this highly unsteady and disorganised flow. The medium grid predictions are reasonably in accordance with those of the finest grid, especially in the top half of the core flow and recirculation region (above
$y=0$
mm). At a further downstream distance
$x^*=20.1$
mm (figure 3
c), the medium grid predictions reasonably match those of the finest grid in the outer jet region. The r.m.s profiles of axial velocity of the medium grid also show good convergence with those of the fine grid. Considering the computational costs, the medium grid with grid size 107
${\unicode{x03BC}} \textrm {m}$
is used in the current study. In this study, we have considered a constant time step
$10^{-5}$
s satisfying the Courant–Friedrichs–Lewy (CFL) condition (
${\lt } 0.1$
) for numerical stability. It will be shown, in subsection 3.5.1, that the grid size employed is of the same order as the Kolmogorov length scale, and hence can capture the lower-order moments of turbulent flow present in BMHVs. We use a constant time step size
$10^{-5}$
s based on the CFL criteria at peak flow rate. The time step employed was found to be one order smaller than the Kolmogorov time scale, and can capture the turnover time of small eddies in the turbulent flow due to BMHVs (shown later).
Finally, we report the validation of the leaflet kinematics in figure 4. Our loosely coupled FSI predictions of the leaflet kinematics are closer to the experimental PIV results of Dasi et al. (Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007). Our model also matches well with the strongly coupled FSI simulation of Borazjani et al. (Reference Borazjani, Ge and Sotiropoulos2008). Our loose coupling scheme gives stable prediction in both opening and closing phases of the leaflets.

Figure 4. Comparison of leaflet angle with experimental PIV observation of Dasi et al. (Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007) and hinge resolved one-way and two-way coupings of Hedayat & Borazjani (Reference Hedayat and Borazjani2019).
We have also compared our leaflet kinematics with Hedayat & Borazjani (Reference Hedayat and Borazjani2019), who have resolved the hinge completely by using an overset mesh, and employed both one-way and two-way interpolation with the base fluid mesh (see figure 4). The opening valve kinematics is unaffected by the hinge model, but the closing kinematics has maximum difference 15–20 ms, depending on the hinge directional interpolation.
2.9. Phases of interest and collection of statistics
In this study, we focus our interest on a few select phases of the cardiac cycle, elucidating the prominent features of the flow field due to different BMHV designs. These phases, which correspond to distinct inflow rates (as governed by the time-dependent inflow rate), are: mid acceleration (MA) at
$t_1=0.1$
s, late acceleration (LA) at
$t_2=0.16 {\textrm{s}}$
, peak flow (PF) at
$t_3=0.21 {\textrm{s}}$
, early deceleration (ED) at
$t_4=0.26 {\textrm{s}}$
, and late deceleration (LD) at
$t_5=0.32\, {\textrm{s}}$
, keeping in view the period of cardiac cycle (
$T=0.86\, {\textrm{s}}$
) such that time instants
$t_i$
and
$T+t_i$
constitute the same phase for all
$i$
from
$1$
to
$5$
. The corresponding instantaneous flow rates and Reynolds numbers are
$Q/Q_{\textit{max}} = 0.53, 0.86, 1, 0.85, 0.31$
and
${\textit{Re}}_{\textit{inst}} = 3074, 4979, 5780, 4887, 1799$
, respectively. These phases are marked on top of the inflow rate profile (figure 1
$d$
).
In this study, both the instantaneous and phase-averaged flow fields are presented. Phase-averaged statistical data are collected over 22 cardiac cycles (
$N=22$
). Due to the requirement of an immensely large number of cycles to collect data for statistically significant values, we create a time window
$T_i$
centred around the phases of interest
$t_i$
having half-width
$T_i/2$
such that data may be collected for all the samples residing within the time window over multiple cycles (
$N=22$
) for the current phase of interest. Hence each phase-averaged quantity is computed by an ensemble average over
$N$
time-averaged fields. Similar techniques have been used by other investigators in experimental (Liu, Lu & Chu Reference Liu, Lu and Chu1999; Haya & Tavoularis Reference Haya and Tavoularis2016) and computational (Yun et al. Reference Yun, Dasi, Aidun and Yoganathan2014c
; Nitti et al. Reference Nitti, De Cillis and de Tullio2022) studies related to BMHV flows. The time window
$T_i$
is small enough that the flow field is statistically steady, but large enough to collect statistically meaningful results. In that regard, the time window is calculated such that the volume flow rate change is negligible (
${\approx}1\,\%$
) over the time window. The phase MA is neglected for phase averaging due to the laminar, structured flow nature persisting in this phase.
The mean field for the
$x$
component of velocity
$\overline {u}$
, obtained by averaging the instantaneous
$u$
fields within the time window
$T_i$
centred around the phase
$t_i$
, for
$N$
cycles, is

Accordingly, the instantaneous fluctuation of
$u$
is

and any second-order statistical quantity is collected as

Statistics for other components of velocity fluctuations are calculated likewise. The phase-averaged data are predominantly used in § 3.5.
2.10. Frequency spectra calculation
The discrete Fourier transform (DFT) of the
$y$
velocity squared signal was obtained at a downstream location of the trailing edge of the leaflet. The signal was sampled at frequency 10 000 Hz in a cardiac cycle and spanning seven such cycles. The DFT converts the discrete data points in the time domain to the equivalent frequency domain. It decomposes the signal into its constituent frequencies, amplitude and phase. Consider a sequence of
$N$
equally spaced complex numbers
$\boldsymbol{x}_n= x_0, x_1,\ldots ,x_{N-1}$
. Its DFT results in a sequence of complex numbers
$\boldsymbol{X}_k= X_0, X_1,\ldots ,X_{N-1}$
via the transformation

which can be expressed as
$X_k = A_k + {\textrm i}B_k$
, where
$A_k$
and
$B_k$
are the real and imaginary parts of
$X_k$
. The magnitude and phase of the
$k$
th frequency bin are obtained from
$M_k=\sqrt {A_k^2 + B_k^2}$
and
$P_k= \tan^{-1} ({B_k}/{A_k})$
. The magnitude corresponds to the amplitude of sinusoids at the
$k$
th frequency, and the phase denotes the phase shift of the sinusoids. The amplitudes are scaled with a factor
$2/N$
, and the frequencies above the Nyquist limit (which is half the sampling frequency) are neglected. The amplitudes corresponding to the different frequency bins are denoted by
$E_{v^2}$
(mm2
$\textrm {s}^{-2}$
) and are plotted on the
$y$
axis, and the frequency bins
$k$
values are plotted in Hz (
$=k \times \text{sampling frequency} / N$
).
3. Results and discussion
In this section, we first present the leaflet kinematics and the performance metrics for the proposed streamlined leaflet heart valve, and compare the results with a standard SJM model. Then we try to explain the behaviour of the flow structures, their fluctuations, and the differences arising due to the geometric modification of the leaflet shape. We further try to explain how the leaflet instabilities and the resulting turbulent fluctuations are alleviated for this streamlined valve model. We have also tested the performance of two other leaflet geometries, respectively with streamlined leading and trailing edges, and the observations are reported in Appendix A.

Figure 5. (
$a$
) Leaflet kinematics of SJM and STE valves for three cardiac cycles (first two cycles are neglected). Letters L and U stand for the lower and upper leaflets, respectively. Cycle to cycle kinematics of the lower leaflet ((
$b$
) opening phases, (
$c$
) closing phases) and upper leaflet ((
$d$
) opening phases, (
$e$
) closing phases). Asynchronous motion between upper and lower leaflets ((
$f$
) opening phases, (
$g$
) closing phases).
3.1. Leaflet kinematics
Figure 5(
$a$
) reports the kinematics of the base (SJM) and proposed (STE) valve models through three consecutive cardiac cycles. The letters L and U represent the lower and upper leaflets, respectively. We can see a nearly repeating behaviour of the leaflet motion over these three cycles. During the leaflet opening phase, there is little variation in the leaflet angle between the two valve models. This similarity persists for most of the opening phase, from the closed position of
$58^{\circ}$
–
$20^{\circ}$
(measured from the horizontal plane). Here, both valve models show very similar behaviour owing to accelerating inflow over them. We further zoom in on the closing and opening phases of different leaflets over these three cycles to understand the differences in their kinematics. Figures 5(
$b$
) and 5(
$c$
) show the opening and closing behaviours of the lower leaflets of the valve models, respectively. It can be seen that the proposed valve design, STE, has more delayed opening and closing compared to the base SJM valve model, and this behaviour is consistent over different cycles. This consistent delayed response in valve opening for the STE model can be attributed to reduced flow obstruction owing to its streamlined shape (smooth variation of leaflet area near the leading and trailing edges). On the other hand, the slower closing of the STE leaflet can be related to the fact that this valve shows less disorganisation in the flow field due to the absence of local recirculation zones and fine-scale vortex structures compared to the base model. We will discuss the flow structures in detail in later subsections. It is also noted that the base (SJM) valve leaflet motion shows overshoots and fluctuations near the end of the opening phase. These fluctuations are not seen in the proposed valve model, STE. Thus the reduced local fluctuations help in smoother opening and closing in the proposed valve model. Similar observations are made for the upper leaflets in the opening and closing phases (figure 5
d,e).
However, we may notice high levels of asynchronous motion between both the leaflets during the closing phase for both valve models (figure 5
$g$
). Flow fields for the different valve models show higher fluctuations due to the large number of small-scale structures during the fully open phase and during the leakage phase (shown in a later subsection). These fluctuations break the symmetry of the valve motion during closure. Hence high levels of asynchronous motion persist between the leaflets in the closing phase for both valve models. A lower degree of asynchronous motion is observed between the bottom and top leaflets during the opening phase (figure 5
$f$
), which could be attributed to the regular, laminar, periodically repeatable, accelerating flow encountered in the opening phase, which corresponds to the phase MA. A stable nature of flow persists in this opening phase (until
${\textit{Re}}_{\textit{inst}}\approx 2700$
).

Figure 6. Averaged leaflet kinematics of SJM and STE valves for (
$a$
) lower and (
$b$
) upper leaflets, respectively.
The averaged leaflet kinematics of base and proposed valve models are shown in figure 6. It can be noticed that the opening and closing for the proposed valve are smoother than those of the base model. Such smooth, gradual motion may be favourable in reducing the stress levels experienced by blood elements closer to the leaflet surface in the opening and closing phases for the proposed valve model.
Table 2 reports the respective time durations for the bottom and top leaflets’ opening, closing and fully open phases for the two valve models considered in this study. These values are averaged over five cardiac cycles. It is observed that the leaflet closure is more rapid than the opening for both the valve models, although leaflet opening and closing are more gradual for the proposed valve compared to the base model. This is on a par with the instantaneous leaflet kinematics discussed above. Table 2 also reports the maximum and average velocities of the leaflet’s tip (averaged over five cycles) for the bottom and top leaflets, respectively, for all the valve models. The average tip velocity of the proposed STE valve is smaller than that of the base SJM model in the opening phase, suggesting a gradual opening of the proposed valve, and a larger time duration of valve opening, compared to the base model. While qualitatively our prediction matches well with Vennemann et al. (Reference Vennemann, Rösgen, Heinisch and Obrist2018) in terms of the leaflet kinematics, quantitatively the values, for example, of maximum leaflet velocity during opening and closing, are different from those reported in Vennemann et al. (Reference Vennemann, Rösgen, Heinisch and Obrist2018) since a different mechanical heart valve was used in their study (Medtronic Advantage) along with a different position of valve placement.
Table 2. Time duration of leaflet opening, closing and fully open, as well as the maximum and average velocities of the leaflet’s tip of top and bottom leaflets, respectively, for the base (SJM) and proposed (STE) valve models, averaged over five cycles.

Note that we have modelled the hinge design by a simple pin joint in the
$z$
direction, as mentioned in table 1, which differs from the actual mechanism consisting of the leaflet ear fitting into the hinge recess housing, also known as a butterfly recess (Simon et al. Reference Simon, Ge, Sotiropoulos and Yoganathan2010; Yun et al. Reference Yun, Wu, Simon, Arjunon, Sotiropoulos, Aidun and Yoganathan2012). This simplification is made due to computational constraint. Even though it is expected that the dynamics of leaflet opening and closing change based on the choice of hinge mechanism, our validation of leaflet angular position with experimental results of Dasi et al. (Reference Dasi, Ge, Simon, Sotiropoulos and Yoganathan2007) (with the resolved hinge of the SJM valve) compares well (shown in § 2.8). Hence we have used the simplified pin joint mechanism in this study.
3.2. Effective orifice area, energy loss and pressure drop
To compare the performances of these two valve models, we have reported the effective orifice area (EOA) and transvalvular energy loss for these valves. The EOA specifies the effectiveness of valve design in utilising the internal orifice area of the valve ring over the entire cardiac cycle (Yoganathan et al. Reference Yoganathan, Chaux, Gray, Woo, DeRobertis, Williams and Matloff1984). The larger value of EOA implies reduced obstruction to the blood flow. We have used the method devised by Clavel et al. (Reference Clavel2011) for computing the EOA based on the continuity equation

where
$d_{\textit{LVOT}}$
is the left ventricular output tract diameter, and
$\textit{VTI}_{\textit{LVOT}}$
and
$\textit{VTI}_{Ao}$
stand for the velocity time integral
$\textit{VTI}=\int _0^{T} u\, {\textrm d}t$
at the LVOT and downstream of the valves, at
$x=-2$
mm and
$x=16$
mm, respectively.
Table 3 reports the EOAs for both the SJM and STE valves. We obtain EOA 3.4
$\textrm {cm}^{2}$
for the SJM valve. Similar values are reported in the literature: e.g. the EOA is 3.18
$\textrm {cm}^{2}$
for a 23 mm SJM valve used in Gray et al. (Reference Gray, Chaux, Matloff, DeRobertis, Raymond, Stewart and Yoganathan1984) . From the table, we can observe that the proposed valve delivers a slightly higher EOA (4.4 %) compared to the base model.
Table 3. Parameters EOA and energy loss for SJM and STE valves.

Another important parameter for quantifying the valve performance is the transvalvular energy loss. Higher energy loss can lead to an increased workload on the left ventricle. Hence lower energy loss is preferred. Using the procedure mentioned in Azadani et al. (Reference Azadani, Jaussaud, Matthews, Ge, Guy, Chuter and Tseng2009), we have calculated the total energy loss during a cycle using

where
$\varPhi _{\textit{TEL}}$
is total energy loss,
$t_o$
and
$t_e$
are the start and end times of a cycle,
$Q_{valve}$
is the instantaneous flow through the valve, and
$\Delta P$
is the instantaneous pressure gradient. We can observe in table 3 that the total energy loss has been reduced substantially (
$\sim\!30$
%) for the proposed valve compared to the base valve model. This reduction in energy loss is due to the flow being much more stable for the STE valve compared to the SJM valve. Furthermore, the losses due to turbulent fluctuations also demand attention, which we will discuss later, in § 3.5.
Figure 7 shows the axial distribution of phase-averaged pressure
$\overline {P}-\overline {P}_0$
, where
$\overline {P}_0$
is the phase-averaged pressure at the inlet at the peak flow rate phase. There is a substantial drop in the pressure in the valvular and aortic root regions following a pressure recovery in the aorta. It can be observed that there is a decrease in the pressure drop for the proposed STE valve as compared to the SJM valve. The streamlined leaflet shape has resulted in a reduction in turbulence-inducing vortices (shown later), which lead to lower pressure drop and energy losses. The transvalvular pressure drop (TPD) is 4.97 mm Hg for the SJM valve, which is close to the value reported by Nitti et al. (Reference Nitti, De Cillis and de Tullio2022) (5.04 mm Hg). The TPD in the case of STE is 3.35 mm Hg.

Figure 7. Axial pressure distribution scaled by pressure at the inlet,
$\overline {P}(x)-\overline {P}_0$
, at the phase-averaged peak flow rate phase, for both SJM and STE valves.
3.3. Three-dimensional features
Figure 8 shows the three-dimensional vortex structures marked by the Q-criterion for select instances during the fully open phase of the valves for both models. Only the structures due to the upper leaflet are reported for brevity in the visualisation.

Figure 8. The Q-criterion iso-surface in the core flow for SJM and STE valve models at various time instances (
$Q=40\,000$
).
During mid-acceleration (
$t=0.11$
s), structures such as vortex tubes or rolls form downstream of the trailing edge of the leaflet of the SJM model due to Kelvin–Helmholtz (KH) type instabilities in the free shear layer between the upper and central jets. Further (
$t=0.13$
s), the axis of the vortex rolls turns and does not remain parallel to the
$z$
axis. At
$t=0.14$
s, the vortex roll makes an S-bend-like shape in the middle region, which then lifts upwards in the
$y$
direction. At
$t=0.15$
s, this central and upward-lifted portion of the vortex roll splits into two, making two upward-facing bumps. This phenomenon continues until
$t=0.15{-}0.18$
s, after which hairpin-like vortex structures are found (
$t=0.19{-} 0.26$
s). Then the late deceleration phases show a chaotic small-scale nature.
However, this complex reformation to hairpin vortices is suppressed for the STE model. Such vortex structures are not observed in the immediate downstream vicinity of the leaflet, suggesting a delayed formation of KH instability due to modification of the cross-sectional profile of the leaflets of the STE valve model.
Topological complexities of the vortex tubes may be quantified by the helicity of the flow field. It is also reported that helicity and fluctuating kinetic energy are related in prosthetic heart valve flows (Gallo et al. Reference Gallo, Morbiducci and de Tullio2022). Quantification of helicity or amount of swirl in flow may indicate phenomena such as mixing or turbulence. It can be used to correlate the amount of blood damage due to unphysiological flow. Hence a quantitative measure of absolute helicity is defined as
$h^*=\boldsymbol{u \boldsymbol{\cdot }} ( \boldsymbol{\nabla } \times \boldsymbol{u} )$
. A normalised helicity is also defined as
$h=\boldsymbol{u} \boldsymbol{\cdot } (\boldsymbol{\nabla } \times \boldsymbol{u} ) / ( | \boldsymbol{u} |\, | \boldsymbol{\nabla } \times \boldsymbol{u} | )$
, where
$h \in [-1, 1]$
indicates the relative orientation of the velocity vector to the vorticity vector. Figure 9 shows the normalised helicity contours for the SJM and STE models, respectively. We can observe negligible normalised helicity values in the central region past the leaflets for the STE valves compared to the SJM model (figure 9
$a$
), as it has a more coherent flow in this region. For the axial cross-sectional planes, in the sinus region (
$x=20$
mm, figure 9
$b$
), the STE model shows considerably lower helicity due to the leaflets compared to the SJM model. In the convergent section also, the STE model suggests lower helicity (
$x=30$
mm, figure 9
$c$
).

Figure 9. Contours of normalised helicity
$h$
at peak flow (PF) phase: (
$a$
)
$z=0$
mm plane, (
$b$
)
$x=20$
mm, (
$c$
)
$x=30$
mm, for SJM and STE models.
Quantitatively, the spatially averaged value of modulus of helicity (within a cuboidal region in the wake of the leaflets) for the STE model is
$39.76\,\%$
lower compared to the SJM model (table 4). Table 4 also reports spatially averaged helicity at other phases, and the STE valve shows considerably lower helicity values than the SJM valve.
To elucidate the flow features in more detail, a two-dimensional
$(x, y)$
cross-sectional plane is extracted, which bisects the valves at
$z=0$
and is illustrated in the following subsections.
3.4. Two-dimensional flow features
3.4.1. Vorticity patterns
Figure 10 shows the instantaneous
$z$
vorticity contours at a plane bisecting the valves corresponding to different phases of interest. We can notice two important vortex structures, namely the vortex sheet rolling up in the sinus due to the expansion after the valve ring, and the counter-rotating Burgers vortices formed due to the instability of the shear layers past the trailing edge of leaflets. A representative Burgers vortex is shown by a circle at phase LA for the SJM valve; see figure 10(
$b$
). At phase LA (figure 10
b), both these structures propagate further downstream, filling up the entire aortic region. The vortex sheet in the sinus sheds periodically due to shear layer instability, and this behaviour is similar in both valve models. In the core flow region, it is seen that the shedding due to the SJM valve is unorganised, with interactions between the Burgers vortices, rendering the flow to be non-symmetric about the centreline (
$y=0$
mm) at the downstream region. Notice that for the STE valve, the formation of Burgers vortices is delayed to a greater downstream distance. It displays superior flow characteristics by exhibiting a stable, straight shear layer past the leaflets until a greater downstream distance in the sinus, and forming the Burgers vortices at its departure, rendering the flow essentially symmetric about the centreline at this phase.
At the PF phase (figure 10 c), the SJM valve model exhibits much disorganisation, with the appearance of smaller-scale structures in the core flow region, while the stable, straight shear layer past the leaflets grows axially for the STE valve model, which is followed by the presence of coherent Burgers vortices (similar to the phase LA). Smaller-scale structures appear at the outer annular flow for all the valve models in this phase. Further, at the ED phase (figure 10 d), due to a reduction in the incoming flow rate, the eddies receive reduced momentum from the mean flow, and rapidly break down to even-smaller-scale structures for the SJM valve model. However, the STE valve shows a straight shear layer with delayed breakdown at this phase. Finally, at the LD phase, the cascading occurs further, with the small-scale structures diffusing in the flow due to viscosity in the core, annular and sinus regions. Notice that even at this phase, the shear layers due to the STE leaflets are much more stable, with a slight separation in the distance between the shear layer structures in the cross-stream direction between the positive and negative shear layers departing from each leaflet.
Table 4. Average modulus of helicity (in
${\unicode{x03BC}} \textrm {m}$
$\textrm {s}^{-2}$
) in the wake of the leaflets for SJM and STE valve models at various phases.

In summary, BMHV flows can be better stabilised by introducing subtle changes in the cross-sectional profile of the valve leaflets, even when other geometrical parameters and flow conditions are kept the same. The proposed STE valves exhibit stable flow structures over the major duration of the cardiac cycle, while the SJM valve shows drastic changes in flow structures within the cardiac cycle. In that regard, the STE valve model may induce favourable haemodynamics in BMHV flows, thereby attempting to reduce the potential thromboembolic complications due to mechanical heart valves.
Cycle-to-cycle repeatability, and hence reliability, of the flow features reported above is ensured by phase-averaged field data (not shown here for the sake of brevity) obtained from six consecutive cardiac cycles, which resembles the instantaneous fields in terms of large-scale structures. In the next subsubsection, we further present the streamtraces at select phases to explain the formation of different flow structures pertaining to the respective valve models.

Figure 10. Out-of-plane vorticity
$\omega _z$
at selected phases of the cardiac cycle for SJM and STE valves.
3.4.2. Streamtraces
Figures 11 and 13 show the streamtraces projected at the
$z=0$
plane for SJM and STE valve models. We focus near the upper leaflet of the SJM model at various time instants in figure 11. The behaviour over the bottom leaflet is largely similar. At the early portions of the cardiac cycle (when the valve becomes fully open,
$t=0.093$
s), a recirculation region is formed at the leeward side of the leaflet (encircled in green). Then it detaches from the solid surface and forms a recirculation vortex, which grows in size and eventually dissipates into the flow (
$t=0.099,0.102$
s). Later, another recirculation region forms at the lower part of the leeward face (
$t=0.106$
s). Similar to the previous one, it grows, separates, and dissipates into the flow. This phenomenon continues until
$t=0.1211$
s, at which stage both the vortices are seen together at the leeward side of the leaflet. Another kind of recirculation region forms at the pressure side of leaflets (corresponding to the inner surface). The genesis of the recirculation bubble is seen at the instant
$t=0.14$
s (marked by a green circle in the same figure). Afterwards, it grows further in intensity and size. Finally, all three of these recirculating regions (recirculation bubble and leaflet trailing edge vortices) appear concurrently before the peak of the pulsation cycle, and continue to exist until the early deceleration phase (
$t=0.1762{-}0.1896$
s; only the representative plot at
$t=0.1762$
s is shown). Figure 12 shows the streamtrace pattern around the entire leaflet. We can see a stretched vortex near the leading edge of the leaflet, and the separation bubble as discussed before. A similar leading edge vortex is reported by Zolfaghari & Obrist (Reference Zolfaghari and Obrist2021), which is attributed to the adverse pressure gradient in the inner surfaces of the leaflets (Zolfaghari & Obrist Reference Zolfaghari and Obrist2019).
The STE valve model attenuates the instabilities associated with the leaflets, as shown above, and no recirculation zones are found on the surfaces of the leaflets with no subsequent detachment and dissipation. Hence the streamtraces are straight, with no curvature, yielding overall flow stability. Such smooth variation in the flow area leads to favourable pressure variation near the valves. The initial oscillation in the streamtraces is due to the impact of the leaflet from opening to the fully open position, which occurs at approximately
$t=0.1$
s.
We have simulated two other leaflet designs, namely sharp leading edge and tapered trailing edge, which help in enunciating subsequent improvements in the flow pattern from SJM to STE; see Appendix A.1.

Figure 11. Streamtrace patterns of the SJM valve model.

Figure 12. Streamtrace patterns of the SJM valve model near the leading edge.
3.4.3. Frequency spectra
Figure 14 shows the frequency spectra of the
$y$
velocity squared at a point downstream of the leaflets (
$x=18$
mm,
$y=2.7$
mm) for both SJM and STE valve models. One can observe the obtained spectra to follow the
$-5/3$
slope of the inertial range of the energy cascade. The procedure for calculating the frequency spectra is provided in § 2.10.
Table 5. The most energy-containing peaks for the valve models considered in the present study.


Figure 13. Streamtrace patterns of the STE valve model.

Figure 14. Frequency spectra of the square of the vertical component of velocity for the (a) SJM and (b) STE valve models. The fundamental frequencies for the SJM valve are denoted by circles.
The fundamental frequencies from the charts (denoted by circles in the figure) are tabulated in table 5. The highest energy-containing frequency for the SJM valve is 96 Hz, with an amplitude of
$O(10^4)$
$\textrm {mm}^{2}$
s–2. This corresponds to the frequency of shedding of a separation bubble formed on the inner surfaces of the leaflet of the SJM model (as depicted in figure 12). This shedding frequency is verified by performing the DFT again in much smaller time window spanning only a few vortex shedding cycles near the peak flow rate phase. A similar shedding frequency is also reported in the literature (Bellofiore, Donohue & Quinlan Reference Bellofiore, Donohue and Quinlan2011). The higher frequencies, 285 and 369, correspond to the vortex shedding frequencies of the leeward corner vortices.
The proposed STE model shows a monotonic decay with no fundamental peaks. The energy levels are orders of magnitude lower than for the SJM mode. Hence the high energy fluctuations are attenuated in this valve model.
3.4.4. Mechanism of various instabilities of a BMHV and its attenuation
From the streamtraces of the SJM valve, we see the recirculation bubble form at the pressure side of the leaflet (figure 11,
$t=0.14$
s). The formation of the recirculation bubble thickens the width of the leaflet, which increases the expansion ratio and promotes the formation of corner vortices. Apart from the recirculation bubble, two recirculating eddies form alternately at the blunt leeward face of the trailing edge (figure 11,
$t=0.12$
s). At later times, these three vortices interact with one another (figure 11,
$t=0.1762$
s). At even later times, these form strong vortices, which then shed from the trailing edge. This mechanism is responsible for the flow disturbances and turbulence production (discussed later).

Figure 15. Contour plots of
${\partial p}/{\partial x}$
(mm
$\textrm {s}^{-2}$
) at the inner surface of the leaflet at time instants (
$a$
)
$t=0.10$
s and (
$b$
)
$t=0.18$
s, respectively, for the SJM valve, and (
$c$
)
$t=0.18$
s for the STE valve.
The formation of the recirculation bubble is related to the vorticity generation at the wall (i.e. at the inner surface of the leaflets), and for a static wall, it is in turn associated with the pressure gradient at the wall surface,
${\partial p}/{\partial x}=-{\unicode{x03BC}}\, {\partial \omega _z}/{\partial y}$
. Hence an adverse pressure gradient
${\partial p}/{\partial x} \ge 0$
leads to vorticity flux into the fluid and thus promotes the generation of recirculation region there. Figure 15(
$a$
,
$b$
) show contour plots of
${\partial p}/{\partial x}$
at the inner surface wall for the SJM valve. At
$t=0.10$
s (figure 15
$a$
), the pressure gradient is close to zero, and the streamtraces are attached to the wall with no formation of a recirculation bubble (figure 11). At
$t=0.18$
s (figure 15
$b$
), however, an adverse pressure gradient
${\gt } 25$
mm
$\textrm {s}^{-2}$
is seen at the downstream side of the surface (red), after which it becomes favourable again near the trailing edge. Thus at this time instant, a recirculation region is seen from the streamtrace plots for the SJM valve.
For the STE valves, the streamtraces are attached to the inner surface, and no recirculation bubble is present (see figure 15
$c$
). Recirculation bubble and corner vortex instabilities for two other valve models, arising from streamlining of the leading and trailing edges, are reported in Appendix A.2.
3.4.5. Vorticity dynamics
We further investigate the vorticity transport terms to understand the differences in three-dimensional vortex structures among these different valve models. The LA phase instantaneous flow field is considered here as the vortex activities are substantially pronounced during this phase. Let us now consider the vorticity transport equation

where
$( {\textrm D}/{{\textrm D}t}) ({\cdot })$
represents the material derivative operator. The
$z$
component of the vorticity transport equation reads

where the first two terms of the right-hand side combined denote turning or tilting of vortex tubes (oriented in the
$x$
and
$y$
directions) to the
$z$
direction, whereas the third term denotes stretching of the vortex tube (oriented in the
$z$
direction) into itself. In three-dimensional flows, these two kinds of terms result in the twisting and curving of the vortex rolls, and lead to the cascading of large-scale eddies into small-scale ones. The last term denotes the diffusion of the
$z$
directional vortex roll in the
$x, y, z$
directions.
Figure 16(
$a$
) represents the tilting of the
$x$
directional vortex roll in the
$z$
direction (the first term in the right-hand side of (3.4)). The SJM valve shows prominent
$x$
directional vortex roll tilting in the
$z$
direction in the wake of the leaflet (marked by red and blue), making the wake three-dimensional. For the STE valve, we see a negligible amount of vortex roll tilting. The behaviour of
$y{-}z$
tilting (not shown for brevity) is antisymmetric to that of
$x{-}z$
vortex roll tilting, as shown above.

Figure 16. (
$a$
) Vortex turning/tilting term 1 (
$x{-}z$
), (
$b$
) vortex stretching in the
$z$
direction, and (
$c$
)
$z$
-vorticity diffusion in all directions, units
$\text{s}^{-2}$
, at phase LA for the SJM and STE valves.
Next, we present the vortex stretching in the
$z$
direction (figure 16
$b$
). This term is prominent in the SJM valve, but shows reduced values in the proposed valve STE. Finally, we depict the distribution of diffusion of
$z$
vorticity in three dimensions (figure 16
$c$
). We can see the enhanced axial length of diffusion in the case of the optimised STE valve, which stabilises the jet, compared to the SJM valve.
3.5. Turbulence features
In earlier subsections, we discussed instantaneous features pertaining to different valve models. In this subsection, we will discuss turbulence features corresponding to the phase-averaged data at a few phases following § 2.9.
3.5.1. The Kolmogorov scales
The Kolmogorov time scale
$\tau _{\eta }$
is calculated from the mean dissipation rate of the fluctuating field components (
$\epsilon$
, shown later) and the kinematic viscosity of the blood, by
$\tau _{\eta }=(\nu / \epsilon )^{1/2}$
. The mean dissipation rate is given by
$\epsilon = 2 \nu\, \overline {s_{\textit{ij}} s_{\textit{ij}}}$
, where
$s_{\textit{ij}}$
is the strain rate tensor of the fluctuating components of velocity, and the overbar represents time averaging (over the PF phase in this case). In the present work, 22 cardiac cycles were considered for obtaining statistically meaningful phase-wise turbulent statistics from the simulation data, following § 2.9. Figure 17 presents the spatial distribution of the Kolmogorov time scale (
$\tau _{\eta }$
) for the SJM valve for the PF phase. In regions far downstream and upstream from the valve and its wake,
$\tau _{\eta }$
is generally high in value (
${\gt}10\ \text{ms}$
), whereas the lowest value of
$\tau _{\eta }$
is obtained just downstream of the valves (
$1{-}2.5$
ms). Considering the lowest value of
$\tau _{\eta }$
in the wake,
$\approx \!1$
ms, the integral time scale of largest eddies (also known as the turnover time scale) is obtained as
$\tau _{\eta _0}=\tau _{\eta }\, {\textit{Re}}^{1/2}$
, which for the PF phase is obtained as
$\tau _{\eta 0}=76$
ms. We may recall that the time window considered for sampling at the PF phase (§ 2.9) is
$T_{i=3}=3.67$
ms, which turns out to be much smaller than the integral time scale at this phase, indicating that the calculated phase-resolved turbulent statistics are free from large-scale fluctuations. Moreover, the time step for the simulations is chosen as
$10^{-5}$
s, which is much smaller than the Kolmogorov time scale. Hence the temporal evolution of the flow is captured completely, resulting in a DNS-type time resolution to that of the scale of the smallest eddies.

Figure 17. Kolmogorov time scales at peak flow rate for the SJM valve model.

Figure 18. Kolmogorov length scales obtained at peak flow rate for (a) SJM and (b) STE valve models.
Now we report the length scales of the smallest eddies, also known as the Kolmogorov length scale, for the SJM and STE valve models at the peak inflow rate. Figure 18 shows the Kolmogorov length scales
$\eta$
for both valve models at the PF phase. This scale is related to the mean dissipation rate
$\epsilon$
, and viscosity of blood
$\nu$
, as
$\eta = (\nu ^3 / \epsilon )^{1/4}$
. Qualitatively, for the SJM model,
$\eta$
reports lower values just downstream of the valves, while for the STE model, the lower
$\eta$
values are delayed to a greater downstream span. Also, the magnitudes of
$\eta$
are higher for the STE valve compared to the SJM valve. This qualitative description in the central plane (
$z=0$
mm) suggests a higher averaged value of
$\eta$
for the STE model rather than the SJM model, which may be linked to the potentially lower averaged blood damage in the case of the STE valve. We may quantify this through the averaged
$\eta$
in the wake of the leaflets due to the core and the lateral jets as 252.1
${\unicode{x03BC}} \textrm {m}$
and 392
${\unicode{x03BC}} \textrm {m}$
, respectively, for SJM and STE valves. A
$55.5\,\%$
higher averaged
$\eta$
is predicted in the wake of the leaflets for the STE valve compared to the SJM model. It is to be noted that the stress due to velocity fluctuations (Reynolds stress) also acts over this length scale (Liu et al. Reference Liu, Lu and Chu1999). Blood cells may be exposed to the fluctuating Reynolds stress terms, and potentially may undergo larger damage if the Kolmogorov length scale is in the same order of their dimensions (RBC diameter
$\approx 8$
${\unicode{x03BC}} \textrm {m}$
) (Ozturk, Papavassiliou & O’Rear Reference Ozturk, Papavassiliou and O’Rear2016). The higher averaged
$\eta$
in the wake of the leaflets for the STE model suggests that a lower mechanical load may be experienced at the scale of individual blood cells. Thus potentially lower blood damage will be seen compared to the SJM model.
The smallest fluctuations are at this scale, therefore stress due to velocity fluctuations (Reynolds stress) also acts over this length scale (Liu et al. Reference Liu, Lu and Chu1999). It is agreeable that higher levels of damage to blood elements (such as RBCs) occur due to the Reynolds stress terms if the Kolmogorov length scale is of the same order as the RBC dimension (
$\approx 8$
${\unicode{x03BC}} \textrm {m}$
in diameter). On the other hand, appreciable damage to blood elements is not expected from turbulent stresses for higher magnitudes of the Kolmogorov scale. In the literature, it is reported to be an order larger than the scale of RBCs (Liu et al. Reference Liu, Lu and Chu1999; Antiga & Steinman Reference Antiga and Steinman2009; Yun et al. Reference Yun, Dasi, Aidun and Yoganathan2014c
). The Reynolds stress influences the true mechanical load experienced by them, with an offset of one order higher than the nominal viscous stresses (Antiga & Steinman Reference Antiga and Steinman2009). Here, we compare the Kolmogorov scales for the SJM and STE valve models to explain the disparity in the levels of blood damage between the models, as reported in the previous subsection.
The Kolmogorov length scale may also be interpreted as the inverse of the fourth root of turbulent dissipation (
$\eta = (\nu ^3 / \epsilon )^{1/4}$
). Hence a smaller Kolmogorov length scale in the core region of the SJM flow indicates a higher rate of turbulent dissipation for it. The reduced dissipation may also be considered a major contributor to smaller energy loss through the STE valve model.
3.5.2. Turbulent kinetic energy and Reynolds stress

Figure 19. Contours of TKE at peak flow rate for (
$a$
) SJM and (
$b$
) STE valve models.

Figure 20. Contours of cross-stream (
$a$
)
$\rho \overline {u'v'}$
, (
$b$
)
$\rho \overline {u'w'}$
and (
$c$
)
$\rho \overline {v'w'}$
components of Reynolds stress at peak flow rate for SJM and STE valve models.
Figure 19 shows the contours of turbulent kinetic energy (TKE) at peak flow rate for SJM and STE valve models. The TKE is calculated as
$\text{TKE}=({1}/{2}) u{'}_i^2$
. Lower values of TKE are observed in the lateral and central jet regions for the STE model compared to the SJM model, indicating an abatement in turbulent fluctuations in the STE model. It is to be noted that TKE correlates with helicity (Gallo et al. Reference Gallo, Morbiducci and de Tullio2022). Thus high TKE zones may show more local circulation and increase residence time, which may lead to higher blood damage.
The cross-components of the Reynolds stress tensor are also negligible in the central region of the STE model compared to SJM (see figure 20 a–c). The value for the cross Reynolds term (
$\rho\, \overline {u'v'}$
) for the SJM valve matches well with Ge et al. (Reference Ge, Leo, Sotiropoulos and Yoganathan2005), with maximum value 83 Pa compared to 61 Pa in Ge et al. (Reference Ge, Leo, Sotiropoulos and Yoganathan2005). Also, the parallel Reynolds term (
$\rho\, \overline {u'u'}$
, contour not shown for brevity) for the SJM valve matches well, with maximum value 255 Pa compared with 226 Pa in Ge et al. (Reference Ge, Leo, Sotiropoulos and Yoganathan2005).
3.5.3. Effect of gradients of Reynolds stresses in phase-averaged velocity field
Considering a low-frequency mean flow field and high-frequency fluctuating velocity field, which are not correlated, the phase-resolved composite velocity field may be decomposed as

where
$\boldsymbol{U}=(U, V, W)$
are the time-averaged components of
$\boldsymbol{u}$
, and
$\boldsymbol{u}'=(u', v', w')$
are the fluctuating components in the
$(x, y, z)$
directions, respectively. Following Perkins (Reference Perkins1970), the mean out-of-plane vorticity (
$\varOmega _z$
) equation can be derived by eliminating the pressure from the
$x$
and
$y$
components of the Reynolds-averaged Navier–Stokes equations, which follows as

where
$\boldsymbol{\varOmega }=(\varOmega _x, \varOmega _y, \varOmega _z)$
represents the average rotation rate about the
$(x, y, z)$
directions, respectively, with

In (3.6), the left-hand-side term denotes the convection of mean
$z$
-vorticity by the mean velocity field. The first term on the right-hand side denotes the diffusion of mean
$z$
-vorticity due to fluid viscosity. The term
${\textrm{T}_{1}}$
,
$(\boldsymbol{\varOmega \boldsymbol{\cdot }\boldsymbol{\nabla }})W = \varOmega _x\, \partial W / \partial x + \varOmega _y\, \partial W / \partial y + \varOmega _z \,\partial W / \partial z$
, represents vortex skewing and vortex stretching in the
$z$
direction, due to the
$z$
component of velocity. The terms
${\textrm{T}_{2}}$
,
${\textrm{T}_{3}}$
,
${\textrm{T}_{4}}$
are present only when the flow is turbulent, where they represent the effect of time-averaged convection of turbulent vorticity by the turbulent flow field (terms
${\textrm{T}_{2}}$
and
${\textrm{T}_{4}}$
) and the time-averaged production of turbulent vorticity (term
${\textrm{T}_{3}}$
). The physical mechanism of these transport terms is given in detail in Perkins (Reference Perkins1970).

Figure 21. Contours of gradient of Reynolds stress terms in (3.6) at peak flow rate, for terms (
$a$
)
${\textrm{T}_{3}}$
, (
$b$
)
${\textrm{T}_{2}}$
, (
$c$
)
${\textrm{T}_{4}}$
, for SJM and STE models.
Figure 21 shows the contour plot of the terms
${\textrm{T}_{3}}$
,
${\textrm{T}_{2}}$
,
${\textrm{T}_{4}}$
for SJM and STE models at peak flow rate. Due to the convection driven by turbulent fluctuations (
${\textrm{T}_{2}}$
and
${\textrm{T}_{4}}$
), the mean vorticity spreads wider in the cross-plane. Therefore, the contribution of the turbulent fluctuations to the spread of large-scale vortices is larger in the SJM case. The term
${\textrm{T}_{3}}$
acts as an added rotational acceleration to the mean flow in the axial direction. Its increased values downstream help the STE flow to form helical Burgers-like vortices there. We can see the reduced production of fluctuating vorticity for the STE model compared to SJM (figure 21
$a$
) at the locations just after the valve where the helicity is reduced, and the axial jet is better stabilised for the STE valve model. Overall, it can be noted that the distribution of fluctuation terms in the STE valve model also contributes to its stabilised physiological favourable behaviour.
3.5.4. The TKE budget
The TKE budget relates the transport of TKE due to the mean flow velocities to the production and dissipation of the TKE. For a statistically steady flow, the expression of the TKE budget is


where the imbalance between the production (
$P$
) and dissipation (
$\epsilon$
) of TKE contributes to the advection of TKE in (3.8). The expressions of the terms are given in (3.9). Figure 22 shows the contour plots of the convection, production and dissipation of TKE at the peak flow phase for the control BMHV (SJM) model and final shape optimised STE valve. Very high magnitudes of convection of TKE are seen for the SJM model, which spatially initiates near the trailing edge of the leaflets (red) and spreads downstream into the sinus and converging part of the aorta. In contrast, the STE model shows negligible advection of TKE compared to the SJM valve. The production and dissipation of TKE are also much higher for the SJM valve compared to the STE valve.

Figure 22. Contours of (
$a$
) convection, (
$b$
) production and (
$c$
) dissipation of TKE at peak flow rate for SJM and STE valve models.
Zooming into the trailing edge of the leaflet for the SJM model (figure 23), we see the action of a recirculation bubble and corner vortices giving rise to the production of TKE. The streamtraces are for the averaged flow at this phase. The produced TKE then advects downstream through the free shear layer into the wake for the SJM model. Also, we see the dissipation of TKE surrounding the regions of the prominence of production of TKE in the SJM model. Due to the absence of a recirculation bubble and corner vortices, the production is reduced for STE.

Figure 23. Contours of production of TKE along with averaged streamtraces at peak flow rate for the SJM valve model.
3.6. Blood damage
In the previous subsections, we have seen that the proposed STE valve model outperforms the standard SJM valve owing to reduced disorganisation of flow structures, which is evident from higher EOA, smaller energy loss, lower TKE and stress levels between the two valves. In this subsection, we report the influence of leaflet shape on blood damage.
The platelet activation and RBC damage are both correlated with the cell loading, which is influenced by the local stress on blood elements and the exposure time (Han et al. Reference Han, Zhang, Griffith and Wu2022). We have used a linear accumulation model (Han et al. Reference Han, Zhang, Griffith and Wu2022) based on the Lagrangian trajectories of the particles to predict the blood damage in BMHV flows. A set of 500 particles, each positioned randomly in the cross-sectional plane at an axial plane
$1d$
upstream of the valve ring, is considered for seeding the initial positions. This seeding process is done for every 20 ms of flow until the end of systole. In total, 9000 particles were tracked until the end of the cardiac cycle, and their exposure levels of fluid stress and time were determined. Such passive fluid particles were considered since tracking individual RBCs or platelets is beyond the scope of the present study.
It is the viscous stresses that contribute to the platelet activation and blood damage (Ge et al. Reference Ge, Dasi, Sotiropoulos and Yoganathan2008). To this end, the viscous stresses are used in this subsection to predict the blood damage.
The viscous shear stress tensor has been interpolated along individual particle trajectories, and an equivalent scalar shear stress is obtained based on the modified von Mises criterion by Han et al. (Reference Han, Zhang, Griffith and Wu2022), given as

where
$\tau _{\textit{ij}}$
are the components of the shear stress tensor. A linear damage model (Han et al. Reference Han, Zhang, Griffith and Wu2022) is used in this study to quantify the blood damage index (BDI) at platelet activation state (PAS) as

where
$\tau$
is the scalar shear stress, and
$\Delta t$
is the numerical time step. All necessary quantities are interpolated along the trajectory of selected particles from their initial seed position by
${{\textrm d} \boldsymbol{x}_{\!p}}/{{\textrm d}t}=\boldsymbol{u}$
.
Other blood damage models stem from experimental tests based on empirical fittings. Giersiepen et al. (Reference Giersiepen, Wurzinger, Opitz and Reul1990) summarised many such models as
$IH=\Delta Hb / Hb = C t^{\alpha } \tau ^{\beta }$
, where the exponents of exposure times
$t$
and scalarised shear stress
$\tau$
give rise to higher-order models. Other models, such as the damage accumulation model, which selectively apply for activated and non-activated particles, are given by Alemu & Bluestein (Reference Alemu and Bluestein2007) and Yeleswarapu et al. (Reference Yeleswarapu, Antaki, Kameneva and Rajagopal1995), and for PAS by Soares et al. (Reference Soares, Sheriff and Bluestein2013).
Figure 24 shows a histogram of the BDI of 9000 particles at the end of the cardiac cycle for the valve models considered in this study. The SJM model predicts a high BDI distribution overall. For STE valves, however, the number of particles is higher at lower damage levels, and the number of higher damage level particles is also lower. This suggests an overall lower blood damage for the STE valve design. Statistical values of BDI for the sample of 9000 particles are reported in table 6. A
$16.7\,\%$
lower mean BDI is reported for the STE valve model compared to the SJM model. Parameters such as average deviation, standard deviation, median, mode and maximum BDI values are also lower in the case of the STE valve. The BDI values for two other valve models are reported in Appendix A.3.
Table 6. Statistics of BDI (in Pa s) for SJM and STE valve models.


Figure 24. The BDI (in Pa s) histogram for 9000 particles obtained after an entire cardiac cycle for SJM and STE valves.
Particle pathline patterns are crucial for suggesting the stress experienced and the exposure times to said stress levels. In that regard, readers are referred to Sarkar et al. (Reference Sarkar, Sharma, Chakraborty and Roy2024b ) for various pathlines typical in BMHV flows.
We mention that the hinge is modelled as a simple pin joint in our work, and utilising the actual hinge mechanism (with the leaflet ear and butterfly recess) may result in a higher damage for the blood elements passing through the hinge region. A recent study, however, reported that the total blood damage through the bulk region is several times higher than due to the hinge flow (Hedayat & Borazjani Reference Hedayat and Borazjani2019). Hence in this study, a simplified pin joint hinge mechanism is used to compare the blood damage effects due to the different valve models.
3.7. Effect of valve design on flow dynamics
In this subsection, we summarise the various mechanisms at play that promote smoother flow features and lower blood damage for the proposed valve design. For the SJM valve, we see that the recirculation vortex (near the trailing edge) plays a very important role in promoting large-scale flow disturbances. The genesis of the recirculation vortex is attributed to an adverse pressure gradient on the inner walls of the leaflets. This recirculation region is associated with high levels of turbulence production and increased helicity. These flow disturbances lead to production of turbulence. Overall, the fluctuations arising in the SJM valve contribute to higher stress accumulation along blood cell pathlines, which may potentially lead to significant blood damage.
The streamlining of the leading and trailing edges results in suppression of the recirculation vortex as well as the trailing edge vortex, which induces lower levels of turbulence, and the jet is much more stabilised. The resulting delayed formation of Burgers vortices gives rise to lower levels of helicity. Such design choices result in stabilised flow and lower blood damage with better leaflet kinematics. The proposed valve design yields a reduced flow energy loss along with smooth leaflet kinematics too.
4. Conclusions
A mechanical heart valve is aimed to act as a durable prosthetic device designed to replace damaged heart valves, ensuring proper blood flow, and preventing complications such as heart failure and stroke. Recent advancements in mechanical heart valve technology include bileaflet designs, featuring two semicircular leaflets that open and close to regulate blood flow. This design aims to optimise flow control and minimise turbulence compared to earlier mechanical valves, but necessitates lifelong anticoagulation therapy to prevent blood clots, which increases bleeding risks and requires regular monitoring. From a fluid dynamics perspective, these risks are associated with the formation of turbulent vortices, which can lead to adverse effects such as haemolysis, increased thrombosis risk, platelet activation, and valve dysfunction due to wear. These issues drive design innovations to mitigate turbulence-induced complications. Excessive turbulence may also cause energy loss, increase cardiac workload, damage endothelial cells, and trigger paravalvular leaks, impairing valve efficiency. Furthermore, turbulence can heighten the risk of infective endocarditis by causing endothelial trauma.
To address these challenges, we investigated a modified design of the St Jude Medical (SJM) valve, a widely used bileaflet mechanical heart valve for aortic and mitral replacements. Using direct numerical simulations (DNS), we evaluated the performance of this redesigned valve, which incorporates streamlined leading and trailing edges. The DNS results revealed that the streamlined edges effectively minimise abrupt changes in blood flow direction, prevent flow separation, and reduce the formation of turbulence-inducing vortices and eddies. The pressure gradient along the valve surface was shown to be significantly smoothed, and the flow structures stabilised, resulting in the elimination of recirculation bubbles and the disappearance of the corner vortices. The resulting geometry yielded a much-stabilised jet with reduced helicity and turbulent kinetic energy production. Our results also revealed that the gradients of the turbulent stresses contribute to the breakdown of the large-scale jets, and once the turbulence is controlled, the jets are better stabilised. Finally, we showed that the blood damage indices are also moderated through the proposed design.
Compared to the conventional SJM valve, the modified design thus promotes smoother flow transitions across valve surfaces, and reduces flow disturbances. The findings emphasise the benefits of edge streamlining in decreasing obstruction, drag, and wake turbulence downstream. This results in smoother blood flow, balanced pressure gradients, reduced shear stress, improved haemodynamic efficiency, and a lower risk of complications such as haemolysis and thrombosis. These advancements represent a significant step forward in mechanical heart valve technology.
Acknowledgements
The authors acknowledge PARAMShakti, IIT Kharagpur and PARAMSiddhi, CDAC, Pune (National Supercomputing Mission, Government of India) for providing computing facilities. S.C. acknowledges SERB, Department of Science and Technology, Government of India, for the Sir J.C. Bose National Fellowship.
Funding
The present research is funded through SERB Project CRG/2019/00265.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The data that support the findings of this study are available from the corresponding author, S.R., upon reasonable request.
Appendix A. Additional valve designs
Figure 25 illustrates the cross-section of the two intermediate valves. The first intermediate leaflet design has features from modifying the SJM leaflet design by streamlining the leading edge (figure 25 a), thus departing from the blunt leading edge of the SJM valve model, which closely resembles the modified leaflet design proposed by Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019) and Zolfaghari et al. (Reference Zolfaghari, Kerswell, Obrist and Schmid2022). This design is the sharp leading edge (SLE) valve. The second intermediate design (figure 25 b) entails a tapered trailing edge (TTE) of the leaflets while maintaining a blunt leading edge similar to the SJM valve.

Figure 25. (
$a$
) The SLE model. (
$b$
) The TTE model.
A.1. Streamtraces
Figure 26 reports the streamtraces due to the SLE valve model. The SLE model is inspired by the works of Zolfaghari & Obrist (Reference Zolfaghari and Obrist2019) and Zolfaghari et al. (Reference Zolfaghari, Kerswell, Obrist and Schmid2022). The profiling of the leading edge to sharp from blunt suppresses the recirculation bubble and thus the formation of the recirculation region at the inner side of the leaflet. Still, because of the bluntness of the trailing edge, the formation of a recirculation region at the upper and lower sides of the leeward face occurs, eventually detaching and dissipating into the flow. Due to this, the streamtraces downstream of the leaflets show undulation. The streamlined trailing edge of the TTE valve model negates the formation of trailing edge vortices but still suffers from a recirculation bubble due to its leading edge bluntness, hence the undulations in the streamtraces are observed (figure 27).

Figure 26. Streamtrace patterns of the SLE valve model.
A.2. Mechanism of various instabilities of the BMHV and its attenuation
The recirculation bubble and associated zone of recirculation are also seen for TTE valves (see figure 28
$b$
), where the zones of positive
$ {\partial p}/{\partial x}$
are separated by a zone of favourable pressure gradient (shown in blue); this is evident from the streamtrace patterns (
$t=0.1375\,\text{s}$
to
$t=0.1820\,\text{s}$
, figure 27), with the advection of the recirculation bubble downwards from the front portion of the lower surface of the leaflet to the trailing edge of the lower surface, and eventual detachment from the valve body.
For the SLE valve, the wall pressure gradient is adverse in the initial region where material has been removed (the sharp portion of leading edge) due to the flow entrance effect (see figure 28
$a$
), but in the distal semi-elliptical region, pressure gradient is negative. Due to the inertia of the flow and the downward favourable pressure gradient, the streamtraces remain attached to the surface in the case of the SLE valve.

Figure 27. Streamtrace patterns of the TTE valve model.
A.3. Blood damage
In table 7, we report the BDI statistics for four different valve models. It may be observed that the SLE valve shows some improvement over the standard SJM model. However, the streamlined STE leaflet model outperforms the other valve models.
Table 7. Statistics of BDI (in Pa s) for various valve models.


Figure 28. Contour plots of
${\partial p}/{\partial x}$
(mm
$\textrm {s}^{-2}$
) at the inner surface of the leaflet at time instant
$t=0.18$
s for (
$a$
) SLE and (
$b$
) TTE valves.