Fully coupled dynamic simulations of bioprosthetic aortic valves based on an embedded strategy for fluid–structure interaction with contact

Abstract

Aims

This work aims at presenting a fully coupled approach for the numerical solution of contact problems between multiple elastic structures immersed in a fluid flow. The key features of the computational model are (i) a fully coupled fluid–structure interaction with contact, (ii) the use of a fibre-reinforced material for the leaflets, (iii) a stent, and (iv) a compliant aortic root.

Methods and results

The computational model takes inspiration from the immersed boundary techniques and allows the numerical simulation of the blood–tissue interaction of bioprosthetic heart valves (BHVs) as well as the contact among the leaflets. First, we present pure mechanical simulations, where blood is neglected, to assess the performance of different material properties and valve designs. Secondly, fully coupled fluid–structure interaction simulations are employed to analyse the combination of haemodynamic and mechanical characteristics. The isotropic leaflet tissue experiences high-stress values compared to the fibre-reinforced material model. Moreover, elongated leaflets show a stress concentration close to the base of the stent. We observe a fully developed flow at the systolic stage of the heartbeat. On the other hand, flow recirculation appears along the aortic wall during diastole.

Conclusion

The presented FSI approach can be used for analysing the mechanical and haemodynamic performance of a BHV. Our study suggests that stresses concentrate in the regions where leaflets are attached to the stent and in the portion of the aortic root where the BHV is placed. The results from this study may inspire new BHV designs that can provide a better stress distribution.

Fluid–structure interaction, Embedded methods, Bio-prosthetic heart valves, Haemodynamic and mechanical performance

What’s new?

• A novel approach to solve fluid–structure interaction with contact problems in bioprosthetic aortic valves.

• Methodology based on a staggered strategy which allows to combine existing state of the art methods and software.

• Estimate of the dependence of the mechanical performance on the anisotropy of the leaflets soft tissue.

• Comparison of the mechanical stresses for different valve designs.

• Assessment of haemodynamic performance of a bioprosthetic heart valve in transitional flow regimes.

Introduction

Aortic stenosis is the most prevalent valvular pathology in Western countries, and the standard medical treatment is the replacement of the native valve with a bioprosthetic heart valve (BHV) or a mechanical heart valve. BHVs are haemodynamically superior to mechanical valves, but the drawback remains their limited durability due to tissue degeneration.¹ Despite the major improvement in the valve design, further investigations are still needed towards an ideal valve replacement.

Computer modelling and numerical simulations can be used to analyse the mechanical and haemodynamic properties of a BHV. Moreover, the full performance of a BHV is directly related to the haemodynamic environment. In this regard, fluid–structure interaction (FSI) simulations can augment our understanding of the haemodynamic and the mechanical response of a BHV and a large number of embedded approaches have been proposed to simulate flow over BHVs.

Fadlun et al.² presented an embedded method based on a simple spring-network model of the deformable structures. Such a method has been recently used to compare the haemodynamics of native and bioprosthetic mitral valves.³ Ge and Sotiropoulos⁴ developed a sharp interface curvilinear immersed-boundary (IB) method, which was generalized to simulate the behaviour of human heart valves⁵ and prosthetic heart valves.^6,7,8 Kamensky et al. introduced an immersed isogeometric approach⁹ that has been recently adopted to design a computational framework for the numerical simulations of BHVs.^9,10 The parametric study¹¹ aimed at investigating the influence of the geometry on the heart valve performance and compare different hyperelastic soft tissue models.^12,13

Experimental validation of computational models of BHV is one of the most significant challenges. In this regard, Quaini et al.¹⁴ validated a simplified model of a flexible heart valve using Doppler ultrasounds. Wu et al.¹⁵ validated leaflet kinematics and the orifice area of an FSI model of a BHV against in vitro measurements. Lee et al.¹⁶ developed a computational FSI model based on the IB method of an experimental pulse duplicator platform for simulating the dynamics of BHVs. A comparison between numerical simulations and experimental data demonstrated good agreement for flow rates, pressures, and valve orifice areas.

Many studies in the literature show that numerical simulations applied to BHV design can improve their overall performance.¹⁷ However, computational models of BHVs have to account for blood–soft tissue interaction as well as the contact occurring among the leaflets during the valve closure.

Fluid–structure contact interaction approaches have been proposed in the literature for very specific scenarios and mainly based on a contact penalty-based strategy.^7,9 Indeed, the penalty method represents the easiest choice to handle contact problems and allows for the use of standard unconstrained minimization strategies.

However, over the last decades, several authors have demonstrated the superior robustness of mortar-based approaches compared to standard penalization strategies.¹⁸

In this work, we present an embedded approach based on the mortar method for the numerical solution of contact problems between elastic structures immersed in a fluid flow. This approach is designed to simulate the full dynamics of a BHV and can be used to analyse its mechanical and haemodynamic performances.

Methods

Geometry setup

A BHV consists of three leaflets made of bovine pericardium mounted on a rigid polymeric stent. Figure 1A shows a schematic representation of a BHV placed in the aortic root which is immersed in the fluid domain.

Figure 1

(A) Schematic description of the leaflets (pink), the stent (grey), the aortic root (red), and the fluid domain (blue). (B) Fibre orientation for leaflets made of bovine pericardium. (C) Staggered approach used to solve the FSI problem with contact. We solve the structural problem with contact mechanics and transfer the velocity field from the solid to the fluid domain. Then, we solve for the Navier–Stokes equations coupled with the velocity constraint and transfer the reaction force from the fluid to the solid. In the panel, we also show the geometric configuration of the model consisting of the aorta and the BHV embedded in the fluid domain. We depict the stresses during the closure of the valve on the right panel, and the velocity field when the valve is open on the left panel. BHV, bioprosthetic heart valve.

Mathematical model

We briefly state the FSI problem used to model the heart valve’s dynamics. The balance equations describing the mechanical dynamics of a BHV and the surrounding aortic root are formulated in Lagrangian coordinates:¹⁹  

$${\rho }_s\ddot{\mathit{U}}-\nabla ·\mathit{P}={\mathit{F}}_{\mathit{fsi}}.$$

(1)

Here, ${\rho }_s$ is the density, $\mathit{U}=\mathrm{ }\mathit{U}(\mathit{X},t)$ is the displacement field, $\ddot{\mathit{U}}$ is the acceleration, $\nabla ·$is the divergence operator in the reference configuration, $\mathit{P}$ is the first Piola-Kirchhoff stress tensor, and ${\mathit{F}}_{\mathit{fsi}}$ is the reaction force exerted from the fluid to the solid. The solid problem includes contact conditions¹⁸ to simulate valve closure. Leaflets, stent, and aorta are connected, and we do not need to explicitly model the aorta-stent, and stent-leaflets tied constraints.

For an hyperelastic material, the stress tensor is related to the strain tensor $\mathit{C}$¹⁹ by the relation:

$$\mathit{P}=\text{2}\mathit{F}\frac{\partial \psi \left(\mathit{C}\right)}{\partial \mathit{C}}$$

with $\psi \left(\mathit{C}\right)$ being a suitable strain–energy function and F is the deformation gradient.

As the volume of fluid in the BVH appears to be tightly bound to the solid structures, we use a nearly incompressible Neo-Hookean material to model the aortic root and the isotropic leaflets, and the St. Venant–Kirchhoff strain-energy to model the rigid stent:¹⁹  

$${\psi }_N\left(\mathit{C}\right)={\mu }_N\left(tr\left(\mathit{C}\right)-3\right)+\mathrm{ }\frac{k_{N\mathrm{ }}}{2}\mathrm{ }({\mathrm{J}-1)}^2,$$

(2)

 

$${\psi }_{SV}\left(\mathit{C}\right)={\mu }_{sv}\mathrm{ }tr\left({\mathit{C}}^2\right)+{\lambda }_{sv}\mathrm{ }[{tr\left(\mathit{C}\right)]}^2.$$

(3)

Here ${\mu }_N$ is the shear modulus, $k_N$ is the bulk modulus, ${\mu }_{sv}$ and ${\lambda }_{sv}\mathrm{ }$are the Lamé parameters, $J$ is the determinant of $\mathit{F}$⁠, and $t\mathrm{r}\mathrm{ }$is the trace operator.

In agreement with experiments¹³ leaflets are also modelled as a nearly incompressible anisotropic fibre-reinforced material:

$${\psi }_H\left(\mathit{C}\right)={\mu }_H\left(tr\left(\mathit{C}\right)-3\right)+{\sum }_{i=1}^2\frac{a_i}{2b_{\mathrm{i}}}\mathrm{ }\left(\exp \left[{b_{\mathrm{i}}\left({\mathit{g}}_i\cdot {\mathit{Cg}}_i-1\right)}^2\right]-1\right)+\frac{k_{H }}{2}\mathrm{ }({\mathrm{J}-1)}^2\mathrm{ }$$

(4)

where ${\mu }_H$⁠, $a_i,\mathrm{ }{b}_i$⁠, and $k_H$ are material parameters, while ${\mathit{g}}_i$ denotes the fiber orientation. Figure 1B depicts the two fibre families, ${\mathit{g}}_1$ and ${\mathit{g}}_2$⁠, which represent the radial and the circumferential directions, respectively.

The fluid subproblem is formulated in Eulerian coordinates.¹⁹ The Navier–Stokes equations for incompressible flow read:

$$\begin{array}{l}{\rho }_f{\dot{\mathit{v}}}_f+{\rho }_f{\mathit{v}}_f\cdot \nabla {\mathit{v}}_f-{\mu }_f\text{Δ}{\mathit{v}}_f+\nabla p_f=\mathbf{0}\\ \nabla \cdot {\mathit{v}}_f=0\end{array}$$

(5)

Here, ${\mathit{v}}_f$ is the fluid velocity vector field, ${\dot{\mathit{v}}}_f$ is its time derivative, $p_f$ is the pressure, ${\rho }_f$ is the density, ${\mu }_f$ is the viscosity, while $\mathrm{\Delta }$ is the Laplacian operator.

Our FSI approach is inspired by the embedded techniques. Thus, in our formulation the fluid domain embeds the solid structure, and the coupling between the two physics is established in the form of a vector constraint enforcing congruent velocities over the entire overlapping region represented by the solid structure:

$${\dot{\mathit{v}}}_f={\dot{\mathit{u}}}_s$$

(6)

with ${\dot{\mathit{u}}}_s$ being the solid velocity-field in the current configuration. The equality constraint is weakly enforced by means of distributed Lagrange multipliers.

For the solution of the highly non-linear FSI system described by the equations (1)–(6), we adopt a staggered approach based on the variational transfer.²⁰ In particular, at each time step, we first solve for the solid problem coupled with the fluid–structure interaction force F_fsi (i.e. the Lagrange multiplier) and the contact conditions and compute the displacement and the velocity fields of the aortic root and the BHV. Then, we transfer the velocity field from the solid structure to the fluid domain and solve for the Navier–Stokes equations (5) coupled with the constraint (6). Finally, we compute the fluid–structure interaction force and transfer it from the fluid to the solid domain. This iterative procedure is repeated until convergence is achieved, i.e. the magnitude of the difference between the force F_fsi of consecutive iterates is below a small numerical threshold (10⁻⁶). Although the contact between structures introduces non-smooth dynamics, the number of iterations required to solve the coupled nonlinear system does not increase dramatically (10–25 iterations per time step). Once, the solution strategy has converged, we update the variable values of the FSI problem and solve it for a new time step. The overall approach is depicted in Figure 1C and we use the finite-element method for discretizing both fluid and structure problems.

Data analysis

Experimental evidence shows that stress distribution and magnitude are associated with the structural deterioration of a BHV.¹ Hence, our first analysis aimed at estimating the distribution of stresses in the leaflets to identify regions with high-stress concentrations. In order to account for the tensile and compressive stresses within the leaflets, the von Mises stress was considered as the indicator of the mechanical performance of a BHV.^19,21 We performed structural simulations to analyse the mechanical performance of isotropic and fibre-reinforced leaflets. In particular, the leaflets were modelled as an isotropic Neo-Hookean material, ${\psi }_N\left(\mathit{C}\right)$⁠, with shear modulus ${\mu }_N=\mathrm{ }{1.34\mathrm{ }·10}^3\mathrm{ }\mathrm{kPa}$ and bulk modulus $k_N=\mathrm{ }{66.67\mathrm{ }·10}^3\mathrm{ }\mathrm{kPa},$²² and as a fibre-reinforced Holzapfel material, ${\psi }_H\left(\mathit{C}\right)$⁠. In the last case, all the material parameters were chosen in agreement with experimental observations.²³ We set ${\mu }_H=20.1 \mathrm{kPa}$⁠, $a_1=a_2=54.62 \mathrm{kPa},\mathrm{ }$and $b_1=b_2=30.86\mathrm{ }\mathrm{and}\mathrm{ }{k}_H=\mathrm{ }{66.67\mathrm{ }·10}^3\mathrm{ }\mathrm{kPa}$⁠. Concerning the remaining material properties, the density was set equal to ${\rho }_s=\mathrm{ }1200 \mathrm{kg}/{\mathrm{m}}^3$ for all the components. The shear modulus ${\mu }_N\mathrm{ }$of the nearly incompressible Neo-Hookean model used for the aortic root was set equal to $340 \mathrm{kPa}$⁠,²² while the bulk modulus $k_N$ was set equal to$\mathrm{ }{16.67\mathrm{ }·10}^3\mathrm{ }\mathrm{kPa}.$²² The Lamé parameters of the St. Venant–Kirchhoff model adopted for the stent were ${\mu }_{sv}=\mathrm{ }{4\mathrm{ }·10}^3\mathrm{ }\mathrm{kPa}$ and ${\lambda }_{sv}=\mathrm{ }{9.3·10}^3\mathrm{k}\mathrm{Pa},$²² respectively. A time step Δt = 0.001 s was employed for all the numerical simulations.

Valve design plays a crucial role in the distribution of the stresses and consequently in the structural deterioration of a BHV. In this work, we considered two geometries denoted by BHV 1 and BHV 2. As one may observe in Figure 2A, leaflets of the valve BHV 2 were more elongated in the axial direction than the leaflets of the model BHV 1. On the other hand, the three posts of prosthesis BHV 2 were thinner and higher than the posts of the geometry BHV 1. As an example, we depict the aggregated geometry consisting of the BHV 1 and the aortic root in Figure 2B. Here, the geometry was cut to show the positioning of the prosthesis inside the aortic root.

Figure 2

(A) Geometries for two BHV designs. The leaflets of geometry BHV 2 are more elongated in the axial direction than the leaflets of valve design BHV 1. Moreover, the geometry BHV 2 has thinner and higher posts than BHV 1. (B) Placement of BHV 1 in the aortic root and boundary conditions for the solid problem. (C) Transvalvular pressure gradient applied on the internal surface of the leaflets. (D–I) Von Mises stress distribution for different material properties and valve designs during systole (D, F, H) and diastole (E, G, I). We observe high-stress values in the central region of the isotropic leaflets in the BHV 1 model (D–E), and close to the attachment in the fibre reinforced valve BHV 2 (H– I). On the other hand, the fibre reinforced BVH 1 model (F–G) reveals the best stress distribution. All the values must be multiplied by 10² kPa. BHV, bioprosthetic heart valve.

In Figure 2B and C, we illustrate the boundary conditions used for the pure mechanical analysis. For each simulation, we applied the transvalvular pressure gradient depicted in Figure 2C on the internal surface of the leaflets and kept the bases of the stent and the aortic root fixed.

As the last contribution, we considered an FSI model of BHV to investigate the haemodynamic performances and their relations with the stress distribution on the leaflets. Blood was modelled as an incompressible Newtonian fluid with density ${\rho }_f=1000 \mathrm{kg}/{\mathrm{m}}^3$ and viscosity ${\mu }_f=0.004 \mathrm{Pa}·\mathrm{s}$⁠. The material properties of the solid components were the same as those adopted in the structural simulation. The fluid mesh consisted of 844 864 tetrahedra elements and 149 567 nodes, while the solid mesh contains 850 136 tetrahedra and 191 523 nodes. A time step Δt = 0.0005 s was adopted for the time integration schemes of the two sub-problems.

Concerning the boundary condition of the fluid problem, we imposed the velocity profile depicted in Figure 3A on the inlet of the fluid domain, and no-slip conditions on the external surface. The use of the Windkessel model on the outlet of the fluid channel, as depicted in Figure 3B, allowed to account for the peripheral artery system. In particular, we adopted a three-elements model:²⁴  

$$\begin{array}{l}Q\left(t\right)=C\frac{d{P}_c(t)}{dt}+\frac{P_c(t)}{R_d}\\ P\left(t\right)=R_{p}Q(t)+P_c(t)\end{array}$$

Figure 3

(A) Velocity profile used as an inflow boundary condition for the fluid problem. (B) Windkessel model used to get a pressure gradient between 80 and 120 mmHg. (C) Fluid flow during systole. One may observe that the jet starts to emerge from the valve and becomes fully developed when the flow rate reaches its maximum value. (D) Recirculation zones formed in the sinuses of Valsalva during diastole.

Here $R_p$ is the proximal resistance, $C$ is the compliance and $R_d$ is the distal resistance. Moreover, $Q\left(t\right)$ represents the blood flow rate, $P\left(t\right)$ is the proximal blood pressure, and $P_c\left(t\right)$ is distal pressure. By calibrating the three parameters of the Windkessel model, one can simulate the experimental physiological pressure profile varying from 80 to 120 mmHg. All the parameters were chosen in agreement with the literature.²⁴

Results

For evaluating the mechanical performances of a BHV, we first performed pure mechanical simulations and analysed the spatial distribution of the von Mises stress in the leaflets. In particular, we compared the stress state of isotropic and fibre-reinforced soft tissues. Then, we considered two different valve designs and estimated the impact of the leaflet shapes on the stress concentration and distribution.

Finally, we presented fluid–structure interaction simulations to analyse the haemodynamic performance of a BHV.

Bioprosthetic heart valve: mechanical simulations

Pure mechanical simulations were performed to find the valve configuration which predicts the best stress distribution. In Figure 2, the von Mises stress patterns are represented for all the simulated scenarios. The stress states associated with the systolic and the diastolic phase are depicted in Figure 2D and E for the isotropic case, and in Figure 2F–I for the fibre-reinforced material. In particular, Figure 2F and G refer to the valve design BHV 1 while the results reported in Figure 2H and I depict the stress distribution in the BHV 2 model.

We observed a stress concentration in the middle zone of the leaflets for the isotropic BHV 1 model and along the junction between leaflets and stent in the BHV 2 fibre-reinforced geometry. On the other hand, the BHV 1 fibre-reinforced model minimizes the stress values compared to the other two configurations presented in this study.

Bioprosthetic heart valve: fluid–structure interaction with contact simulations

In this section, we present the results obtained by performing a fluid–structure interaction simulation with contact mechanics. As depicted in Figure 1A, we immersed the geometry BHV 1 and the surrounding aortic root into a fluid channel. Results were analysed in terms of stresses and fluid flow.

Figure 3C shows the snapshots of the flow field at the selected time instants of the same heartbeat. During the systolic acceleration, we observed a fully developed jet with a maximum velocity of about $2.0\hspace{0.17em}\mathrm{m/s}$ confined to the core region of the aortic root. Figure 3D depicts a snapshot of the diastolic blood deceleration. We observed flow recirculation and fluctuations along the sinus and commissure sides of the aorta.

Figure 4 reports snapshots of the spatial distribution of the stresses in the leaflets (A–F) and the stent (G–N). The stress state is associated with six successive BHV configurations. We show the stresses undergone by leaflets and stent on the top and the bottom, respectively. We observed high-stress values in the central region of the leaflets during systole. On the other hand, stress concentrated on the coaptation of the leaflets during the closure of the valve. High-stress values were also computed in the stent posts during the systolic phase.

Figure 4

Snapshots of the von Mises stress in the leaflets (A–F) and the stent (G–N). A stress concentration is observed in the central region of the leaflets when the valve is fully open, and on the coaptation of the leaflets during the valve closure. High-stress values are found in the three stent posts during systole. All the values must be multiplied by 10² kPa. Stress values are reported in the legend on top of figure (B) for the leaflets, and on the bottom of figure (M) for the stent.

For a complete analysis of the stress state, Figure 5 depicts the spatial distribution of the von Mises stress on the external and internal surface of the aorta at systole and diastole. We observed a stress concentration along with the sinuses of Valsalva where the stent is attached to the aortic root.

Figure 5

Spatial distribution of the von Mises stress in the aortic root during systole (top) and diastole (bottom). We observe a stress concentration in the region where BHV is attached to the aorta during the entire heartbeat. All the values must be multiplied by 10² kPa. BHV, bioprosthetic heart valve.

Discussions

Clinical failures of BHV are generally related to structural deterioration and calcification. One of the most common hypotheses is that stress concentration and abnormal haemodynamic may promote the progressive deterioration of the leaflets soft tissue.^1,26 The aim of this work is to present a framework for analysing the performance of a BHV, identifying regions affected by high-stress values, and relating them to the BHV design and leaflets material model. Our numerical results suggest that stresses are substantially reduced if leaflets are modelled as fibre-reinforced material, and that the valve design may influence the stress patterns. In this regard, we observed that elongated leaflets experience higher stresses along with the attachment to the stent. The use of an FSI model allows the assessment of the mechanical and haemodynamic characteristics of a BHV. Our analysis shows a stress concentration in the central region of the leaflets in the systolic phase and close to the attachment between leaflets and stent during the valve closure. All our numerical results agree very well with previous studies^{11,12,21,25,27} where the authors observed a stress/strain concentration in the central region of the leaflets during systole and close to the stent attachment during diastole.

The overall stress state analysis shows that the leaflet-stent attachment, the three stent posts, and the portion of the aortic root attached to the stent are the most solicited regions, and careful attention should be paid to them at the stage of valve design and placement. Moreover, we point out that the compliance of the aortic root reduces the level of the stresses and strain occurring at the point of the valve closure, as already discussed in a previous study.²⁸

Haemodynamic characteristics of a BHV are also related to their clinical performance. In this regard, the fully coupled approach adopted in our computational framework allows the estimation of abnormal flow patterns and recirculating areas.

In our study, we analysed the haemodynamic performance of the valve design BHV 1 which revealed an optimal stress distribution compared to the other configurations. The fully developed forward flow observed in the opening phase, as well as the flow recirculation zones occurring in diastole at the sinus side of the aorta, agree with experimental observations. ²⁵ Moreover, as expected, the physiological conditions do not introduce stagnant flow regions persisting throughout the whole cardiac cycle. This represents an important feature of the BHV design since the persistence of vortices or recirculating zones may be responsible for the early deterioration in a BHV.^1,26

Conclusions

In this work, we show that our FSI approach can be employed to identify regions of high-stress concentration and flow reversal zones. We point out the optimization of the leaflets stress distribution and its relation with abnormal flow patterns are extremely useful for improving the clinical performance of a BHV. Future work could account for the soft-tissue prestress, the non-Newtonian rheological characteristics of blood flow, and mechanical damage under cyclic loading. Moreover, validation against experimental data represents an essential step towards the full integration of advanced simulations into the clinical practice and medical device manufacturing. Indeed, FSI numerical simulations of BHV could support biomedical engineers for optimizing the valve design and clinicians for improving prosthesis patient matching.

Acknowledgements

The authors thank Barna Becsek (ARTORG Center for Biomedical Engineering Research, https://www.artorg.unibe.ch) for providing the BHV geometries.

Funding

This work was financially supported by the Theo Rossi di Montelera Foundation, the Metis Foundation Sergio Mantegazza, the Fidinam Foundation, the Swiss Heart Foundation, and the PASC projects FASTER (Forecasting and Assessing Seismicity and Thermal Evolution in geothermal Reservoirs) and HPC-Predict (High-Performance Computing for the Prognosis of Adverse Aortic Events). This paper is part of a supplement supported by an unrestricted grant from the Theo-Rossi di Montelera (TRM) foundation.

Conflict of interest: none declared.

Data availability

No clinical data has been used in the present study.