Overset meshing in combination with novel blended weak-strong fluid-structure interactions for simulations of a translating valve in series with a second valve

Abstract

Mechanical circulatory support (MCS) devices can bridge the gap to transplant whilst awaiting a viable donor heart. The Realheart Total Artificial Heart is a novel positive-displacement MCS that generates pulsatile flow via bileaflet mechanical valves. This study developed a combined computational fluid dynamics and fluid-structure interaction (FSI) methodology for simulating positive displacement bileaflet valves. Overset meshing discretised the fluid domain, and a blended weak-strong coupling FSI algorithm was combined with variable time-stepping. Four operating conditions of relevant stroke lengths and rates were assessed. The results demonstrated this modelling strategy is stable and efficient for modelling positive-displacement artificial hearts.

Keywords:

• Computational fluid dynamics
• fluid-structure interactions
• mechanical heart valve
• positive displacement
• total artificial heart

1. Introduction

Heart failure (HF) currently affects 64.35 million patients worldwide (Lippi and Sanchis-Gomar 2020). Treatments for HF range from healthy lifestyle changes and medication to heart transplant surgery, which is considered the gold standard in treating end-stage HF. However, the number of available donors is limited. Since 2011, the number of patients on the heart transplant waiting list in the UK has increased by 85%, but there are only half the number of donors available (NHS BaT 2021). To alleviate the shortfall, doctors can turn to mechanical circulatory support devices to bridge the gap to transplant. These devices include ventricular assist devices, used in cases of ventricular failure to partially assist in cardiac function, and total artificial hearts (TAHs), which fully replace the native heart in cases of biventricular HF.

The Syncardia (Tucson, AZ, USA) TAH, a pneumatically driven pulsatile TAH (Slepian et al. 2013), is currently the only FDA-approved bridge-to-transplant TAH (Melton et al. 2019). However, it suffers from complications including infection and thromboembolic events (Torregrossa et al. 2014) that newer TAHs are seeking to overcome. The CARMAT (Paris, France) Aeson TAH, a positive-displacement device that aims to improve biocompatibility by using bioprosthetic materials and mimicking the shape of the heart (Mohacsi and Leprince 2014), has recently been granted the CE mark in Europe and has gained approval for early US feasibility studies (Han 2021).

The Realheart TAH, developed by Scandinavian Real Heart (Västerås, Sweden), is a novel TAH that aims to improve biocompatibility by mimicking a key component of the mechanics of the native heart: the vertical translation of an atrioventricular (AV) plane (Szabo et al. 2018). When the native heart contracts during systole, the vertical displacement of the base of the ventricle towards the apex accounts for 60% of the stroke volume (Carlsson et al. 2007). The Realheart TAH employs this positive-displacement pumping mechanism to generate pulsatile blood flow, using a pair of bileaflet mechanical heart valves (BMHVs) in the atrioventricular and semilunar positions to govern the forward motion of blood through the device. As such the BMHVs are a key component of the TAH and understanding their motion is vital to predicting fluid flow and pressure characteristics. Mechanical heart valves (MHVs) have been associated with complications such as haemolysis and coagulation due to non-physiological flow patterns that are generated (Dasi et al. 2009) and therefore knowledge of shear stresses and exposure times around the valves are also required.

Computational fluid dynamics (CFD) has been widely used to simulate blood flow through MHVs to gain an insight into these problems. Simulating the interaction between MHVs and the blood is a fluid-structure interaction (FSI) problem and can be solved using different numerical techniques including: the arbitrary Lagrangian-Eulerian (ALE) or moving grid method, the immersed boundary (IB) or fixed grid method, and the overlapping domain (OD) method, known as overset or chimera meshing. The ALE method is the most common FSI technique, as it is widely available and in theory more accurate for near wall modelling when compared to the IB method, as the ALE method explicitly captures the fluid-solid interface (Bavo et al. 2016). However, recent high-order IB methods that implement sharp interfaces between the fluid and solid result in improved shear stress capture (Kolahdouz et al. 2021). IB methods have been developed with large displacements and elastodynamic behaviour of structures in mind (Tian et al. 2014; Nestola et al. 2019) and have been applied to bioprosthetic and trileaflet prosthetic heart valves (Hsu et al. 2015, Borazjani 2013, Gilmanov et al. 2015) as well as bicuspid aortic valves (Gilmanov and Sotiropoulos 2016). Early approaches of the ALE method simulated 2D leaflet sections (Dumont et al. 2004), but as computational power improved, the models were extended to full 3D MHVs in realistic aortic geometries (Nobili et al. 2008; Mirkhani et al. 2016). The IB method can handle large solid displacements efficiently for MHV FSI problems (Borazjani et al. 2010), as well as providing improved stability over the ALE method as mesh motion and deformation are not required (Bavo et al. 2016). The OD method features less frequently for MHV applications, possibly due to a lack of commercial software support until recently. It offers accuracy, like the ALE method, without the need for mesh deformation, like the IB method (Spühler and Hoffman 2021). Compared to the ALE approach, the OD method also simplifies mesh generation of more complex devices such as a TAH, where each component can be individually meshed and combined. This approach also allows for new components to be included or other parts suppressed with ease, whilst also allowing for the swapping of parts to analyse geometry variations, all without the need for a complete regeneration of the mesh. The OD approach has been used to obtain flow structures in a fully open valve configurations (Ge et al. 2003), as well as prescribed motion valves in an anatomically accurate human aorta (Nowak et al. 2020). Nowak et al. (2020), concluded the OD method was stable and efficient, and produced typical BMHV fluid flow phenomena.

FSI methods have also been applied to full device simulations. Luraghi et al. (2018) investigated the CARMAT TAH using a combined IB and ALE approach. The initial IB approach determined valve and membrane displacement, which was then used as prescribed motion in a subsequent ALE study to determine local velocity fields and shear stresses. Al-Azawy et al. (2016) used the OD method to model an LVAD with prescribed motion of the monoleaflet valves and pusher plate, and returned good agreement with experimental data.

For the Realheart TAH, Kelly et al. (2022) investigated four computational methods for the AV plane and valve motion using an IB approach. These were: prescribed idealised motion; prescribed motion using motion data obtained through video capture of the valves in vitro; prescribed motion of the mitral valve but fluid-driven motion of the aortic valve; and fluid-driven motion of both valves. The motion prescribed from video analysis of experiments returned the most accurate solution, whilst the fluid-driven approach could not account for the instabilities of the complex valve motion. However, the video data approach was restrictive as the range of in silico studies was limited by the experimental analysis available. Additionally, the IB approach resulted in poor shear stress resolution on the leaflet surface – an area of great interest with regards to blood damage.

This study sought to produce a new modelling approach for the motion of positive-displacement MHVs, by creating a fluid-driven strategy which would eliminate the need for prior experimental data and provide accurate shear stress capture in the vicinity of the valve leaflets. The modelling strategy was developed and tested on a simplified positive-displacement pump geometry and the valve kinematics, fluid flow characteristics and shear stresses were analysed to assess model stability and accuracy. The modelling strategy can be readily transferred to the Realheart TAH, as well as other positive-displacement artificial hearts with MHVs.

2. Methods

2.1. Computational domain

To accelerate development of the modelling approach, a simplified positive displacement pump embodying the critical features of the Realheart TAH was used. The computational domain consisted of a 3D fluid cylinder containing two BMHVs based on the On-X valve, aligned in series, shown in Figure 1(a), and was generated using Ansys DesignModeler V2021.R1 (ANSYS Inc., Canonsburg, PA, USA). The model had two degrees of symmetry, in the yz and yx planes, and as such a quarter model was considered, shown in Figure 1(b).

Figure 1. Outline of the computational domain (a) Full domain geometry with fluid cylinder and two BMHVs (b) 1/4 domain geometry: yz-symmetry (red), yx-symmetry (green) (c) Angle of leaflet in closed and open positions (d) Monitor point location with inlet, outlet, and wall boundaries.

The two valves represented the mitral (upstream) and aortic (downstream) valves. The valve leaflets were positioned at 40^° to the horizontal when closed and 84^° to the horizontal when fully open, shown in Figure 1(c), similar to Mirkhani et al. (2016). The motion behaviour of the valves, which was the same as the full Realheart TAH, was as follows: the mitral valve underwent a forced vertical sinusoidal translation, approximating the translation of the AV plane. As the mitral valve travelled upwards, a pressure gradient across the valve leaflets caused the mitral valve to open. When the mitral valve travelled downwards, a pressure gradient was created in the opposite direction and the mitral valve leaflets closed. This then increased the pressure between the valves and opened the aortic valve. The valve leaflets of both the aortic and mitral valves could open and close freely due to the fluid-structure interaction, whilst the sinusoidal translation of the mitral valve was prescribed at the pivot point.

Two key simplifications were made that improved computational efficiency – the valve housing was not explicitly modelled, and the peripheral gap around the outside of the leaflet, as well as between the leaflets, was not considered. Instead, the cross section of the fluid cylinder was the same shape as the valve leaflets in the fully closed position. This approximated a perfect seal between the valve leaflets and fluid cylinder when the valves were fully closed. This removed the need for regions of high cell density around the leaflet edges to explicitly model the gaps, reducing the overall cell count and improving model efficiency.

A monitor point was placed downstream of the aortic valve, representing the aorta, to measure pressure variation, shown in Figure 1(d).

2.2. Boundary conditions

The left side of the Realheart TAH was used as a reference for boundary condition values. At the outlet, a two-element Windkessel model approximated the haemodynamic pressure response of the downstream vasculature: resistance $R=$3e8 kg m⁻⁴ s⁻¹ and compliance $C=$1e-10$$ m⁴ s² kg⁻¹. The value of $R$ was chosen to generate a peak systolic pressure of 120 mmHg and $C$ was then varied to produce a smoothed pressure response. The inlet was assigned a constant pressure of 15 mmHg approximating pulmonary venous pressure (Kelly et al. 2022). Symmetry conditions were applied to the two symmetry faces, and no-slip walls were used on the remaining outer surfaces of the fluid cylinder, as well as the surfaces of the leaflets.

2.3. Meshing

Overset meshing was used to incorporate the BMHVs into the fluid cylinder. The overset meshing approached maintained a consistent and refined mesh around the surface of the valve leaflets, enabling accurate shear stress capture. Five fluid mesh zones were generated: the background zone for the fluid cylinder, two component zones for the mitral and aortic valves, and two component collar meshes, shown in Figure 2(a). Due to some differences in mesh resolution between the background mesh and leaflet meshes, collar meshes were included to improve the resolution of the contact between the outer faces of the valve leaflets and the background walls and symmetry faces.

Figure 2. Fluid mesh definition (a) Fluid mesh zones: mitral component zone (blue), aortic component zone (green), background zone (grey), collar zones (red) (b) before and after initialisation of whole mesh (c) before and after initialisation of mitral mesh zone.

The background zone mesh was split into coarse and fine element regions of 2 mm and 0.75 mm element size respectively and joined together with conformal interfaces. The fine regions overlapped with the valve component zones, which had a body element sizing of 0.75 mm and face sizing of 0.35 mm on all valve leaflet surfaces. Prism layers were used on all walls and the leaflet surfaces to capture boundary layers. Upon initialisation of the solution, the separate mesh zones combined to form a single continuous mesh, shown in Figure 2(b) and (c).

Initial meshes were generated using Ansys Meshing V2021.R1 (Ansys Inc.). The collar zones were meshed with hexahedral elements. The background zone and valve component zones were meshed with tetrahedral elements, which were later converted to polyhedral elements within Ansys Fluent V2021.R1 (Ansys Inc.), reducing the total cell count by a factor of approximately 3x.

The mesh convergence study used three different meshes of increasing element counts, shown in Table 1. Each mesh was solved for two cycles, and the valve kinematics and wall shear stresses were monitored. The focal point of mesh refinement was in the immediate region around the valve leaflet, where the leaflet face element size, valve component zone body size, the fine mesh body size of the background zone and collar component zone body size were varied to produce the three meshes.

Table 1. Mesh densities of the three meshes used in the mesh convergence study.

	Mesh Density	Tetrahedral Element Count	Polyhedral Element Count
1	Coarse	0.84M	0.28M
2	Medium	1.2M	0.4M
3	Fine	2.1M	0.7M

2.4. Valve motion

Motion of the valves was achieved through the dynamic mesh utility and intrinsic six-degrees of freedom (6DoF) solver within Ansys Fluent. Pre-existing User Defined Functions (UDFs) were adapted and used to assign the motion properties to each of the moving zones.

The Define_CG_Motion UDF was used for the vertical translation of the mitral valve pivot, denoted as the AV plane, by assigning a sinusoidal $y$ velocity to the mitral component zone (Figure 2(a)) with Equation (1)(1) $$v_y=-\mathit{\text{Afcos}}\left(\mathit{\text{ft}}\right)$$(1) , (1) $$v_y=-\mathit{\text{Afcos}}\left(\mathit{\text{ft}}\right)$$(1) where $A$ (mm) was the amplitude, or half the stroke length of the translation, $f$ (rad s⁻¹) was the frequency or stroke rate (bpm) of the translation, and $t$ (s) was the time. Four studies were undertaken at $A$ = 12.5, 25 mm (equivalent to a stroke length of 25 mm and 50 mm) and $f$ = 80, 100 bpm, and the response of the model to these parameters was measured.

The Define_SDOF_Properties UDF was employed to achieve fluid-driven rotation of the valve leaflets. For the valve leaflets, the motion was constrained to rotate around the $x$ axis at the valve pivot point. Contact between leaflets was not explicitly modelled, instead the motion was constrained by including limits on the angle of rotation, $\theta$ (^o), between 40° (fully closed) and 84° (fully opened). The rotation of the valve leaflets was calculated by the 6DoF solver. The mass of each leaflet was $m=$0.35 g, and the moment of inertia $I=$2.3e-8 kg m².

2.5. Numerical procedure

Ansys Fluent was used to solve the Navier-Stokes and 6DoF equations. The solution strategy used an angle-dependent blend of weak and strong coupling algorithms between the fluid flow and structural domain (Figure 3). Strong coupling, where the position of the valve leaflets was updated at the end of every iteration, was used to maximise stability and accuracy of the solution. Weak coupling, where the position of the valve leaflets was updated at the end of every time step, was used to maximise computational efficiency. With this blended approach, it was possible to choose which phases of valve motion were strongly and weakly coupled depending on threshold valve leaflet angles, thus achieving a balance between stability, accuracy and efficiency. In this study, the focus was minimising computation time, whilst still maintaining a minimum level of solution stability. This was achieved by only using strong coupling when one valve closed and the other opened, activating when valve angles were less than 7^° (scheme denoted as B1). Additionally, a scheme which utilised strong coupling for both small (<7^°) and large (37^°- 44^°) valve angles was also considered (denoted as B2), the larger threshold related to leaflet rebound and flutter, and both schemes were compared to a fully strongly (denoted as FS) coupled solution at a stroke rate of 80 bpm and stroke length of 25 mm.

Figure 3. Solution strategy flow diagram. Rectangular boxes represent processes, diamond boxes represent decisions, parallelograms represent input/output. Arrows are the flow of logic through the strategy. n is the time step number, k is the iteration number, and $T_{\mathit{\text{END}}}$ is the end time of the solution.

No turbulence model was used, as the peak bulk flow Reynolds number (Re) was 900, well below the threshold for turbulence in a pipe (Darbyshire and Mullin 1995).

The coupled algorithm was used to solve the Navier-Stokes equations, using a second-order pressure discretisation scheme and a second-order upwind scheme for the momentum. A least-squares cell based gradient scheme was used throughout. A first order pressure interpolation scheme using a least-squares method was used at overset interfaces. Solution relaxation was used for both the pressure and momentum variables, as well as for higher order terms, to improve convergence behaviour and reduce the risk of solution divergence. The default value of 0.75 was used for these relaxation parameters.

The model was initialised with the mitral valve at $y=$0 mm and moving downwards with a velocity $v_y=-\mathit{\text{Af}}$ m s⁻¹ (eq. (1)(1) $$v_y=-\mathit{\text{Afcos}}\left(\mathit{\text{ft}}\right)$$(1) with $t=$0 s). The aortic valve was in a fully open state so that there was no large displacement of the aortic valve at the beginning of the solution. This allowed the fluid flow solution to converge faster during this time, reducing computational cost. The fluid was initialised in a quiescent state, with a uniform pressure of 90 mmHg. The outlet Windkessel model was also initialised at 90 mmHg so that there was no initial pressure jump across the outlet interface. Fluid properties were density $\rho =$1060 kg m⁻³, and fixed viscosity $\mu =$0.0035 N m⁻¹ s⁻², approximating blood as a Newtonian fluid.

A first-order implicit transient formulation scheme was used, as the second-order scheme was not compatible with overset meshing. A custom angle-dependent variable time stepping scheme was developed and implemented to improve accuracy and stability where needed, whilst improving computational efficiency. The scheme was divided into two main conditions, outlined in Table 2. A maximum $\mathit{\text{dt}}$ of 2.5 ms was used as solution instabilities were experienced with larger values.

Table 2. Variable time stepping scheme conditions. Two conditions depend on moving and stationary valves. When valves are moving, two sub conditions interrogate the leaflet angle, and adjust the time step accordingly.

Number	Condition		Timestep Size
1 (Moving)	$\left|{\omega }_{\mathit{\text{mitral}}}\right|>0$ or $\left|{\omega }_{\mathit{\text{aortic}}}\right|>0$	${\theta }_{\mathit{\text{mitral}}}>2°$ or ${\theta }_{\mathit{\text{aortic}}}>2°$	$\mathit{\text{dt}}=\mathrm{\text{min}}\hspace{0.17em}\left(\frac{0.5°}{{\omega }_{\mathit{\text{mitral}}}},\frac{0.5°}{{\omega }_{\mathit{\text{aortic}}}},2.5\mathit{\text{ ms }}\right)$
		${\theta }_{\mathit{\text{mitral}}}<2°$ and ${\theta }_{\mathit{\text{aortic}}}<2°$	$\mathit{\text{dt}}=\mathrm{\text{min}}(\frac{0.1°}{{\omega }_{\mathit{\text{mitral}}}},\frac{0.1°}{{\omega }_{\mathit{\text{aortic}}}})$
2 (Stationary)	${\omega }_{\mathit{\text{mitral}}}=0$ and ${\omega }_{\mathit{\text{aortic}}}=0$		$\mathit{\text{dt}}=2.5 \mathit{\text{ms}}$

All simulations were performed in parallel on an Intel Xeon workstation with 4 cores at 3.50 GHz and 128 GB of memory.

3. Results

All results shown use the small angle blended scheme (B1), except for the 'Blended Weak Strong Coupling’ section which compares two blended schemes with a fully strongly coupled scheme.

3.1. Mesh study

The three meshes were compared using the leaflet angle and area average wall shear stress ($\overline{\mathit{\text{WSS}}}$) of the mitral and aortic valve leaflets over the course of two cycles at a stroke length of 25 mm and stroke rate of 100 bpm.

The mitral and aortic valves in the fine and medium meshes opened and closed at similar times, but the valves in the coarse mesh opened at slightly different times (Figure 4). The profile of the mitral and aortic valve $\overline{\mathit{\text{WSS}}}$ was similar for both the fine and medium meshes, whilst the coarse mesh $\overline{\mathit{\text{WSS}}}$ for both valves reflected the differences in the opening and closing times (Figure 5). Peak mitral valve $\overline{\mathit{\text{WSS}}}$ decreased for an increase in mesh density, whilst peak aortic valve $\overline{\mathit{\text{WSS}}}$ increased for an increase in mesh density. The root mean square error (RMSE) between the medium and fine mesh was 2.0^° (4.6% of the peak vale) and 2.9^° (6.6% of the peak value) for the mitral and aortic valve leaflet angles respectively, and 4.3 Pa (9% of the peak value) and 8.5 Pa (6.8% of the peak value) for the mitral and aortic valve $\overline{\mathit{\text{WSS}}}$ respectively. In all cases the coarse mesh returned greater RMSE values. The total CPU hours for the medium and fine meshes was 192 and 360 respectively. Therefore, the medium mesh was used for the study as it balanced computational time and accuracy.

Figure 4. Variation in $\theta$ over time for the coarse, medium, and fine meshes for the (a) mitral and (b) aortic valves.

Figure 5. Variation in $\overline{\mathit{\text{WSS}}}$ over time for the coarse, medium and fine meshes for the (a) mitral and (b) aortic valves.

3.2. Variable time stepping

The change in ${\theta }_{\mathit{\text{mitral}}}$ and ${\theta }_{\mathit{\text{aortic}}}$ varied linearly with $\mathit{\text{dt}}$ size, limited to $0.5°/\mathit{\text{dt}}$ above $2°$ and $0.1°/\mathit{\text{dt}}$ below $2°$ (Figure 6). The maximum $\mathit{\text{dt}}$ of 2.5 ms occurred when both valves were not moving, or when the change in angle per time step for either valve was less than $0.5°$/$\mathit{\text{dt}},$ as seen between the 300^th and 350^th time step. Compared to a fixed $\mathit{\text{dt}}$ scheme of 50,000-time steps for four cycles, the variable $\mathit{\text{dt}}$ scheme required 2,100-time steps to complete the four cycles, a reduction of 96%.

Figure 6. Variability in time step size and $\theta$ for the mitral and aortic valves, showing the effect of the variable time stepping scheme outlined in Table 2.

3.3. Flow distribution

Time dependant velocity fields for the four cases tested are shown via videos in the supplementary files. The distribution of the fluid velocity was investigated over the last cycle, at time points corresponding to both valves fully open ($\theta =44°$) (Figure 7).

Figure 7. Velocity magnitude and pressure contours when leaflet angle is $\theta =44°$ (fully open) for a change in stroke length and stroke rate. Note the scaling of the velocity magnitude contours is 0-0.3 m s⁻¹ for the 25 mm stroke length cases, and 0-0.6 m s⁻¹ for the 50 mm stroke length cases. Velocity magnitude contour plots are clipped to the maximum values of 0.3 m s⁻¹ and 0.6 m s⁻¹ in each of the cases respectively.

3.4. Leaflet kinematics

${\theta }_{\mathit{\text{mitral}}}$ and ${\theta }_{\mathit{\text{aortic}}}$ was captured for each operating condition and is shown alongside the sinusoidal displacement of the AV plane (Figure 8). The opening of one valve and the closing of another, and vice versa, occurred at the same time step. Mitral valve closure and aortic valve opening happened when the AV plane was translating downwards, whilst mitral valve opening and aortic valve closing occurred when the AV plane was translating upwards.

Figure 8. Valve motion over final cycle alongside AV plane displacement for the (a) mitral valve and (b) aortic valve.

There was an increase in ${\omega }_{\mathit{\text{opening}}}$ for both the mitral and aortic valve for an increase in stroke rate and stroke length (Figure 9). However, ${\omega }_{\mathit{\text{closing}}}$ only increased for both valves with an increase in stroke length but did not increase with an increase in stroke rate.

Figure 9. Change in (a) ${\omega }_{\mathit{\text{closing}}}$ and (b) ${\omega }_{\mathit{\text{opening}}}$ for the mitral and aortic valve leaflets with a change in stroke length and stroke rate.

3.5. Pressure variation

The variation in pressure in the aortic section was measured over time to assess the effect of the two-element Windkessel at the outlet (Figure 10). When the aortic valve was closed (before approximately cycle number = 1.9), the aortic pressure was a constant 90 mmHg for all cases. When the aortic valve opened, the pressure in the aorta increased as the flow rate through the outlet also increased, where the peak aortic pressure occurred at the peak outlet flow rate. The peak aortic pressure increased with an increase in stroke length and stroke rate; stroke length had the largest influence. When the flow began reversing, the pressure in the aorta dropped below the outlet operating pressure of 90 mmHg, until the aortic pressure reached its minimum as the aortic valve closed. The minimum aortic pressure at this time decreased with an increase in stroke length, however there was little difference when varying stroke rate. After aortic valve closure, the aortic pressure in all four cases returned to the operating pressure due to the compliance element of the Windkessel.

Figure 10. Variation in aortic pressure with a change in stroke length and stroke rate. Time frame limited between cycle number 1.75 and 2.75 to highlight variation during full range of aortic valve motion. Vertical lines correspond to instantaneous pressure spikes.

3.6. Shear stress

Peaks in both ${\overline{\mathit{\text{WSS}}}}_{\mathit{\text{mitral}}}$ and ${\overline{\mathit{\text{WSS}}}}_{\mathit{\text{aortic}}}$ occurred at instantaneous valve opening and closing; the magnitude of the peak was generally greater for valve opening (Figure 11). Time-averaged $\overline{\mathit{\text{WSS}}}$ ($\overline{\mathit{\text{TAWSS}}}$) was calculated for each valve leaflet and is shown in Table 3, where $\overline{\mathit{\text{TAWSS}}}$ increased for both valves for an increase in stroke rate and stroke length. The distribution of wall shear stress over the valve leaflet surfaces when the valves were fully open shows the largest values are along the leading edges (Figure 12).

Figure 11. Variation in $\overline{\mathit{\text{WSS}}}$ over the course of the last cycle for the (a) Mitral valve and (b) Aortic Valve.

Figure 12. Contour plot of wall shear values on surface of (a-d) aortic valve leaflet and (e-h) mitral valve leaflet when leaflet angle is $\mathrm{\theta }=44°$ (fully open), for a change in stroke length and stroke rate. Note a difference in scaling between individual contour images.

Table 3. $\overline{\mathit{\text{TAWSS}}}$ values for a variation in stroke length and stroke rate.

	25 mm/ 80 bpm	25 mm/ 100 bpm	50 mm/ 80 bpm	50 mm/ 100 bpm
${\overline{\mathit{\text{TAWS}}}}_{\mathit{\text{mitral}}}$	0.80	1.08	1.89	2.31
${\overline{\mathit{\text{TAWS}}}}_{\mathit{\text{aortic}}}$	1.16	1.50	1.98	2.32

3.7. Volume flow rate

The volume flow rate, $Q$ (L min⁻¹), was measured through the inlet, $Q_{\mathit{\text{inlet}}}$ (Figure 13(a)). Forward flow is indicated by positive volume flow rate and backflow is negative. An increase in stroke rate and stroke length both returned an increase in peak $Q_{\mathit{inlet}}$. The cumulative output fluid volume, $V$ (mL), defined as the fluid moved forward through the domain less the backward fluid flow, was calculated (Figure 13(b)). $V$ increased for an increase in stroke length, however after four cycles, there was little difference in $V$ with an increase in stroke rate.

Figure 13. (a) Variation in $Q_{\mathit{\text{inlet}}}$ for all four cycles with a variation in stroke length and stroke rate (b) Cumulative outlet flow over the four cycles.

Comparing the forward and backward fluid volumes ($V_{\mathit{\text{forward}}}$ and $V_{\mathit{\text{backward}}}$) in Figure 14(a), $V_{\mathit{\text{forward}}}$ increased for an increase in stroke length but decreased for an increase in stroke rate. $V_{\mathit{\text{backward}}}$ decreased with both an increase in stroke rate and stroke length. Peak volume flow rate, $Q_{\mathit{\text{peak}}}$ (L min⁻¹), increased with an increase in stroke rate and stroke length (Figure 14(b)). The change in average flow rate, $Q_{\mathit{\text{averge}}}$ (L min⁻¹), was marginal for a change in stroke rate, but a doubling in stroke length increased $Q_{\mathit{\text{average}}}$ by a factor of approximately three.

Figure 14. Variation in (a) volume flow and (b) volume flow rate through the outlet for a change in stroke length and stroke rate.

3.8. Blended weak strong coupling

A fully strongly (FS) coupled solve was compared to the two blended schemes to assess any loss of solution accuracy when using a blended approach. The first blended scheme (B1) considered small (<7^°) valve angles, whilst the second scheme (B2) considered both small (<7^°) and large (37^°- 44^°) valve angles.

Largest differences in mitral valve angle were observed during valve closure, where B1 closed 0.016s before B2, and 0.024s before FS (Figure 15(a)). The aortic valve reached fully open in B1 first, 0.032s before B2 and 0.061s before FS (Figure 15(b)). Timings of peak $\overline{\mathit{\text{WSS}}}$ occurred at the same time as valve opening and closing of the three schemes. FS returned a smoother $\overline{\mathit{\text{WSS}}}$ profile for both valves during valve opening, with both B1 and B2 displaying small oscillations (Figure 15(c) and (d)).

Figure 15. Comparison of blended scheme 1 (B1), blended scheme 2 (B2) and fully strongly coupled (FS) schemes for (a) mitral valve leaflet angle, (b) aortic valve leaflet angle, (c) mitral valve $\overline{\mathit{\text{WSS}}}$ and (d) aortic valve $\overline{\mathit{\text{WSS}}}.$

The RMSE between B1 and FS was 8.6^° and 11.3^° (19.5% and 25.7% of the maximum value) for the mitral and aortic valve leaflet angles respectively, and 2.9 Pa and 4.3 Pa (13.4% and 16.4% of the maximum value) for the mitral and aortic valve leaflet $\overline{\mathit{\text{WSS}}}$ respectively. Between B2 and FS, the RMSE was 3.6^° and 3.9^° (8.1% and 8.8% of the maximum value) for the leaflet angles, and 4.7 Pa and 5.3 Pa (11.2% and 8.9% of the maximum value) for the leaflet $\overline{\mathit{\text{WSS}}}$. The total number of CPU hours for one cycle was 96, 144 and 192 for B1, B2 and FS respectively.

4. Discussion

In this study, a CFD model that used a custom FSI coupling algorithm to simulate the fluid-driven motion of the BMHVs in the Realheart TAH has been developed and tested across a range of relevant operating conditions, highlighting the ability to extract applicable and important device results, and showing the robustness of the model to a change in these conditions.

The model used an overset modelling approach to both facilitate the valve motion and maintain a consistent and refined mesh on the surface of the leaflets and constant mesh quality throughout the flow domain. This yielded spatially well resolved wall shear stresses on the surface of the leaflets, highlighted in Figure 12.

The approach used to blend strong and weak coupling was a novel application of the two existing fluid flow and structural coupling algorithms, resulting in a new simulation method for artificial heart valves. Most artificial heart valve simulations have used strongly coupled approaches, as methods using solely weakly coupled schemes have had mixed success regarding solution stability (Abbas et al. 2022). For this study, the blended weak-strong coupling and variable time stepping schemes were tuned to maximise computational efficiency whilst providing the minimum level of stability. Solution stability was established as solving until completion, rather than the solution diverging. The initial model used only weak coupling, and divergences in the solution arose at instantaneous valve opening - a problem also encountered by Borazjani et al. (2008). When strong coupling was introduced for small valve leaflet angles (<7^°), stable solutions could be obtained. The threshold angle of 7^° was the smallest threshold by which stability was maintained and maximised the computationally efficient portions of the cycle that were weakly coupled. Accuracy of the solution could be further improved by increasing the thresholds by which strong coupling was active, as was shown when comparing to a fully strongly coupled scheme and a second blended scheme with an additional threshold. B2 was included to capture leaflet rebound and fluttering, although this was not observed during the operating condition tested. B2 did however improve the accuracy of the solution whilst at a computational expense. Operating conditions, particularly those found in a full TAH, may induce leaflet fluttering and rebound due to higher fluid velocities and subsequent fluid forces. Therefore, the second threshold may be required in these cases.

Using a variable time stepping scheme is common when simulating MHV motion when using the ALE method, as it limits mesh motion, and subsequent remeshing between time steps (Bang et al. 2006; Choi and Kim 2009; Annerel et al. 2012). For the OD method used in this study, limiting mesh motion leads to improved interpolation of data between the component meshes of the valve leaflets and the background mesh of the fluid cylinder – improving stability of the solution (Ansys 2021). In addition, by decreasing the time step size at valve opening and closing, the temporal resolution of the leaflet wall shear stress improved, which varied greatest during this time.

The accuracy of this specific solution could not be validated, as there was no physical representation of the simplified geometry to test against. However, the implementation of the method into the full Realheart TAH model, and experimental validation of that model, is the subject of future work. Historically, CFD studies on BMHVs have investigated valves that simply open and close, where the housing is fixed in position, which is the same as the function of the aortic valve in this study. The range of values of ${\omega }_{\mathit{\text{opening}}}$ and ${\omega }_{\mathit{\text{closing}}}$ for the aortic valve in this model (15 to 35 rad s^-1) compares well against literature values of 11 to 50 rad s^-1 for ${\omega }_{\mathit{opening}}$, and 30 to 150 rad s^-1 for ${\omega }_{\mathit{\text{closing}}}$ (Choi and Kim 2009; Borazjani et al. 2010; Su et al. 2015; Mirkhani et al. 2016), noting that these studies use different operating conditions to those used in this investigation. There are no studies that consider a vertical translation of a BMHV in addition to its opening and closing, and so the values for $\omega$ cannot be compared directly to any of these studies, but it is noted that the values lie within the range of the literature for a stationary BMHV (10 to 40 rad s^-1 in this case).

There was a small rebound for both valves upon reaching fully open, a feature also observed by Borazjani et al. (2010). This feature was present in this model at 50 mm stroke length, but did not occur at the 25 mm cases, perhaps owing to the lower velocities and forces experienced at 25 mm compared to 50 mm. An increase in stroke rate and stroke length increased the translational velocity of the mitral valve, increasing fluid forces on the mitral valve leaflet, returning an increase in ${\omega }_{\mathit{\text{closing}}}$ and ${\omega }_{\mathit{opening}}$. For the aortic valve, only an increase in stroke length brought about an increase in ${\omega }_{\mathit{closing}}$ and ${\omega }_{\mathit{opening}}$, with the stroke rate changing very little. This may be attributed to the relatively small change in bulk fluid velocity between the 80 bpm and 100 bpm cases, as well as small variations in the cycle times in which the aortic valve opened and closed in each case.

The flow characteristics around the aortic valve were in line with what is normally observed, with the one central and two peripheral jet structures that correspond to the orifices when the valve is fully opened (Dasi et al. 2009), shown in Figure 7. The peripheral jets developed first and persisted for longer due to the fluid flowing around the outside of the opening valve leaflets. The central jet developed later, as upstream bulk fluid is pushed through the open valve by the closed mitral valve. When the mitral valve is fully open and moving upwards, the leaflets pierce through the stationary fluid in the valve’s vicinity and create a trailing flow structure in its wake. At 25 mm stroke length, the wake is short and surrounds the surfaces of the valve leaflet. At 50 mm however, the valve moves much further upwards causing the wake to detach.

These flow features translate to wall shear on the surface of the valve leaflet. In Figure 12(a)–(d), the leading edge of the aortic valve experienced elevated shear stresses, as the central jet moved through the open valve. The level of shear stress increased with an increase in stroke length and stroke rate, as these parameters equated to an increase in peak fluid velocities past the valve. For the mitral valve in Figure 12(e)–(h), at 25 mm stroke length, the trailing edge experienced the elevated shear stresses due to the wake surrounding this area, whereas at 50 mm, this occurred on the leading edge as the valve moved upwards through stationary fluid. It was observed that an increase in stroke length and rate returned an increase in $\overline{\mathit{\text{TAWSS}}}$ for both valves, attributed to the higher velocities caused by an increase in either of the two parameters.

The values of $Q_{\mathit{\text{peak}}}$ increased for an increase in stroke rate and stroke length, as the value of $Q$ at all times correlated directly with translational velocity of the mitral valve when it was closed. $V_{\mathit{\text{backward}}}$ remained relatively consistent between cases, indicating that it is the volume of fluid that moves back past the valve that determines its closure.

Pressures within the aorta varied as per the two-element Windkessel model that was applied at the outlet. The pressure variation followed the variation in $Q_{\mathit{outlet}}$, as the value of $R$ in the Windkessel model was the driving force behind the change in pressure. The effect of the value of $C$ can be seen at the ‘rebound’ pressure phase, where the aortic valve closes, and the pressure returns to the ambient pressure assigned to the Windkessel model – in this case 90 mmHg. The time taken for this rebound to return to ambient pressure was a function of the product of $R$ and $C$. The inclusion of the two-element Windkessel model highlights the possibility for further lumped parameter or 1D representation of downstream, or indeed upstream, vasculature to obtain a more realistic physiological pressure response. For full device simulation, these parameters can easily be altered to represent real patient data, and the response of the device to that specific patient could be analysed.

4.1. Limitations

In previous uses of overset meshing for evaluating artificial valve motion, authors have either chosen to include a gap model (Al-Azawy et al. 2016), where the cells in the gap between valve and valve housing are numerically ‘ignored’, or to explicitly model the gaps with a high-density region of cells (Nowak et al. 2020). In this study, there was no peripheral gap defined between the valve leaflet and valve housing, as a gap model was not available within the software, and a defined gap approach would have required approximately five times the number of elements, making model development challenging and time consuming. Peripheral gaps around the outside of mechanical heart valves do create high shear leakage flow, giving rise to the possibility of blood damage (Min Yun et al. 2014). Future studies where blood damage is the focus can include explicit gap modelling to account for this region.

The exclusion of a peripheral gap meant the fluid domain was divided into three separate closed sections at the point of either valve closing – the atrial, ventricle, and aortic regions. This is not a recommended approach, as the incompressibility of the fluid resulted in large pressure spikes during the period of instantaneous valve opening and closure, as seen in Figure 10 (Ansys 2021). These pressure transients did not affect the motion of the valves however, as the timeframe over which they occurred was very small (<1ms), and the subsequent transient pressure became ‘smooth’ after the initial phase of valve opening. Future studies will seek to overcome these unrealistic pressure transients, either through the implementation of a gap model or a small amount of fluid compressibility.

The choice of stroke length and stroke rates considered in this study were the same as would be found within the Realheart TAH. However due to the smaller geometry of the fluid tube in this study, the flow rates that were produced were much smaller than that found in a clinical setting and produced flow profiles that were in the laminar range. Adding a turbulence model is straight forward however and could be included in future studies that consider flow regimes that are turbulent.

5. Conclusion

This study presented a new methodology for the motion of bileaflet heart valves, which we have demonstrated on a simple positive-displacement pump using such valves. The method used overset meshing to both enable the fluid-driven motion of the valve leaflets, as well as maintain a consistent and refined mesh around the leaflets. This resulted in efficient capture of the valve leaflet kinematics across all ranges of motion, as well as high-quality wall shear stress capture on the surface of the leaflets. The proposed variable time-stepping scheme ensured computational efficiency, whilst improving the stability of the solution through the minimisation of mesh motion. The blended weak-strong coupling algorithm between the structural and fluid equations provided additional stability benefits, as it ensured only the most relevant range of motion required maximum computational stability, allowing for a more efficient computational approach for all other ranges of motion. This study demonstrated that the new methodology can be applied to the full range of physiologically relevant operating conditions, through the modification to the stroke length and stroke rate of the positive displacement mechanics. Velocity contours showed classical three jet flow features associated with bileaflet valves, angular velocities of the valves during opening and closing aligned with those found in the literature, and the variation in volume flow rate was typical of a pulsatile flow system. We have successfully identified a modelling approach that will allow for the efficient and accurate characterisation of full positive-displacement TAHs from both a haemodynamic and hemocompatibility perspective.

Disclosure statement

AN and ILP are employees of, or consultants to, and/or shareholders of, Scandinavian Real Heart AB.

Additional information

Funding

JB PhD funded 50/50 from EPSRC and Scandinavian Real Heart AB (Reference: 2426107)