A numerical study of hemodynamic effects on the bileaflet mechanical heart valve

The work is focused on calculating hemodynamically negative effects of a flow through bileaflet mechanical heart valves (BMHV). Open-source FOAM-extend and cfMesh libraries were used for numerical simulation, the leaflet movement was solved as a fluid-structure interaction. A real model of the Sorin Bicarbon heart valve was employed as the default geometry for the following shape improvement. The unsteady boundary conditions correspond to physiological data of a cardiac cycle. It is shown how the modification of the shape of the original valve geometry positively affected the size of backflow areas. Based on numerical results, a significant reduction of shear stress magnitude is shown. The outcome of a direct numerical simulation (DNS) of transient flow was compared with results of low-Reynolds URANS model k-ω SST. Despite the limits of the two-dimensional solution and Newtonian fluid model, the suitability of models frequently used in literature was reviewed. Use of URANS models can suppress the formation of some relevant vortex structures which may affect the BMHV's dynamics. The results of this analysis can find use in optimizing the design of the mechanical valve that would cause less damage to the blood cells and lower risk of thrombus formation.


Introduction
Heart valves regulate the direction of blood flow in heart chambers. During severely pathological states of the valves, such as aortal stenosis or mitral regurgitation, their replacement with artificial implants may be unavoidable. With the current trend of population aging, it is estimated that the current number of 280 000 heart valve transplantations per year will triple over the next 30-year period [1]. It can also be assumed that the requirements for the lifetime span of the artificial valves will further increase.
The artificial replacements can either use biological tissues or rely purely on mechanical design. Although the tissue-engineered heart valves generally reach better hemodynamic properties, their usage is limited by a relatively short lifetime. Mechanical valves have therefore a stable role in cardiovascular surgery for patients with a longer life expectancy.
The current design of mechanical heart valves suffers from several flaws. Interventions in the physiological flow result in the occurrence of high shear stress in the fluid, which can damage the fragile membranes of blood cells. With the associated effects of hemolysis or thrombogenic platelet activation, the presence of shear stress can be considered a major issue [2,3]. Other phenomena such as the formation of backflow areas or flow stagnation can also have a negative impact as they promote platelet aggregation and prolong their occurrence near the biologically foreign surface. Consequently, anticoagulant therapies are necessary, which again brings additional threat and health complications.
In the past, aortic implants were a more frequent subject of scientific research and CFD modeling than mitral replacements. The apparent reason was that aortic stenosis is a more common cause of valve replacement (53% of surgical procedures on heart valves in years 1993-2006) than mitral regurgitation (17%) [4]. Such disparity is also reflected in the current range of valve models for the mitral position [15]. In current studies, the hemodynamic effect of individual parts of both the aortic and mitral artificial valves is usually investigated using the platelet activation model on selected trajectories and the evaluation of maximum shear stress values [2,5,6]. In another work, the impact of leaflet dysfunction was investigated [7]. The influence of implantation technique on thromboembolic potential and the formation of shed vortices was shown in [8]. Overall, the CFD results of the current widely used models show the possibility of improving hemodynamic properties at spots of high shear stress occurrence and boundary layer instabilities on the inner side of the leaves [6,8].
The cardiac cycle represents a vastly complex process in terms of both experimental measurement and CFD modeling. A pulsatile flow changes from laminar to turbulent state, although a fully turbulent whirl cascade cannot be developed in such a short time [1]. Rotation of the valve leaflets and the movement of the surrounding domain responds to blood flow, muscle contractions and wall elasticity. Measurement of the three-dimensional flow field with the use of advanced laboratory techniques (such as laser Doppler anemometry and particle image velocimetry) finds its limits even in relatively simplified in vitro conditions [2]. Although CFD allows to obtain information at any scale and point of a domain, the impact of simplification at the level of the boundary conditions and numerical model remains discernible.
As the comparatively realistic CFD modeling requires quite complex numerical model as well as extensive computational resources, the introduction of some model simplifications is generally desired. Specifying the leaflet movement in order to simplify the interaction model between the fluid and solid phases was found inappropriate due to the formation of additional vortex structures [1,9]. Despite the distinctly lower computational costs, the two-dimensional solution must be used with consideration that some vortex structures can not develop due to a symmetry condition [10,11]. The significant issue of turbulence modeling on a coarser grid will be described later in this paper.
In order to handle the large deformation of a computational domain when valve leaflets move, a dynamic mesh has to be implemented. Although an arbitrary Lagrangian-Eulerian (ALE-FVM) approach is suitable for weakly coupled Fluid-Structure Interaction (FSI) problems, the mesh deformation leads to a reduction of cell quality and may result in an occurrence of invalid cells. Hence, usually some adaptive remeshing technique has to be used and tetrahedral cell zones employed. Other common methods such as Immersed Boundary Method (IBM) or Overset Grid (OG) imply non-conforming grids without the need of a remeshing process, but do not provide suitable options for variable mesh refinement according to the leaflet position (and therefore may require larger grids). Since their implementation for FSI problems with a rigid body model is not much prevalent in Open-FOAM library and its branches, the ALE-FVM approach was linked with automatic mesh generator cfMesh in this work and thus a hexahedral-dominant mesh was utilized.

Methods
Numerical simulations capture the forward flow through mitral valve during the cardiac cycle. The opening phase begins in 50ms from the start of diastolic phase and the valve remain opened until the beginning of the systole. For the valve geometry, a Bicarbon Fitline Mitral 23 model by Sorin (LivaNova) was employed. Since the valve leaflets may be considered inflexible under the given conditions, their movement was solved using Six Degrees of Freedom (SDOF) approach with weakly coupled fluid-structure interaction solver. The real hinge mechanism of the valve model has been substituted by rotational axis within a fixed point, which neglects the minor forward movement of leaflets in the opening phase ( Figure 1). Fig. 1. The hinge mechanism of the Sorin Bicarbon valve [15].
To handle the movement of the leaflets, a custom workflow algorithm was created for linking the ALE-FVM dynamic mesh approach with automated mesh generation. Generator cfMesh utilizes input geometry in STL format, which has proven to be an effective way for meshing multiple similar geometries while preserving almost identical grid parameters. Before each remeshing task, the exact position and orientation of the input geometry is updated according to the actual SDOF solution. The remeshing task is performed only when grid deformation is certain (e.g. each 500-1000 time steps). A dominant-hexahedral pseudo-2D grid is then created with the spatial density prescribed in absolute or boundary-relative manner. Values of the variable fields are transferred on the generated mesh using the mapFields utility with default cellVolumeWeight mapping method.
Even though extracted blood plasma can be considered nearly Newtonian fluid, blood itself displays as a pseudo-plastic fluid because of the presence of blood elements. Particularly at low shear rates, the tendency of erythrocytes to create formations becomes important. Within a high shear stress rates, however, such formations cannot sustain and blood tends to behave Newtonian. With regards to the nature of the flow through BMHV, a Newtonian fluid model with kinematic viscosity 4·10 -6 m 2 s -1 and density 1060kg·m -3 was employed.
The inlet boundary conditions ( Figure 2) were derived from the physiological cardiac cycle with output of 5,62l·min -1 . Reynolds number reaches around 9500 at the peak flow rate, but since the turbulence is not fully developed, the flow may be considered as transient. Estimation of the Kolmogorov length scale predicts the size of smallest eddies approximately 0,02mm. It is worth mentioning, though, that the estimate itself is based on assumption of a fully developed turbulent flow and the actual smallest vortices in the observed part of the domain may remain larger. Although many authors present the results of Direct Numerical Simulation (DNS) [5,6,9] or Large Eddy Simulation (LES), there are still frequent attempts to lower the computational costs by solving Reynoldsaveraged Navier-Stokes equations (URANS) with a corresponding turbulence model [8,[12][13][14].
DNS resolves the flow structures using Navier-Stokes equations on a very fine computational grid. Considering the smallest turbulent eddies, the grid resolution is strongly tied to the Kolmogorov length scale of the flow. Computational costs of such solution may be enormous, but the flow is captured precisely.
Using URANS approach brings certain limitations. The averaging of Navier-Stokes equations is based on the idea that instantaneous quantity can be decomposed into its time-averaged and fluctuating components. Immediate fluctuations are then represented by a turbulent stress tensor (Reynolds stress tensor) that is given by a model. Such models are based on the certain assumption of flow turbulence properties and can only provide the approximate time-averaged solution. Most of the commonly used models were originally designed to approximate a simple, fully developed and laboratorymeasurable turbulent flows [9]. Unfortunately, such characteristics do not quite correspond to the transient and vastly variable nature of the flow through heart valves.
Regardless to some negative implications, URANS models provide a huge advantage in the sense of lowering the computational costs. In the less demanding two-dimensional study, results of common low-Reynolds URANS model k-ω SST-2003 will be compared to DNS solution. In order to investigate the influence of the turbulence resolving methods, two computational grids of different density were introduced. The resolution of the computational mesh for DNS was set to 0,025mm (1,25mill cells) in the observed part of the domain and further refined near the surfaces. Despite the usual standards of the DNS approach, it was shown that the calculation of flow through BHMV on much coarser grid (with resolution in the order of Taylor microscale) provided good results when compared to experimental measurement of transvulvar pressure drop and accurate leaflet displacement [9]. To review the impact of the grid resolution, a simulation with the same settings was performed again on coarser grid of resolution 0,04mm (520k cells) that was also used with the URANS model.
Given the use of k-ω type model, boundary layer was resolved on the grid resolution that corresponded to nearwall treatment standards.
The spatial discretization was combined with second-order numerical schemes. For gradient terms, leastSquares scheme was used in order to reduce the effect of the cell non-orthogonality in some regions. While the maximum non-orthogonality values varied for different time steps and reached 34-43°, the average values ranged within 1,2-2° only. For the temporal discretization, first-order Euler scheme with a time step corresponding to a Courant number of less than 0.3 was selected.

Evaluation of RANS approach
The computational costs of remeshing tasks proved to be negligible when compared to the execution time of the simulation. On the other hand, the field mapping between deformed and freshly created grids introduced some moderate disruption of the flow field. The most affected outcome was the course of force effects on the boundaries. Since a volume-weighted mapping method was used for internal cells (with face-area weighted AMI for boundaries), the only more accurate field mapping appears to be achievable through a higher grid resolution of the boundary area.
To compare the methods of turbulence resolving, attention was paid to results with a direct impact on the hemodynamic properties of the valve. An apparent difference was seen in velocity profiles behind the valve, where the peak values differs extensively (see Figure 3). Limitation of instantaneous velocity deviations caused by Reynolds averaging also led to overall vorticity decrease in the entire domain. Since several works resolves the platelet activation along chosen trajectories [2,5,6], both the shear stress values and the trajectory shape would have been affected and the results devaluated if URANS modeling was used. The turbulence formation is also crucial in terms of disturbing platelet and thrombus deposits on the walls. In addition, more powerful instabilities are associated with increased wall shear stress occurrence. Even though the boundary layer was handled using the grid resolution (and the temporal resolution was fine enough to capture any possible vortex formation), URANS model averaged the flow near the boundary layer and no significant instabilities could develop in this area. Even though wall shear stress reached lower values with URANS modelling than for DNS simulation (as seen in Table 1), the grid resolution was much more pronounced. Vortex separation behind the both trailing edge and the inner side of leaflets causes certain oscillation of the acting moments on fully opened leaflets. By using URANS approach, the dynamics of leaflet movement was not perceptibly affected during the opening phase. In fully opened state of the valve, the acting force effects differed between DNS and URANS approach in both the oscillation frequency and character. Since URANS model did not capture the formation of eddies on the inner side of leaflets, corresponding effects on momentum oscillation were not observed. As a result, only significant waves produced by the instability behind the trailing edges remained and, moreover, their frequency has also decreased by more than 10%.

Shape modification
Formation of large vortices behind mechanical valves represents a serious change in physiological flow of the left ventricle and hence potential risk. With heavily swirled areas, there is also an increase in shear stress presence and the frequency of the blood cell collisions. The simulation of the original geometry has shown that most vortices are formed inside the middle flow channel between valve leaflets. Leading edges on the inner side of leaflets cause sudden flow separation, which breaks down to medium-sized periodic eddies ( Figure 5, top). Unfortunately, their propagation and growth seem to be further supported by already rich and even extending distance between the compact jet flow and leaflet walls. Shifting the jet separation spot further beyond the leading edge, while simultaneously reducing the distance between the jet and the wall, could therefore help to depress either size and intensity of the eddies.
Adjusting the camber of leading edges has shown a significant improvement of the jet stability. Moreover, the peak values of wall shear stress dropped noticeably from 787,2Pa with the original geometry to 521,5Pa with the final modification ( Figure 4). The shape of both the inner side of leaflets and trailing edge was changed in order to narrow the mentioned gap between jet stream and leaflet walls. Such adjustments perceptibly reduced the size of the sealing area of the valve, most manifested near the leading edges with nearly half the thickness decrease. Using new production techniques, it is likely that the mechanical durability of the leaflets would not be adversely affected; the opposite state would lead to extension of the area on the outer side of leaflet tip. In addition to the improved jet stream stability, a better uniformity of flow behind the valve was acquired and spreading whirl structures reduced ( Figure 5). The reduced size of eddies in the middle flow channel induced an increased frequency of the whirl separation. Since regions of the flow stagnation or backflow may bring an increased danger of thrombus formation, better adhesion of the current to the inner walls may be an important advantage too since it reduced the size of such areas. Despite encouraging numerical results, presented work was resolved as simplified two-dimensional flow. Some vortex structures hence could not be truly captured and the mode of their disintegration was affected as well. Even though usage of the incompressible Newtonian model of fluid was justified, it could still influence the nature of the flow (especially near the boundary layers).
The subject of the future work will be to improve the existing computational model in sense of eliminating the introduced simplification. It is important to mimic the real shape of the computational domain, which will allow flow structures to undergo a more realistic propagation. Attention should be focused on extending the hinge model by the translational motion component. That would possibly facilitate the future analysis of three-dimensional flow around the hinge mechanism and its impact on the valve hemodynamic properties.
Solution of some other problems may nevertheless remain unclear. A well-known difficulty in solving rapid valve closure using finite volume method makes it impossible to assess the permeability of the closed valve. To verify the results and answer certain questions, the calculations need to be followed by future experimental verification.

Conclusion
Present study is focused on analysis of turbulence modelling techniques and shape improvement of 23mm Sorin Bicarbon Mitral valve design. Numerical results were carried out using open-source libraries FOAMextend and cfMesh that showed good robustness and advantageous properties for optimization deployment.
Through the comparison of various turbulence modeling approaches, URANS models were found to be unsuitable for an accurate prediction of hemodynamic properties of mechanical heart valves. DNS simulations have shown similar outcome in flow behaviour when computed on both fine and coarser grids, which is in good agreement with the conclusions of the literature. However, alongside the overall good testimony of such solution, lower peak wall shear stress values as well as stronger tendency to maintain flow symmetry were registered on the coarser grid.
The outcome of the simulations showed that the valve leaflet design leaves room for an improvement and that through an adjustment of both the leading edge and inner leaflet side, a significantly better overall valve hemodynamics can be achieved. However, minor possibly negative impacts were observed in shifting the values of moments acting on the leaflets and thus the risk of their possible oscillation during the opened phase. As the present leaflet hinge model doesn't allow accurate assessment of this phenomena, it will be a subject of future work.
Faculty of mechanical engineering BUT, within the project FSI-S-17-4615, is gratefully acknowledged for support of this work.