Heart valve diseases like valvular stenosis or valvular insufficiency are a significant cause of global mortality, especially in Western countries. The human heart has two types of valves: mitral and aortic valves in its left compartment, which regulate unidirectional blood flow during the cardiac cycle. The mitral valve controls the blood flow into the ventricle and has two leaflets, and the aortic valve controls the blood flow into the aorta and has three leaflets. These valves operate passively by their leaflet movements (deformations) merely by the pressure gradient in the blood flow during a cardiac cycle. Because of various diseases caused by different agents such as bacterial infection, congenital heart defects, and rheumatic fever, the natural heart valve is subject to replacement by animal or prosthetic heart valves. Animal originated valves, which are called biological heart valves, are mainly bovine or porcine valves which are chemically treated to adapt human body (Liao et al., 1992). These valves are susceptible to tissue degeneration and calcification and, as a result, subject to failure 10–15 years after implantation (Siddiqui et al., 2009). Another problem with these valves is the fact that they have no uniform structure and properties like thickness and mechanical strength due to the creatures’ differences. Mechanical prosthetic heart valves solve this shortcoming of biological valves, and because of their durability, some studies (Zakaria et al., 2017) show that almost 60% of worldwide annual heart valve replacements use mechanical heart valves.

Mechanical heart valve designs have progressed in recent decades. The early designs of the mechanical valves, caged ball valves and tilting disc valves, caused several problems such as hemolysis of red blood cells and platelet activation because of their blood flow patterns (Black et al., 1983). Polymeric heart valves have gained the researchers’ interest because of their similarity with the natural heart valve in their mechanical behavior and controllability on the biological and biocompatibility characteristics. Different polymers are used as the base material for synthetic heart valves such as polyester (Vaesken et al., 2014; Meddahi-Pelle et al., 2021; Zhou et al., 2021), polyurethane (Gasparotti et al., 2018; MohmadSaberi, 2020; Ataee et al., 2021; Kazeroni et al., 2021), thermoplastic polyurethane (Fallahiarezoudar et al., 2017a; Fallahiarezoudar et al., 2017b; Zuke, 2017), and ultrahigh–molecular weight polyethylene (Thomas and Jayabalan, 2008; Zhou et al., 2021).

Designing an appropriate synthetic heart valve requires assessing its performance during a cardiac cycle from different contexts. From the mechanical design point of view, both leaflet deformation/stress and blood flow pattern should be studied. As the leaflet is under fluctuating stress, large fatigue life requires lower stress values in addition to the flexibility to open and close the flow path as the main function of the heart valve when we consider the deformation requirements. Compared with in situ investigations, numerical simulations give in-depth, high-resolution values in parameters as functions of time and position and, as a result, complete information on the performance of the synthetic heart valves under operation. Because of the nature of the valves, fluid–structure interaction (FSI) simulations with two-way coupling are the best analysis methods for this scenario (Bailoor et al., 2021; Geronzi et al., 2021). Leaflets and the vessel are modeled as solid parts with Cauchy equations of motion, while the blood is modeled as the fluid part with Navier–Stokes equations of motion, having deformable boundaries (the solid part) with two solid and fluid equations being solved simultaneously. Borazjani (2013) developed an immersed boundary-finite element model to study the FSI of a heart valve and studied the blood flow pattern in the mechanical and bioprosthetic heart valves. Gilmanov and Sotiropoulos (2015) compared the blood flow pattern in an aorta with bicuspid and trileaflet valves.

One of the early works on the analysis of the bioprosthetic heart valves is the study by Arcidiacono et al. (2005). They determined the stress distribution on the leaflets made with the bovine pericardium during a full cardiac cycle. They assumed a linear elastic orthotropic material model for the leaflets and concluded that the orthotropy affects the leaflet performance and its stress values. Driessen et al. (2007) studied the deformation of tissue-engineered heart valve leaflets under pressure and culturing conditions. They concluded that stresses increase with the culture time, while strains decrease. Kim et al. (2007) studied the dynamic response of a bioprosthetic heart valve using a shell model with the Fung constitutive model for leaflets. Lee et al. (2020) studied a porcine- and a bovine-based bioprosthetic heart valve using FSI simulations and compared the results with a model under a pulse duplicator. They showed a good agreement between experimental and numerical results for the leaflet deformation and dynamic valve area. Kerr and Gourlay (2020) investigated an expandable valve design using FSI simulations. Park et al. (2021) reviewed various study techniques of the prosthetic heart valves and categorized them. Nestola et al. (2021) proposed an FSI model to study the dynamic behavior of the bioprosthetic heart valves with the goal of optimizing the stress distribution on the leaflets. Recently, a review on the numerical studies published so far is presented by Abbas et al. (2022). Hsu et al. (2015) studied FSI of a prosthetic heart valve using parametric design and compared the results of the FSI simulation with those of structural dynamic simulation of the valve under pressure loading.

Gulbulak et al. (2020) studied the effect of the generating curve of the leaflet geometry on the orifice area of the polymeric heart valve. They used a three-point generating curve of the leaflet geometry and obtained their best parameters using the Taguchi method.

In this article, a prosthetic valve design is presented by studying the effects of the geometric parameters on the overall performance of the valve. To this end, first, the effects of the leaflet curvature and its inclination on the valve opening area and maximum von Mises stress of the leaflet are obtained by some 2D simulations. Then, a 3D model is constructed based on the results obtained from the 2D analyses. At this stage, some incremental changes are made on the design variables, and the simulation of the valve is done for a fraction of the total time of the cardiac cycle. The final value is fixed based on the obtained results of the incremental studies. Finally, the full design is analyzed, and the performance parameters are compared with the published articles.

Mechanical performance of the prosthetic heart valves depends on the design parameters of the leaflets such as the curvature, height, and the slope at the skirt location. To come up with an appropriate design, first, we conduct a preliminary 2D study on the effects of the leaflet slope at the skirt on the mechanical parameters such as the valve opening and stress values on the leaflets. It is noted that fully coupled FSI analyses are carried out using the Comsol Multiphysics finite element software package. Moving mesh on the fluid domain is controlled by arbitrary Lagrangian–Eulerian (ALE) formulation. Comsol couples the solid and fluid domains on the common boundaries by transferring displacements, velocities, and pressure loads at every step taken by the solver to solve the governing equations of the solid and fluid domains. According to the loading conditions, we have applied the blood velocity at the inlet of the vessel, as shown in Figure 1. At the outlet, the pressure is set to be zero. The valve material is chosen to be the same as chosen by Luraghi et al. (2017), poly(styrene–ethylene–propylene–styrene) (SEPS) block copolymers with 22% polystyrene fraction for comparison purposes. The polymer has a density of 830 kg/m³, Young’s modulus equal to 3.2 MPa, and Poisson’s ratio equal to 0.4. The blood is modeled as the Newtonian fluid with a density of 1,060 kg/m³.

Four models of the leaflets with different slopes and curvatures are considered, as shown in Figure 2. Leaflet designs L1 and L2 have linear profiles with different slopes at the skirt: 40° for L1 and 50° for L2. L3 design has a circular curvature, and L4 has a parabolic curvature to give insights on the effects of the curvature and slope on the performance metrics of the prosthetic heart valve.

Maximum opening of the leaflets is shown in Figure 3. It follows that the leaflet design L2 results in a larger opening, compared with other designs. In addition to the temporary value of the valve opening, the averaged opening area of the valves is shown in Figure 4. The results show that a larger slope at the skirt of the valve results in a larger valve opening area. In addition, for leaflets to have large opening area, the profile should be linear, at least at the skirt area of the valve.

Maximum von Mises stress on the leaflets as a function of time during a cardiac cycle is illustrated in Figure 5. From the stress point of view, it is observed that leaflets with linear profile have smaller stress values.

From these 2D simulations, we can conclude that the leaflet should have a larger slope at the skirt area of the valve and as linear as possible at that location so that the valve has a larger opening area and smaller stress values during the operation of the valve.

Based on the published studies and the 2D simulation results in the previous section, a prosthetic heart valve with the leaflets having parabolic curvature at the center and linear profile at the skirt area with the parameters listed in Table 1 is constructed, as shown in Figure 6. This way, the advantages of the linear profile with a large slope at the skirts are employed, while the rest of the profile with the parabolic curvature fills the cross-section of the vessel to yield the needed geometry of the valve as close to the natural valves as possible. Because of the difficulty of running a full simulation of the models, incremental design changes are applied to the model and a fraction of the cardiac cycle is simulated to get the results. So, all the parameters are not available for the intermediate designs.

For the 3D simulations, similar boundary conditions discussed in the previous section are applied to the model. Blood is modeled as a Newtonian fluid. The large eddy simulation (LES) turbulence model is taken into account for turbulence modeling of the flow. Fluid region mesh is smoothed and moved by using a Yeoh mesh smoothing model during the FSI solution steps. The solution step is considered to be a transient solution with inertial effects activated. The fully coupled explicit solver is employed to get the results for a full cardiac cycle.

In LES, the larger three-dimensional, unsteady eddies are resolved, whereas the effect of the smaller eddies is modeled. The velocity and pressure fields are divided into resolved and unresolved scales. Denoting the fields containing all scales with capital letters, the decomposition can be expressed as follows: U = u + u ′ , P = p + p ′ . (1)

Inserting Eq. 1 into the incompressible form of the Navier–Stokes equation and the continuity equation yields the following equation: ( v , ρ ∂ u ∂ t ) Ω − ( ∇ ⋅ v , p ) Ω + ∇ v , μ ( ∇ u + ( ∇ u ) T ) Ω + ( q , ∇ ⋅ u ) Ω − ( v , f ∂ Ω = ( ∇ ⋅ v , p ′ ) Ω + ( ∇ q , u ′ ) Ω + ( ∇ v , ρ ( u ⊗ u + u ⊗ u ′ + u ′ ⊗ u + u ′ ⊗ u ′ ) ) ) Ω . (2)

In Eq. 2, f is the applied traction force on the boundary ∂Ω of the spatial domain Ω, and the stresses in the last term on the right-hand side are the resolved non-linear advection term, the two cross-stresses, and the Reynolds stress. In the residual-based variational multiscale (RBVM) method, the unresolved velocity and pressure scales are modeled in terms of the equation residuals for the resolved scales as follows: u ′ = − τ m res m ,   p ′ = − τ c r e s c (3) and the momentum and continuity equation residuals are given as follows: res m = 9 ( ∂ u ∂ t + u ⋅ ∇ u ) + ∇ p − ∇ ⋅ [ μ ( ∇ u + ( ∇ u ) T ) ] ,   res c = 9 ∇ . u (4) and the intrinsic timescales are given as follows: τ m = 1 ( C 1 9 ∇ t ) 2 + 49 2 u ⋅ G u + C 2 μ 2 G : G . τ c = C 3 ( 9 τ m trace ( G ) ) − 1 (5) Here, C₁, C₂, and C₃ are constants depending on the temporal scheme, the shape of the element, and the order of the shape functions, respectively. G is the covariant metric tensor.

The fluid–structure interaction (FSI) couplings appear on the boundaries between the fluid and the solid. The physics interface of the Comsol Multiphysics uses an ALE method to combine the fluid flow formulated using an Eulerian description and a spatial frame with solid mechanics formulated using a Lagrangian description and a material (reference) frame. The fluid flow is described by the Navier–Stokes equations which provide a solution for the velocity field u _fluid. The total force exerted on the solid boundary by the fluid is the negative of the reaction force on the fluid as follows: f = n ⋅ { − p I + ( μ ( ∇ u fluid + ( ∇ u fluid ) T ) − 2 3 μ ( ∇ ⋅ u fluid ) I ) } , (6) where p denotes pressure, μ the dynamic viscosity for the fluid, n the outward normal to the boundary, and I the identity matrix. Because the Navier–Stokes equations are solved in the spatial (deformed) frame while the structural mechanics interfaces are defined in the material (undeformed) frame, a transformation of the force is necessary. This is done according to the following equation: F = f ⋅ d v d V , (7) where dv and dV are the mesh element scale factors for the spatial frame and the material (reference) frame, respectively. The coupling in the other direction consists of the structural velocity, where dv and dV are the mesh element scale factors for the spatial frame and the material (reference) frame, respectively.

For every finite element simulation, mesh independency or convergence should be inspected. To this end, a series of simulations with various element sizes (Table 2) are performed. It is noted that since full simulations are very time-consuming because of the nature of the problem, only a fraction of the step time (0.01 s) is taken into account for mesh convergence study.

The maximum von Mises stress is obtained for the models with the element sizes listed in Table 2, and the results are shown in Figure 7. It is observed that the last two models show nearly identical results for the maximum stress on the leaflets. So, we adopt the element size equal to 0.03 mm for the full analysis.

The PC used for the simulation is a workstation with an Intel Xeon CPU with 16 cores and 64 GB of RAM. Each full cycle simulation took 12 h to complete. As mentioned earlier, the design of the leaflet is obtained by incremental changes in the geometric parameters of the leaflets and observing the effects on the overall performance of the valve. The stress distributions of the leaflets are shown in Figure 8 during a full cardiac cycle. The contour legend bar is common for all time steps with the unit of MPa for the von Mises stress. It is observed that the regions with maximum stress vary as a function of time, except for the corners of the leaflets as they experience large bending moment during the valve operation.

The maximum von Mises stress on the leaflets is equal to 0.36 MPa in the present valve design. Luraghi et al. (2017) reported the stress value of 0.85 MPa for their valve design. Note that for comparison purposes, we have used the same polymer as them in our simulations.

The flow field of the blood during the valve opening is shown in Figure 9 at different time steps. It is observed that while small vortices form from the beginning of the valve opening, considerable vortices appear as the valve approaches its full opening.

Another important parameter of the performance of the prosthetic heart valves is the opening area percentage which is obtained as follows: O A = v a l v e   c r o s s − s e c t i o n   a r e a   a t   t h e   l a r g e s t   o p e n i n g b l o o d   v e s s e l   c r o s s − s e c t i o n   a r e a × 100.

For the present design, the valve cross-section area at the largest opening is equal to 224.3 mm², which is obtained using image-digitizing software. Since the vessel radius is 10 mm, the OA is obtained to be equal to 71%. This parameter is compared with the similar published studies in Table 3. It follows that the present design yields better results, compared with the similar studies. As we discussed earlier, the slope and curvature of the leaflets play a crucial role in the opening area of the valve.

Fluid-structure interaction simulations are promising for design cycle of the prosthetic heart valves because of the detained information that are obtained by them. From the mechanical point of view, a good prosthetic heart valve should have a large opening area to allow the blood flow to be as close to the natural valve as possible, while with as small stress as possible so that the operational life of the valve is as high as possible. To come up with the proper design, first, we conducted a series of 2D simulations to assess the effects of leaflet curvature on the performance of the valve. It is observed that the leaflets with linear profile and larger slope at the skirt area give the required features. Using this outcome, a 3D model is constructed and its geometrical parameters are varied so that the best possible performance parameters are obtained. The final design is simulated using fully coupled FSI analysis with results compared with the published data. It is shown that the present model gives better results for the valve opening area and the stress developed on the leaflets.