A 3D computational model similar to that in our previous work (Chen and Luo, 2018; Chen and Luo, 2020) is adopted and shown in Figure 1. The aorta is simplified to a cylindrical tube of diameter D = 2.1 cm and length L = 19 cm. It includes three-lobed dilation to represent the aortic sinuses, of which the geometry and dimensions are based on physiological measurements of the human aortic root (Swanson and Clark, 1974; Reul et al., 1990). A tri-leaflet aortic valve is positioned within the sinus region, and its three flexible leaflets can deform independently from each other. Despite the variability of the anatomy of the human aorta, simplified computational domains similar to this are often used for the FSI study of the native aortic valve and its prostheses (Borazjani, 2013; Hsu et al., 2015; Kamensky et al., 2015; Griffith, 2012; Marom et al., 2013; Mao et al., 2016). Similar to previous FSI studies (Borazjani, 2013; Kamensky et al., 2015), a transient transvalvular pressure load is applied at the inlet of the aorta tube to drive the blood flow (see Figure 2). The exit pressure at the outlet is 0 kPa. This pressure drop is consistent with previous studies of the human aorta (Hole, 1996; Kim et al., 2008). Since the aortic wall is assumed to be rigid in the model, which is a limitation to our study, the specific reference pressure at the outlet does not matter, and the FSI is purely driven by the pressure difference between the two ends.

For spatial discretization, the aorta wall is divided into 20,735 triangular surface elements with local refinement in the sinus region, where the element size is about 0.3 mm. A uniform thickness of 0.3 mm is assumed for each leaflet, which consists of a total of 539 FE serendipity (20-node Hexahedron) elements and 4,021 nodes. The meshes of the aorta wall and leaflets are separate, and they intersect each other without necessarily sharing the nodes. The aorta wall is assumed to be rigid, while the leaflets can undergo free deformations. The discretization of 20,735 triangular surface elements of the aorta wall serves as a physical boundary in the immersed boundary method–based simulation, where the no-slip and no-penetration conditions are applied when solving the Navier–Stokes equations. On each leaflet, 830 nodes (red markers) located in both the commissure region of two neighboring leaflets and along the base are fixed, while 903 nodes (blue markers) consisting of the prescribed contact region are used to prevent leaflet inter-penetration, as shown in Figure 1C.

Aortic valves are known to be nonhomogeneous (Kuang et al., 2017; Oveissi et al., 2020), and different constitutive models have been proposed to investigate the mechanical behavior and the failure mechanisms of the aortic valve (Weinberg and Kaazempur-Mofrad, 2005; Mao et al., 2016). In this work, the hyperelastic Saint Venant–Kirchhoff model is adopted to represent the tissue behavior of the valvular leaflets. The constitutive relationship of the Saint Venant–Kirchhoff model can be expressed as follows: σ i j K = D i j m n E m n (1) where D i j m n is the elastic matrix with 12 non-zero elements that depend on Young’s modulus E and Poisson’s ratio ν s , and E m n is the Lagrangian strain tensor (Tian et al., 2014).

The leaflet dynamics is governed by the following equation: ρ s d 2 u i d t 2 + η d d u i d t = ∂ σ i j ∂ x j (2) where u i is the nodal displacement, η d is the damping coefficient representing structural damping of the tissue, and σ i j is the Cauchy stress tensor. The density of the leaflets is ρ s = 1 g/cm³, and η d is chosen to be 100 g/cm³·cs to ensure the reasonable time scale for valve opening and closing (Chen and Luo, 2018; Chen and Luo, 2020). Three different values of Young’s modulus, E = 375 , 750 , and 1500  kPa are selected to represent the difference in the stress–strain relationship between the more compliant PVP and the stiffer pericardium measured from planar uni-axial tests (see Figure 3). The Poisson’s ratio ν s = 0.4 is used for all three cases (Hsu et al., 2015; Kamensky et al., 2015).

The blood is assumed to be Newtonian and incompressible. The governing equation of the flow is the unsteady Navier–Stokes equation as follows: ∂ v i ∂ x i = 0 (3) ρ ( ∂ v i ∂ t + ∂ v i v j ∂ x j ) = − ∂ p ∂ x i + μ ∂ 2 v i ∂ x j ∂ x j (4) where v i is the velocity component, p is the pressure, ρ is the blood density, and μ is the dynamic viscosity. No-slip and no-penetration boundary conditions are imposed on the aorta wall and the leaflet surface. The fluid domain is a 19 × 4.4 × 4.4 cm³ rectangular bounding box and is divided by a 400 × 130 × 130 non-uniform Cartesian grid. A fine resolution with ∆x = 0.025 cm and ∆y = ∆z = 0.034 cm is used in the region around the aortic valve. A mesh refinement study has been conducted to validate the accuracy of the current mesh, in which the flow domain is divided by a 500 × 260 × 260 non-uniform Cartesian grid with a resolution around the leaflets of ∆x = 0.015 cm and ∆y = ∆z = 0.017 cm, and the mesh for the finite element model of the valve is increased by approximately five times. Comparisons of the leaflet dynamics, hemodynamic resistance of the aortic valve, transient flow rate, and pressure distribution on the leaflet surface have shown that the results between the baseline and the refined mesh are in excellent agreement (Chen and Luo, 2020). The density and dynamic viscosity of the blood are ρ = 1 g/cm³ and μ = 0.005 Pa s, respectively (Chen and Luo, 2018; Chen and Luo, 2020).

We applied an in-house computational approach that was previously developed for simulating biological systems involving large deformations to solve the heart valve FSI problem (Luo et al., 2012; Tian et al., 2014). In this approach, the flow and solid solvers are arranged in a partitioned manner such that the flow is simulated using an accurate direct-forcing IB method based on a Cartesian grid, and the solid is solved using a nonlinear FE method. Strong FSI coupling is achieved by iterating the two solvers while communicating the transient position and velocity information at the fluid–structure interface until convergence. Detailed explanations of the numerical algorithms and their parallel implementation are included in our previous work (Chen and Luo, 2018; Chen and Luo, 2020). Each cardiac cycle has a duration of T = 0.86 s, which corresponds to a heart rate of 70 beats per minute. For convenience, we used centisecond (cs) as the time unit thereafter. To ensure numerical stability of the FSI coupling, the time step used for the flow solver is ∆t = 4.0 × 10^–3 cs. The time step for the structural simulation is smaller, ∆t = 5.0 × 10^–5 cs, so that each FSI step contains 80 sub-steps for the solid.

To accelerate the costly FSI simulation, a parallel computing technique based on a domain decomposition strategy has been implemented on the flow side. On the structure side, multiple OpenMP threads are forked on a single CPU processor, and the data communication between the fluid and solid sides is achieved via message passing interface (MPI) calls. Details of the implementation have been discussed in our previous work (Chen and Luo, 2018; Chen and Luo, 2020). For the FSI simulations in this study, the rectangular flow computation domain is divided into 13 subdomains in each of the y- and z-directions, which requests a total of 169 CPU cores. A separate CPU core handles the computation of solid mechanics, and the FE solid solver is parallelized using 16 OpenMP threads. Each heart valve FSI simulation takes approximately 50    h for one cardiac cycle on Stampede 2 partitions at the Texas Advanced Computing Center (TACC).

The contact force was calculated to prevent the inter-penetration between leaflets, especially during the closing phase. At each time step, the contact distance d of the nodes within the prescribed contact detection region (blue markers in Figure 1C) was first calculated by projecting it onto the surface of its neighboring leaflets. d 0   is a prescribed threshold of distance. When d < d 0 , the contact algorithm is activated, and the contact force is calculated as follows: f c n t = − f e x t − k ( d − d 0 ) (5) where f c n t is the nodal contact force, f e x t is the hydrodynamic force from the blood flow on the node, and k is the contact stiffness. It should be noted that the external load f e x t   is canceled out, and a net force of magnitude − k ( d − d 0 ) is added to prevent inter-penetration of colliding leaflets when d < d 0 . It is noted that the contact force vanishes outside of the contact distance (when d > d 0 ). Here, we set k = 0.04  g/cs² and d 0 = 0.08  cm so that the leaflets are stopped without collision, and the gap between them is below one fluid cell width (Chen and Luo, 2018; Chen and Luo, 2020).

In Figure 4, we presented the transient flow rate Q and valve geometric opening area (GOA) for the three FSI cases. The flow rate is calculated at the outlet. The peak flow rate reduces from 527 and 492 to 459  ml/s as Young’s modulus increases from 375 to 750 and 1500  kPa. Integrating the transient flow rate in time, we obtained the corresponding stroke volume, and they were 112 , 99.6 , and 91.4    ml, respectively, which are within the physiological range for normal adults (Griffith et al., 2009; Borazjani, 2013). During the valves’ rapid closure, the negative dip flow rates were − 144 , − 131 , and − 102  ml/s for E = 375 , 750 , and 1500  kPa, respectively. Valve reverberation during rapid closure corresponds to the clinical second heart sound (S2) and has an indication of heart function. Stronger reverberation is found for the more compliant aortic leaflet valve. We also calculated the regurgitation volume, which represents the blood leakage backward through the aortic valve toward the inlet during diastole. For E = 375 , 750 , and 1500  kPa, the regurgitation volume is less than 0.3  ml per cardiac cycle.

The GOA from the FSI simulation is presented in Figure 4B for all three cases. This area is calculated by projecting the valve in the axial direction and finding the opening area. From the GOA history, we observed that all the aortic valves opened rapidly in 5 to 6  cs and experienced rapid closure within 3 to 4  cs. The more compliant aortic valve with a smaller E = 375    kPa opened slightly faster while shutting down slightly delayed by scrutinizing the GOA profile. Also, a larger GOA value after reaching full open is observed for the more compliant leaflets. The peak GOA values were about 2.3 , 2.2   , and 2.0    cm² for E = 375 , 750   and 1500  kPa, respectively, at the peak systole. In addition, a small oscillation is observed for the more compliant valve ( E = 375  kPa) during the peak systole. This oscillation is related to the fluttering motion of the leaflets’ free edge after it reached the fully open stage, and it is negligible when compared with the more compliant aortic leaflet valves with thicknesses of 0.05, 0.08, and 0.1 mm, as shown in Ref. 16. All three valves closed properly without leakage since their GOAs reached zero during diastole.

Figure 5 shows the opening and closing phases of the aortic valves for all three cases. The opening process for the more compliant aortic valve ( E = 375 and 750    kPa) is faster than its stiffer counterpart E = 1500  kPa (see t = 4.4  cs). At peak systole t = 22.0  cs, all three aortic valves reached the fully open stage. The closing phase of the more compliant leaflet valve is slightly delayed (see t = 34.8  cs), which can also be seen from the GOA history shown in Figure 4B.

In Figure 6, we plotted the normalized hemodynamic force along the x-direction experienced by the aortic valve. During the systolic phase, the more compliant aortic leaflet valve ( E = 375  kPa) experienced a lower axial hemodynamic force when compared with its stiffer counterpart. The momentum balance study presented in our previous work (Chen and Luo, 2018; Chen and Luo, 2020) has revealed the following approximation: Δ P A ≈   ρ L Q ˙ + F v a l (6) where Δ P is the transvalvular pressure drop shown in Figure 2, A = π D 2 / 4 is the cross-section area, Q ˙ = d Q / d t is the time derivative of the flow rate, and F v a l is the total axial hemodynamic force on the leaflet surfaces (both aortic and ventricular sides included). Given the pressure loading Δ P A , the smaller hemodynamic force F v a l of the more compliant aortic valve in systole corresponds to a greater acceleration of the fluid column ρ L Q ˙ inside the aortic tube, which resulted in a higher transient flow rate Q   for the higher compliance aortic valves, as presented in Figure 4A. This result means that the more compliant aortic leaflet valve takes more advantage of the transvalvular pressure drop for blood delivery. Meanwhile, the dip of negative hemodynamic force during the rapid closure is larger for a higher compliant valve to stop the reversal of flow. The magnitudes of the forces are − 5.8 , − 5.7   , and − 4.8 for E = 375 , 750   , and 1500  kPa, respectively.

Figure 7 shows the maximum principal stress (MPS) of the aortic leaflet for all the three cases during the opening phase ( t = 4.4  cs), at peak systole ( t = 22.0  cs), and at the instant of peak “water-hammer” which corresponds to the maximum negative hemodynamic impact shown in Figure 6A. The stress distributions on the leaflet surfaces are qualitatively similar among the three valves, while their magnitudes are significantly different. The leaflet stress for the more compliant aortic valve is substantially lower during both the opening and peak systole phases (the first two columns in Figure 7 can be referred). For example, at peak systole, the maximum value of stress is around 271  kPa for E = 375  kPa, while it increases to 383 and 540  kPa for E = 750 and 1500  kPa, respectively. The “water-hammer” effect emerges because of the aortic valve’s rapid shutdown and the ensuing impingement of fast-stopped blood on the leaflets, which creates high stress on the aortic side of the leaflets (the third column of Figure 7 can be referred). In contrast to the systole phase, at the instant of peak “water-hammer”, the strength of “water-hammer” is stronger for the more compliant valve, which produces a larger hemodynamic force on the valve, and the leaflet stress for the more compliant valve is higher than its stiffer counterparts. The maximum value of stress is about 550 kPa for E = 375  kPa, and it reduces to about 450 and 400 kPa for E = 750 and 1500    kPa in the belly region, respectively. In addition, the more compliant aortic leaflet valve is also found to be pushed more upstream toward the inlet due to this strong impact, which also generates higher stress along the fixed attachment.

In the present study, we simulated the stress distribution of the leaflets made of bovine pulmonary visceral pleura (bPVP). The bPVP is the layer of the serous membrane overlying the lungs and is composed of abundant collagen and elastin. Due to the larger proportion of elastin, PVP is more elastic and resilient when compared with porcine heart valves and bovine pericardium, which may render bPVP potentially advantageous for mitigation of calcification due to its much smaller Young’s modulus than that of porcine heart valves and bovine pericardium, which are collagen-dominant biomaterials and have been broadly utilized to construct bioprosthetic heart valves (BHVs) (Lindberg and Badylak, 2001; Maestro et al., 2006; Manji et al., 2012; Gauvin et al., 2013; Kheradvar et al., 2015; Dasi et al., 2017). Our simulation demonstrates that the mechanical stresses in bPVP leaflets are much smaller than those in the biomaterial where the Young’s modulus is similar to that of the bovine pericardium during the systole, even though some increase is also found during the valves’ rapid close and the ensuing “water-hammer” effect. The overall smaller mechanical stresses in bPVP leaflets may mitigate the degradation and calcification of leaflets in the bPVP prosthetic valve. This hypothesis remains to be validated in preclinical large animal models and clinical studies.

The aortic valve performance is determined by both the blood flow and valve structure since the instantaneous leaflet shape and motion are a result of the two-way interaction between them. In our previous work (Chen and Luo, 2020), a dimensionless parameter named normalized bending rigidity E B ∗ = E B / ( Δ p R 4 / L ) is defined with an optimal range for proper valve opening performance. E B = E h 3 / 12 ( 1 − ν 2 ) is the bending rigidity, R is the aorta radius, and Δ p / L is the pressure gradient applied to drive the blood flow through the aortic valve. Given the driven pressure gradient Δ p / L , Young’s modulus E , and Poisson’s ratio ν , the optimal range of E B ∗ ∈ [ 0.003 ,   0.4 ] is determined by investigating the aortic valve performance with leaflet thickness h varying from 0.005 to 0.8    mm. In the current work, the leaflet thickness is fixed with h = 0.03  mm, while Young’s modulus E is varied from 375 and 750 to 1500  kPa to represent the material property difference between the more compliant PVP and its stiffer counterpart of porcine/bovine pericardium. Following the definition of E B ∗ proposed in Ref. 16, we obtained E B ∗ = 8 × 10 − 3 , 1.6 × 10 − 2 and 3.1 × 10 − 2 for E = 375, 750, and 1500 kPa, respectively. We thus observed that for all the three cases studied in this work, their E B ∗ fall within the optimal range of E B ∗ between 0.003 and 0.04 . Thus, the three aortic valves in this work can perform properly without experiencing either evident leaflet fluttering ( E B ∗ too small) or opening difficulty ( E B ∗ too large). This is also confirmed by our results on the flow rate Q , valve opening area GOA, and leaflet deformation in the previous section.

In our previous work (Chen and Luo, 2020), the optimal range of E B ∗ is established based on the proper valve opening area (GOA). Excessively small and large E B ∗ ∉ [ 0.003 ,     0.4 ] are excluded to prevent severe leaflet fluttering and valve opening difficulty. In this work, we extended our previous study by providing an additional dimension of investigating the relationship between the aortic leaflet stress and normalized bending rigidity E B ∗   within the optimal range. This dimension is important since a recent FE study has pointed out that the significantly elevated mechanical stress is strongly correlated with the regions with a high risk of calcium buildup (Qin et al., 2020). Even though E B ∗ for all three cases in this work fall within the optimal range E B ∗ ∉ [ 0.003 ,     0.4 ] , the difference in the magnitudes of the aortic leaflet stress (maximum principal stress) is evident (as shown in Figure 7). During the systole phase, the stress is substantially lower for the more compliant leaflet valve, which suggests that PVP could serve as a potentially advantageous biomaterial over traditional pericardium for the construction of BHVs in the systole period. However, during the fast-closing phase, the more compliant aortic valve experiences a stronger “water-hammer” effect with a high peak hemodynamic force impinging on the leaflets, which pushes the aortic valve more toward the left ventricle side and results in an increase of stress. Thus, the potential advantage of PVP might be compromised due to the “water-hammer” effect during the fast-closing phase. Future studies will be conducted to investigate the performance of PVP in our follow-up effort.