1 Introduction

Congenital heart defects are the most prevalent type of birth defect, occurring in approximately 9 out of every 1,000 live births (1). Among these, 9%–12% are classified as single ventricle physiologies (2). Patients with single-ventricle physiologies often have atrial and ventricular septal defects, leading to the mixing of arterial and venous blood within the heart chambers. A portion of the venous blood enters the systemic circulation without being oxygenated, causing organ hypoxia. Typically, patients undergo three-stage palliative surgeries: the Norwood, Glenn, and Fontan procedures (3). After the Fontan procedure, patient's superior and inferior vena cava (SVC and IVC, respectively) are connected to the pulmonary artery (PA), forming a total cavopulmonary connection (TCPC) structure. This setup forces the single ventricle to pump blood through both systemic and pulmonary circulations in the absence of a functional right ventricle. Over time, this can result in high central venous pressure and low pulmonary artery pressure (known as Fontan failure) (4). Approximately 40% of patients eventually suffer from single ventricle pumping failure (5).

Assistive devices play a critical role in augmenting circulatory power for these patients. In 2010, Throckmorton et al. (6) introduced an axial blood pump into the IVC to aid blood flow to the lungs. However, as the outlet of the blood pump is oriented towards the SVC, the blood flow pressurized by the pump collides with the blood flow within the SVC. This collision leads to an increase in blood pressure inside the SVC and hinders blood return. Then in 2013, Throckmorton et al. (7) implanted axial blood pumps in both SVC and IVC to provide dual-lumen assistance. In 2014, Wang et al. (8) used a paired umbrella double-lumen cannula (Maquet Cardiovascular, Wayne, NJ) and extracorporeal blood pump CentriMag (Thoratec, Pleasanton, CA) for total cavopulmonary assist. Animal studies revealed that this setup effectively improved hemodynamics in failing Fontan sheep. Subsequently, in 2019, Wang et al. (9) combined a double-lumen cannula with a two-valve extracardiac conduit and CentriMag to achieve effective support. In 2016, Gandolfo et al. (10) designed a T-shaped mechanically assisted TCPC model using the Jarvik Child 2,000 axial pump (Jarvik Heart, New York, NY, USA), which showed promising results in vivo. Although these techniques can provide total cavopulmonary assist, they have drawbacks such as the need for multiple devices, limitations on patient activity due to extracorporeal blood pumps, and the necessity for vascular reconfiguration.

To address the issues mentioned above, this study proposed a novel right heart assist device, as shown in Figure 1. The main components of this device are made of flexible materials. During minimally invasive surgery, the device is in a contracted state with a smaller diameter. Once implanted in the human body, the device expands and begins to function. It is anticipated that the deformation of flexible components during operation will impact the device's performance. Consequently, ensuring rigidity in its expanded state while maintaining flexibility during folding poses the principal challenge in developing this device. The material selection for each component is crucial. This research focuses on a key component—the blades, investigating their optimal elastic modulus. A one-way FSI simulation was employed to study blades’ deformation during rotation, and their impact on blood pressure increase, hemolysis risk, and thrombosis risk. The optimal elastic modulus is determined through a comprehensive analysis of these indicators.

Figure 1. Contraction state and post-implantation state of the flexible right heart assist device. (1. bearing seat; 2. shaft; 3. protective stent; 4. impeller; 5. constrictor ring; 6. motor; 7. cable; 8. inferior vena cava; 9. superior vena cava; 10. pulmonary artery).

2 Materials and methods 2.1 Flexible right heart assist device and computational models

Figure 1 shows the proposed right heart assist device, consisting of a protective stent, a bearing seat, a bearing, a shaft, an impeller, a constrictor ring, a motor, and cables. The protective stent is made of Nitinol. One end of the stent is fixed to the bearing seat, while the other end connects to the constrictor ring and can move axially. The impeller is fabricated from flexible materials such as polyurethane through an integrated molding process. It consists of two parts: the hub and the blades, with the base of the blades affixed to the surface of the hub. The hub's upper end is mounted on the shaft, and its lower end can move axially. The blades can adopt two configurations: curled and unfolded. When the assist device contracts, the lower end of the hub moves downward, and the blades curl. As the device expands, both the hub and the blades revert to their initial shape. The motor is powered by cables, and drives the rotation of the shaft and impeller. During surgery, the entire device is compressed within a delivery catheter and implanted into the TCPC via the femoral vein. Once released, the protective stent expands to tightly fit against the inner walls of the vena cava, stabilizing the device. The rotating impeller aids in directing blood within the SVC and IVC into the PA. Additionally, the hub is uniquely designed in a spindle shape to prevent blood from colliding. The design goal of this assist device is to provide a pressure rise of 2–25 mmHg for the total systemic blood flow (3l /min 11) in Fontan patients when operating at 4,000–8,000 rpm (12).

This paper focuses on flexible blades, therefore adopts simplified computational models. The TCPC model is an ideal circular tube model derived from patient's data (13). According to the principle of equal area, the diameters of the SVC and IVC are set at 13 mm, while those of the LPA and RPA are set at 9 mm. In CFD simulation, only the impeller is retained. The base of the blades is fixed to the hub's surface to ensure the model's consistency with reality. The minimum and maximum diameters of the impeller are 0.4 mm and 10.9 mm respectively, the number of blades is 12, and the thickness of the blades is 0.2 mm. In FEA simulation, only the blades are retained. The deformation of the hub is ignored. The impeller and blade models presented in this paper are the optimized versions after structural optimization. Figure 2 illustrates the final computational models used in this study.

Figure 2. CFD and FEA simulation model. (A) Ideal TCPC structure; (B) blades of right heart assist device.

2.2 One-way FSI methodology, boundary conditions and solving settings

In this study, one-way FSI simulation was used to analyze the flexible blades’ deformation and its effect on the device performance, as shown in Figure 3. First, ANSYS Fluent 2020R2 (Canonsburg, USA) was used to model the fluid dynamics in TCPC and obtain the blade surface pressure. Subsequently, these pressure data were utilized in ANSYS Structural 2020R2 (Canonsburg, USA) to compute blades’ deformation. The deformation of the blades was quantitatively described by the displacement of nodes on the blades. After the deformed shape of the blades was exported, the fluid domain was remeshed. Subsequently, the CFD simulation was performed again. The cycle was ended when the difference in blood pressure rise between two consecutive CFD simulations was less than 5%. Finally, the last CFD result was used to calculate the hemolysis risk and thrombus risk.

Figure 3. One-way FSI simulation process.

In CFD simulation, flow boundary conditions were defined at the SVC and IVC inlets, with flow rates of 1.2l /min and 1.8l /min respectively (11). The LPA and RPA outlets were set as pressure boundary conditions, with a pressure value of 1,750 Pa, corresponding to the average pulmonary artery pressure during a patient's cardiac cycle (13). Blood was modeled as a Newtonian fluid with a density of 1,060 kg/m3 and a dynamic viscosity of 3.5 cP. This article analyzes the operation of the assist device at its lowest and highest rotational speeds. The flow domain is divided into rotating and stationary regions. The rotational speed of the rotating region is set at 4,000 rpm and 8,000 rpm. Furthermore, rotational boundary conditions are applied to the surfaces of the hub and blades, with the rotational speed settings matching those of the rotating region. The remaining surfaces are treated as stationary boundaries. Additionally, the SST k-ω model was used to simulate turbulent conditions, as this model can accurately resolve the flow fields both at the wall and in the bulk flow regime (12, 14, 15). The fluid domain was meshed as a polyhedral grid using Fluent Meshing 2020R2 (Canonsburg, USA). After grid independence study, the fluid mesh consisting of approximately 1.35 million elements was used for CFD simulation. Inflation layers were also used in the fluid meshes to achieve the SST k-ω turbulence model requirement (y + less than 2). A steady-state solver was employed, with the convergence criteria was set to a maximum residual of below the 1 × 10−4 threshold.

In FEA simulation, the blade material was set based on the polyurethane material properties. There is a wide variety of linear elastic polyurethanes, with elastic moduli ranging from 6.2 MPa to 27,600 MPa (16). A pre-simulation (a single round of one-way FSI) was conducted to narrow down the range of elastic moduli for the study, as shown in Figure 4. At 4,000 rpm and 8,000 rpm, the maximum deformation of the blades showed a consistent pattern. When the elastic modulus is within 50 MPa, the maximum deformation decreases rapidly as the elastic modulus increases. Between 50 MPa and 200 MPa, the decrease in maximum deformation is more gradual. Above 200 MPa, both the absolute value and the variation in maximum deformation are extremely small. Based on these results, the elastic moduli selected for the study were 10 MPa, 20 MPa, 30 MPa, 40 MPa, 50 MPa, 100 MPa, and 200 MPa. In the actual configuration, the base of the blades is fixed to the hub's surface, leading to the application of a fixed constraint at the base of the blades. The structural mesh of the blade, generated using HyperMesh (Altair, USA), is a hexahedral mesh and consists of two layers along the thickness direction. After the grid independence study, the structural mesh consisting of 36,096 elements was used for FEA simulation. A transient solver was employed, with a total computational duration of 1.1 s and a time step of 0.001 s. The convergence criteria was set as changes of the maximum deformation, the minimum deformation, and the average deformation are all less than 1%.

Figure 4. Maximum deformation of blades with different elastic moduli at rotational speeds of 4,000 rpm and 8,000 rpm. (These results were obtained from one round of one-way FSI simulation).

2.3 Hemolysis risk assessment model

During the high-speed rotation of the impeller, red blood cells are destroyed due to shear forces, resulting in hemolysis. The hemolysis risk can be quantitatively evaluated using the Normalized index of hemolysis (NIH). Numerical simulations and hemolysis experiments employ different calculation methods to derive this index. In numerical simulation, NIH is commonly derived through a power law model, as detailed in Equations 1–3 (17–19).

σ=(16∑(σii−σjj)2+∑σij2)12(1)

D=∑inletoutletC⋅σα⋅ΔTβ(2)

NIH=0.00015∗D(3)

where σ is the scalar stress, ΔT is the exposure time, and C,α,β are constants obtained from hemolysis experiments. Mitamura et al. (20) derived a set of parameters based on the experimental data from axial blood pumps: C=3.62×10−5,α=2.416,β=0.785. Heuser et al. (21) used experimental data from a Couette viscometer to fit another set of parameters: C=1.8×10−6,α=1.991,β=0.765. The working principle of the impeller in this paper differs considerably from that of the axial blood pump and Couette viscometer, and is closer to that of the centrifugal blood pump. Therefore, this paper uses experimental data of centrifugal blood pumps, provided by the U.S. Food and Drug Administration (FDA), to determine the parameters of the hemolysis risk assessment model (22).

In addition, Blackshear et al. (23) conducted a hemolysis experiment and proposed the relationship between shear stress experienced by red blood cells and their exposure time:

(τrbc2)(t)=C1(4)

In the above equation C1 is a constant of proportionality. The results of Blackshear et al. indicate that the hemolysis risk is proportional to the exposure time and the square of the shear stress. Combined with Equations 1–4, the final hemolysis risk assessment model used in this paper is as follows:

σ=(16∑(σii−σjj)2+∑σij2)12(5)

NIH=∑inletoutletA⋅σ2B⋅ΔTB(6)

In Equation 6, A represents the overall scaling factor, set at 2.7 × 10−10. This value was derived by multiplying the scaling factor C, provided by Heuser (21), with the scaling factor of 0.00015 from Equation 3. The parameter B, set at 0.55, was determined from experimental data, as detailed in the Supplementary Document.

2.4 Thrombus risk assessment model

Thrombus formation is another common complication associated with rotary assist devices (24). It involves the generation, transport, and cascade reactions of numerous components. The adhesion and activation of platelets are critical steps in thrombus formation. Cell adhesion only occurs in areas where the wall shear strain rate (γ˙) is lower than 100 /s (25). Therefore, the shear strain rate can be used to identify areas prone to thrombus growth (26). Its calculation method is shown in Equation 7. An accelerated thrombus formation model is utilized to quantitatively assess the thrombus risk (27). This model includes three convection-diffusion equations (Equations 8–10), describing the generation and transport of non-activated platelets (NP), activated platelets (AP), and adenosine diphosphate (ADP). The concentration of AP is ultimately used to represent the thrombus risk.

γ˙=[2((∂u∂x)2+(∂v∂y)2+(∂w∂z)2)+(∂u∂y+∂v∂x)2+(∂u∂z+∂w∂x)2+(∂v∂z+∂w∂y)2]12(7)

∂(ϕa)∂t+(u⋅∇)ϕa=Da∇2ϕa+(8)

∂(ϕn)∂t+(u⋅∇)ϕn=Dn∇2ϕn−(9)

∂(ADP)∂t+(u⋅∇)ADP=DADP∇2ADP+RADP(10)

class="mb15">In the above equations, Dx is the respective diffusivity coefficient, AM and AC denotes mechanical activation and chemical activation, and RADP is the amount of ADP contained in a platelet. Platelet activation (change from NP to AP) occurs either with mechanical cues (Equation 11) or chemical cues (Equation 12).

AM(ϕf,τ)=(1−ϕf)C11β1β1ϕfβ1−1β1τα1β1(11)

AC(ADP)={ADPADPt⋅tADP,ADP≥ADPt0,ADP<ADPt(12)

where ϕf is the ratio of the number of activated platelets to total platelets count, α1, β1 and C1 are the power law model parameter, ADPt is the threshold for chemical platelet activation, tADP is the characteristic time of platelet activation, and τ is the shear stress calculated from Equations 13–15.

τ=13σxx2+σyy2+σzz2−σxxσyy−σyyσzz−σxxσzz+3(σxy2+σyz2+σxz2)(13)

σij=μ(∂ui∂xj+∂uj∂xi)−ρui′uj′¯(14)

−ρui′uj′¯=μt(∂ui∂xj+∂uj∂xi)−23ρkδij(15)

In the above equations μ is the dynamic viscosity, μt is the eddy viscosity, ρ is the fluid density, k is the turbulent kinetic energy, δij is the Kronecker delta, ui is the velocity component and ui is the spatial component.

3 Results 3.1 Variation of blades’ deformation and pressure rise during FSI simulation iterations

Figure 5; Table 1 illustrate the changes in blades’ deformation and blood pressure rise during FSI simulation iterations. Figures 5A,C,D sequentially show the blades’ deformation during three rounds of FSI simulation. It is evident that as the iterations progress, the solid body gradually aligns with the dashed frame. Additionally, the maximum deformation of the blades gradually decreases, indicating that the blade shape is stabilizing. According to Table 1, the pressure rises in both the SVC and IVC also tend to stabilize during the iteration process. The pressure rises of the last two iterations differing by less than 5%. These results validate that the one-way FSI simulation used in this paper can achieve a stable blade shape and the corresponding pressure rise.

Figure 5. Blades' deformation during the one-way FSI simulation iterations. (A) Blades' deformation obtained from the first round of FSI; (B) Blades' deformation under stress; (C) Blade deformation obtained from the second round of FSI; (D) Blade deformation obtained from the third round of FSI. (Dotted frame represents the blade shape before deformation, solid body represents the blade shape after deformation).

Table 1. Maximum blade deformation and blade pressure rise during three rounds of one-way FSI simulation.

Moreover, stress exists within the blade in the second and subsequent rounds of FSI simulation. Figures 5A,B were compared with each other to verify the accuracy of applied stress. Figure 5A shows the blade bending under surface pressure. Figure 5B depicts the blade returning to vertical shape after the removal of surface pressure, while stress is retained. The maximum deformations in Figures 5A,B are 1.55 mm and 1.505 mm, respectively. The difference is only 3%, indicating that the stress is sufficiently precise.

3.2 Blood pressure rise corresponding to blades with different elastic moduli at two rotational speeds

Figure 6 displays the blood pressure rise corresponding to flexible blades and rigid blades. DP-Total represents the total pressure rise, obtained by the weighted sum of DP-SVC and DP-IVC, based on the flow rates in SVC and IVC (DP-Total = DP-SVC × 0.4 + DP-IVC × 0.6). At a speed of 8,000 rpm, blades with an elastic modulus of 10 MPa deform excessively and lose their functionality entirely. Consequently, data related to this condition are not included in the analysis.

Figure 6. Blood pressure rise corresponding to blades with different elastic moduli at two rotational speeds. (A) At 4,000 rpm; (B) at 8,000 rpm.

From Figure 6A, it can be observed that at 4,000 rpm, the pressure rises corresponding to flexible blades show an increasing trend. As the elastic modulus increases, DP-SVC, DP-IVC and DP-Total all gradually increase. The pressure rises corresponding to flexible blades with elastic moduli of 50 MPa, 100 MPa, and 200 MPa are relatively close to those of rigid blades. Figure 6B reveals that at 8,000 rpm, the pressure rises corresponding to flexible blades also trends upward. However, when the elastic modulus of the blades increases from 40 MPa to 50 MPa, the pressure rise decreases instead. The pressure rises corresponding to flexible blades with elastic moduli of 40 MPa, 50 MPa, 100 MPa, and 200 MPa are relatively close to those of rigid blades. Notably, the pressure rise corresponding to the 40 MPa elastic modulus flexible blade is higher than that of the rigid blades. Comprehensive analysis indicates that in most cases, blade deformation leads to a decrease in functionality. However, at the highest rotational speed, there are instances where blade deformation results in an increase in functionality.

3.3 Hemolysis risk corresponding to blades with different elastic moduli at two rotational speeds

Figure 7 shows the hemolysis risk associated with flexible blades and rigid blades. After completing the streamline number independence test, 1,100 streamlines were released at both the SVC and IVC inlets. The NIH values for each streamline were calculated and average values were used to represent the hemolysis risk in the blood flow of the SVC and IVC (NIH-SVC and NIH-IVC). The overall hemolysis risk (NIH-Total) was derived by the weighted sum of NIH-SVC and NIH-IVC (NIH-Total = NIH-SVC × 0.4 + NIH-IVC × 0.6). The results indicate that at both rotational speeds, the hemolysis risks associated with flexible and rigid blades are relatively similar. At 4,000 rpm, the NIH-Total for flexible blades ranges between 7.26 × 10−4 and 7.67 × 10−4 g/100 L, with a maximum difference of 5.6%. The NIH-Total for rigid blades is 7.03 × 10−4 g/100 L. Therefore, the difference in NIH-Total between flexible and rigid blades ranges from 3.3% to 9.1%. At 8,000 rpm, the NIH-Total for flexible blades ranges between 15.01 × 10−4 and 16.65 × 10−4 g/100 L, with a maximum difference of 10.9%. The NIH-Total for rigid blades is 15.48 × 10−4 g/100 L.As a result, the difference in NIH-Total between flexible and rigid blades ranges from 0.2% to 7.6%.

Figure 7. Hemolysis risk corresponding to blades with different elastic moduli at two rotational speeds. (A) At 4,000 rpm; (B) at 8,000 rpm.

Under the same rotational speed, flexible blades with different elastic moduli have varying shapes, yet the corresponding hemolysis risks are similar. This suggests that the impact of rotational speed on hemolysis risk is greater than that of blade deformation. Moreover, literature data indicate that the NIH values for HeartMateII (Thoratec Inc., Pleasanto, California, USA) and CH-VAD (CH Biomedical, Inc., Suzhou, China) at medium speeds are 5.83 × 10−3 g/100 L and 1.35 × 10−3 g/100 L, respectively (28). The highest NIH-Total in this study is 1.665 × 10−3 g/100 L, falling between the NIH values of HeartMateII and CH-VAD. This result demonstrates that the assist device proposed in this paper is reliable in terms of hemolysis.

3.4 Thrombus risk corresponding to blades with different elastic moduli at two rotational speeds

Figure 8 presents the results of thrombus formation simulation, corresponding to a rotational speed of 4,000 rpm and a blade elastic modulus of 10 MPa. As mentioned earlier, the shear strain rate (SSR) is utilized to determine areas prone to thrombus growth. The concentration of AP is used for quantitatively assessing the thrombus risk. To enhance readability, the AP concentration is scaled and represented as the SAP index. This index has a minimum value of 0, corresponding to the initial AP concentration (25 × 1012 /m³). Each unit increase indicates an increase of one millionth of the initial AP concentration. Figure 8A shows the distribution of wall SSR. It is observed that regions with SSR less than 100 /s are only present on the walls of the SVC, LPA and RPA. There are no low SSR areas on the blades’ surfaces. Figure 8B, based on Figure 8A, displays grids where the SAP exceeds 23 in a dark red color. Detailed examination through close-up and cross-sectional views reveals that all four vessels exhibit wall grids with higher SAP. However, IVC lacks regions with low SSR, and the high SAP at the LPA and RPA outlets is due to the transport of AP. Only the wall of the SVC displays areas with both low SSR and high SAP, indicating that thrombus is more likely to form on the SVC wall. Figure 8C shows streamlines released from the IVC. Due to the larger flow volume in the IVC, a portion of the blood enters the SVC, leading to a collision within the SVC. This collision area has a low SSR, which facilitates the deposition of AP and thus the formation of thrombus. At both rotational speeds, the thrombus formation scenarios for different blades are similar to those depicted in Figure 8B. Thrombus are more likely to growth on the wall of SVC, hence separate illustrations are not provided.

Figure 8. Simulation results of thrombus formation. (Corresponding to rotational speed of 4,000 rpm, blade elastic modulus of 10 MPa. (A) Distribution of wall shear strain rate; (B) thrombus distribution; (C) IVC streamline diagram).

A statistical analysis of the mean SAP values in the SVC for different blades provides the following insights: At 4,000 rpm, the mean SAP value for flexible blades with various elastic moduli ranges between 2.60 and 2.68, with a maximum difference of 3%. The mean SAP value for rigid blades is 2.60. As a result, the maximum difference in the thrombus risk between flexible and rigid blades is 3%. At 8,000 rpm, the mean SAP value for flexible blades with different elastic moduli ranges from 14.92 to 15.88, with a maximum difference of 6%. The mean SAP value for rigid blades is 15.24. Therefore, the maximum difference in thrombus risk between flexible and rigid blades is 4.2%. These results indicate that, in terms of thrombus formation, the impact of rotational speed is also much greater than that of blade deformation.

4 Discussion

The Fontan surgery successfully diverts venous blood into the lungs by constructing the TCPC structure. The issue of mixing between arterial and venous blood has been resolved. However, long-term reliance on a single ventricle to power both systemic and pulmonary circulations can lead to various complications, such as elevated central venous pressure, increased pre-load on the single ventricle, and insufficient pulmonary perfusion. Patients are at high risk of developing single ventricle heart failure. Numerous scholars have proposed a variety of solutions to provide additional power for the Fontan circulation. Based on the location of the assistive device, these solutions can be categorized into extracorporeal and intracorporeal assistance. Trusty et al. (29) used pulsatile blood pumps, conduits, and baffles to explore efficient extracorporeal assistance methods in an in vitro loop. Their findings revealed that without partitions, blood pumped out of the pump's outlet tube could be reabsorbed by the inlet tube, causing recirculation. Effective power assistance was only possible by completely separating the inlet and outlet tubes. Wang et al. (8, 9) used dual-lumen catheters and extracorporeal blood pumps for Fontan support. The successful assistance achieved was also due to the use of valves or umbrella-like membranes that separated the inlet and outlet of the dual-lumen catheters. Therefore, providing extracorporeal assistance to patients requires the additio