Computational Fluid–Structure Interaction Study of a New Wave Membrane Blood Pump

Wave membrane blood pumps (WMBP) are novel pump designs in which blood is propelled by means of wave propagation by an undulating membrane. In this paper, we computationally studied the performance of a new WMBP design (J-shaped) for different working conditions, in view of potential applications in human patients. Fluid–structure interaction (FSI) simulations were conducted in 3D pump geometries and numerically discretized by means of the extended finite element method (XFEM). A contact model was introduced to capture membrane-wall collisions in the pump head. Mean flow rate and membrane envelope were determined to evaluate hydraulic performance. A preliminary hemocompatibility analysis was performed via calculation of fluid shear stress. Numerical results, validated against in vitro experimental data, showed that the hydraulic output increases when either the frequency or the amplitude of membrane oscillations were higher, with limited increase in the fluid stresses, suggesting good hemocompatibility properties. Also, we showed better performance in terms of hydraulic power with respect to a previous design of the pump. We finally studied an operating point which achieves physiologic flow rate target at diastolic head pressure of 80 mmHg. A new design of WMBP was computationally studied. The proposed FSI model with contact was employed to predict the new pump hydraulic performance and it could help to properly select an operating point for the upcoming first-in-human trials.


INTRODUCTION
As the prevalence of heart failure continues to rise over time with aging of the population, 51 left ventricular assist devices (LVADs) are life-sustaining therapeutic options that offer mechanical circulatory support in end-stage patients. 2 In particular, LVADs could act as a temporary bridge to heart transplantation 5 or as a destination therapy in patients that are not eligible for transplantation. 3,49 LVADs have to meet many requirements in terms of device implantability, 11,15,32 durability, 17,53 hydraulic power 29,36 and hemocompatibility. 6,39 All currently available LVADs are rotary blood pumps. They are small and reliable devices equipped with an internal impeller, which propels blood by rotating at elevated velocities, exerting high stresses on blood cells. 24,48 A new frontier in the LVAD applications is represented by the wave membrane blood pumps (WMBPs), developed at CorWave SA. 12 In this innovative pump model, blood propulsion is induced by the undulations of an immersed elastic membrane. 42 This wave propagation technology has the potential to pump blood at low stress conditions and consequently with a reduced risk of hemocompatibility-related adverse events. In addition, the pump hydraulic output of WMBPs can be regulated by manipulating the membrane wave frequency and amplitude, potentially adapting to patient specific demands, possibly mimicking the pulsatility of the native heart. 12 In Reference 34, a computational work aimed at studying WMBP performance was performed, based on the numerical solution of the fluid-structure interaction (FSI) problem arising between blood and wave membrane. The developed numerical model was em-ployed to perform three-dimensional simulations of the pump system at different pressure conditions and the results were validated against in vitro measurements.
In this study we aimed at taking further steps in the computational analysis of WMBPs, in the direction of the clinical application of the device. First, we examined a revised geometry design of the WMBP (J-shaped), that differs from that studied in Reference 34 primarily in the shape of the membrane, demonstrating improved performance in experimental testing, both in terms of hydraulic power and hemocompatibility. On the modeling side, a contact model 14 was added to study the membrane deformation at high oscillations, which may potentially lead to contact of the membrane with the pump housing. Finally, we provided a preliminary investigation of the potential hematologic impact of WMBPs, studying the distribution of the blood shear stress in the device, 4,41 and identifying non-physiological flow patterns, such as recirculation regions. 10,37

Wave Membrane Blood Pumps
The aim of LVADs is to reduce the work load of the left ventricle of failed hearts by taking over the blood pumping function. Thus, the wave membrane technology in WMBPs should overcome the pressure difference between the left ventricle (inlet) and the aorta (outlet) to pump the oxygenated blood into the circulatory system.
The structure of WMBPs comprises the following components (see Fig. 1): -an inflow cannula, implanted at the apex of the left ventricle; -an outflow cannula, placed at the opposite side of the device and in fluid communication with the ascending aorta via an outflow graft; -an actuator system, located at the center of the pump chamber, composed of a stator, equipped with two levels of electromagnetic coils, and of a mobile magnet ring; -a membrane assembly, that is suspended within the flow channel, displaced concentrically along the pump longitudinal axis; it includes a rigid membrane holder, which is directly connected with the magnet ring, and an elastomeric wave membrane, suspended in the so-called pump head, i.e., the region delimited by the pump head flanges.
The membrane design has been revised to improve pump efficiency and reduce the potential for blood trauma. In this paper we specifically considered two different pump geometries, showed in Fig. 2.
In the first pump design (A), extensively studied from a computational point of view in Reference 34 the membrane assembly shows a planar discoidal geometry with an annular membrane holder running along the external edge of the membrane disc (Fig. 2,  left). Hence, we refer to it as flat membrane pump design (or, concisely, flat design).
The new pump geometry (B) shows a longer membrane holder that bends towards the inlet up to the magnetic ring (see Fig. 2, right). This design is called Jshape membrane pump design (or J-shape design), due to the shape of the cross-section of the membrane assembly. The additional vertical portion of the membrane holder should work as a flow separator, distributing blood flow across the lower and upper sides of the membrane, and increasing the resistance to retrograde flow to improve hydrodynamic output. Thus, the J-shape design should reduce flow recirculation around the edges of the membrane holder and the exposure time to high shear rates. Other differences with respect to design (A) involve the wider clearance gap between the magnet and the actuator, to reduce local wall shear stress, and the smoother geometry of the housing and of the ring.
WMBPs are activated by an alternating current applied to the electromagnetic coils of the actuator, resulting in an oscillatory motion of the magnet ring and in a deformation wave propagating along the elastic membrane, which is damped by the effect of the surrounding viscous fluid. This results in an energy transfer from the membrane to the blood that ends up in pressure buildup in the flow channel, allowing blood FIGURE 1. Cross-sectional view of the main components of the implantable wave membrane blood pump, including inlet and outlet channel, the actuator assembly (stator, electromagnetic coils and magnet ring) and the membrane assembly (wave membrane and membrane holder).
to be pushed towards the outlet channel. Notice that the membrane approaches very closely to the pump flanges (see Fig. 1), isolating fluid portions in the flow channel, termed fluid pockets, and avoiding potential backflows.

Experimental Setup
Experimental in vitro data of the hydraulic performance of WMBP were collected in a re-circulating fluid loop, where the response of the system to different hydraulic conditions is measured for any operating point of the membrane. The pump characterization bench consists of a system of reservoir and tubing in which the hydraulic resistance can be modified by centrifugal pumps placed in a series circuit with the WMBP. A glycerin-water solution at 39% concentration in volume (37 C) was used to mimic blood. HQ curves (H head pressure, Q flow rate) combine the measurements of the head pressure and the corresponding flow at the outlet, measured with a pressure sensor (PendoTECH, France) and an ultrasonic flowmeter (Sonotec, Germany), respectively. The membrane was equipped with a sensorless position control system that makes the oscillation amplitude closer to a sine function via a feedback loop based on the actuator current. 47

Mathematical Formulation
The dynamics in WMBP were described by means of an FSI problem with contact, solved in the domain represented in Fig. 3.
For the fluid problem, we assumed the blood to be an incompressible, Newtonian and viscous fluid 34,44 with density q f and viscosity l f , and the incompressible Navier-Stokes equations were considered in the fluid domain X f . We remark that, since the Reynolds number in the pump varies between 600 and 1800 in space and time, there might be effects of transition to turbulence in the fluid dynamics that are not considered in this model, but that will be included in future works. The fluid domain changes in time due to the movement of the structures formed by two distinct bodies: the membrane assembly X m and the magnet ring X r . Following Reference 34, we assumed the membrane made of a homogenous and isotropic material with density q m , whose deformation is described by a linear infinitesimal elasticity model, with Lame´parameters k m and G m . On the contrary, the magnetic ring X r was modeled as a rigid structure, whose vertical movement has an impact on determining the fluid domain. The mechanical system connecting the magnet ring to the membrane, which consists of few supporting pins and screws, was omitted from the computational domain, as done in Reference 34. However, its effect on the flow dynamics will be object of future studies.
The fluid and structure problems were coupled at the interface R to guarantee geometric consistency and physical continuity of velocities and tractions. Accordingly, the FSI problem reads as follows: for each time t 2 ð0; T, find the fluid velocity u, the fluid pressure p, and the structure displacement b d, such that: T f ðu; pÞ n À T s ðdÞ n ¼ vn w on RðdÞ; ð1eÞ Dð b dÞ is first Piola-Kirchhoff structure tensor, coinciding with the Cauchy one due to the small displacement regime, and T f ðu; pÞ ¼ ÀpI þ 2l f DðuÞ, with DðuÞ ¼ 1 2 ðru þ ru T Þ, is the fluid Cauchy stress tensor. Quantities withbrefer to a Lagrangian framework, the other ones to an Eulerian framework. We have also highlighted the dependence of the fluid domain X f on the structure displacement d (geometric coupling) and indicated by n the external fluid normal.
Notice that in the continuity of traction (1e) we added a term v which modeled a repulsive contact force, acting in the direction of the wall normal n w , to take into consideration possible contacts between the wave membrane and the pump flanges. Specifically, we considered the relaxed contact model from Reference 22, that is a surrogate penalization method based on the assumption of an infinitesimal layer of fluid that always separates the fluid-structure interface from the wall boundary, as originally proposed in Reference 14. Therefore, we consider a contact layer of width e c , which extends inwards from the wall boundary see Fig. 4, left). In this framework, when the membrane penetrates into the contact layer, an artificial repulsive force is applied in the opposite direction to prevent collision and it is assumed to be proportional to the degree of penetration of the structure in the contact layer (see Fig. 4, right). In particular, given the distance gðxÞ between a membrane point x and the pump wall, we have   where c c is a positive constant and e c is a suitable a function of the local fluid mesh size h. Notice that the contact parameters e c and c c have not a physical meaning and they need to be calibrated on the specific application to correctly reproduce the physics of the contact dynamics, as done in other works. 22,25 In conclusion, we need to prescribe proper initial and boundary conditions. For the latter, we have: -the head pressure DP between the pump outlet C out and inlet C in , indicating the hydraulic resistance over the pump; -a sinusoidal displacement with frequency f and amplitude U=2 (with U being the oscillation stroke) on C m and C r , corresponding to the boundaries of the rigid complex formed by the membrane holder and the magnetic ring; -no-slip wall condition at the boundary C w , indicating the internal pump housing and the stator boundary.

Hemocompatibility
Hemodynamic quantities were derived from the FSI results to give a preliminary insight on the hemocompatibility of WMBPs. Specifically, we computed a scalar index of the blood shear stress, r scalar , based on the shear stress tensor r ¼ l f DðuÞ, used to provide indications about the risk of blood damage in the pump, as done in References 6, 33. In this work, it is defined as with k ¼ 1=3, in analogy with the Von Mises criterion for blood cells, as proposed in Reference 9 and used in many numerical works. 6,16,20,33 Alternative choices for factor k have been recently studied in Reference 21. In particular, the volumetric distribution of r scalar was investigated to identify potential critical regions in the pump domain. Moreover, the blood shear rate field were calculated to compare the risk of shear-induced blood trauma in the J-shape pump design with respect to the results in the flat design presented in Reference 34. We specifically computed the maximum shear rate in space and time value c max and its average in the pump head c avg PH . Finally, the wall shear stress (WSS) was considered to evaluate the risk of thrombosis at the pump contacting surfaces, as done in References 28, 45. In particular, we computed the percentage of the wall surfaces that presents values of WSS smaller than reference thresholds found in literature, that are rep-resentative for potential thrombus deposition. In particular, we used 0.1 and 0.3 Pa as thresholds, that are reference values for potential thrombus deposition taken from References 45 and 28, respectively.

Numerical Methods
For the numerical solution of the FSI problem we employed the extended finite element method (XFEM), 30 successfully used for several FSI problems with immersed structures. 1,26,35 In particular, we specifically referred to the numerical approach proposed in Reference 54 and applied to WMPB in Reference 34, in which a discontinuous Galerkin (DG) mortaring was employed to weakly couple the fluid and the solid subproblems at the interface. The XFEM-DG technique is an unstructured mesh method using a fixed background mesh with an imposed moving immersed structure mesh which cuts the fluid mesh generating general polygonal elements. The numerical solution on such cut elements is built by employing the standard basis functions defined on the uncut elements and doubling the degrees of freedom. 1,30 For this reason, it is particularly suitable for the problem at hand preventing large distortion of the fluid elements in the pump head. For details on the complete discretized problem, we refer the reader to Reference 34.
The novelty with respect to Reference 34 is the integration of the contact model, see (1e)- (2). XFEM framework allowed to easily integrate such a model in the corresponding weak formulation, see Reference 14.
The problem was solved in LIFEV, 8 a C++ finite element library that has been employed in a wide range of applications, particularly in hemodynamics, and validated against experimental data 19,40,43 and analytical solutions. 38 Specifically, the implementation of the XFEM-DG technique in LIFEV was validated and compared with standard fitted mesh methods in References 34, 54.
The linear system has been solved using a generalized minimal residual method (GMRES), with a block Gauss-Seidel preconditioner. In particular, we employed 7 cores (Intel Xeon E5-4610 v2) with 2.3 GHz frequency and a RAM of 252 GB.

Geometric Aspects
Exploiting the symmetry of the pump, for the numerical experiments we considered a reduced geometry of the pump system, restricted to a 120-degree section (see Fig. 5, left), in order to reduce the computational costs.
The construction of the reduced geometry introduces the artificial cut surfaces C f cut , C m cut and C r cut on the fluid, membrane, and ring domains, Here, we prescribed typical symmetry boundary conditions on the fluid and membrane unknowns, i.e., while on the rigid magnet ring we imposed the same sinusoidal displacement as on C r .

Space and Time Discretization
Starting from CAD geometry files, unfitted meshes were generated using GMSH, 27 see Fig. 5: a fluid mesh (520k elements), a membrane mesh (400k elements) and a magnet ring mesh (140k elements). In all cases, meshes were locally refined in regions of interest for the pump dynamics, such as the pump head, or with high curvature. Such meshes showed positive convergence results with respect to flow and pressure solutions computed on further refined meshes. Specifically, mesh convergence was studied across four different grids with decreasing mesh size h. The refinement ratio r i ¼ h iÀ1 =h i (h i being a representative fluid mesh discretization parameter for the ith mesh) was computed for each grid i. The indices used for the convergence analysis were the time-averaged flow rate Q at the outlet, the mean pressure px and the mean scalar stress rx at pointx 2 X f . Specifically, we took x ¼ ½0:0; 0:0; 1:4, which is located along the pump axis at the center of the pump head (see Fig. 5, left). Table 1 reports the values of such quantities and the corresponding approximated relative differences.
The time discretization is based on a uniform timestep Dt, that is chosen depending on the frequency of oscillation f. We studied the effect of the timestep on the solution for f ¼ 60 Hz by halfing it from 0.0004 s (42 time points per cycle) to 0.0002 s (84 time points per cycle). Analogously to what done for the grid sensitivity, Table 2 reports the convergence indices and the corresponding relative errors for flow rate, pressure and scalar stress. The results, featuring a maximum relative difference of 6:6%, indicate that it is sufficient to select the timestep so that 42 time points are taken for each period of oscillation.

RESULTS
In WMBPs the operating conditions of the device are defined by the following parameters: (a) the head pressure DP between the outlet and the inlet ports; (b) the frequency f of membrane oscillation; (c) the oscillation stroke U, corresponding to twice the amplitude of the membrane vibration.
In particular, we want to: (1) compare the J-shape membrane pump design (Design B) performance with that of the planar membrane (Design A), see In all the simulations, the fluid parameters were q f ¼ 1050 kg m 3 and l f ¼ 0:0035 Pa s, whereas, unless differently specified, the contact parameters were e c ¼ 0:1 mm and c c ¼ 0:04.

Design Comparison
We considered the same operational conditions used in Reference 34 for the planar membrane (Design A), detailed in Table 3. In this case, the oscillation parameters were fixed and we just varied the head pressure. We refer to this set of operating points as high-frequency (HF) points, where f ¼ 120 Hz. The pump system was simulated for T ¼ 0:025 s, corresponding to three periods of oscillation, using a fixed timestep Dt ¼ 0:0002 s. The elasticity parameters were q m ¼ 1125 kg m 3 , k m ¼ 2772 MPa and G m ¼ 56:6 MPa.
In Fig. 6 we report, at the instant of maximum membrane oscillation, the results for the operating point HF1 (DP ¼ 50 mmHg), in terms of velocity, pressure, and shear rate fields. From these results, we observe physiological values for all the quantities, in particular a favorable pressure gradient generated by the membrane undulations between the pump flanges and the outlet, that overcomes the adverse pressure gradient between inlet and outlet. In addition, the magnitude of the blood velocity in the J-shape design exceeds the values that were found in the flat design in Reference 34 (see Fig. 6d). Moreover, we notice that in the new pump design there are not the recirculation The values of convergence indices w ¼ Q i ; p x i ; rx i È É are presented for each mesh i, together with the corresponding approximated relative difference e w i ¼ jw i À w iÀ1 j=jw i j. Level 3 was selected for the rest of the paper since the errors are less than 2% on hydraulic quantities (e Q 4 , e p 4 ) and 5% for hemocompatibility index (e r 4 ). areas nearby the edges of the membrane holder and the magnet ring, observed in the old design in Reference 34. Indeed, the flow coming from the inner vein of the membrane (region A) is drawn into the outer vein (region B) following the direction of the pressure gradient. The vortex located in the outlet channel (region C) is continuously disturbed at high velocities, thus prohibiting the local formation of thrombus.
Also, we observe that the maximum shear rate c max is equal to 17,156.2 s À1 , in proximity to the membrane tip, while the volumetric mean value c avg PH at the same time instant is 369.5 s À1 . Both these values are significantly lower than proposed critical thresholds (e.g., 42; 000 s À1 for hemolysis 4,6 ) confirming what was found for Design A (c max ¼ 5640 s À1 , c avg PH ¼ 2209 s À1 ) in Reference 34. We remark that the higher value of c max in Design B is justified by the higher velocity  magnitude; nonetheless, the mean value of shear rate is lower than in Design A. Furthermore, the WSS was computed in the proximity to the magnet ring and shown in Fig. 7 for both the pump designs. The extension of the clearance gap between the magnet and the wall in the new pump design allows to reduce the average of the local WSS from 3.32 Pa in Design A, to 2.67 Pa in Design B.
In Table 4 we compared the hydraulic performance of the flat design (Design A), taken from Reference 34, and of the J-shape design (Design B), considered in this work. In particular, we computed the outlet flow rate Q and the hydraulic power W ¼ QDP. From these results we can observe the improved hydraulic outputs featured by Design B, in terms of both outflow rate and generated power.
We noticed that the contact model was never activated for the HF operating conditions considered in this section, because the membrane oscillations were not wide enough to enter inside the contact layer.

Analysis of the Contact Model
The contact model described in ''Numerical Methods'' section was introduced to reproduce possible contact dynamics between the wave membrane and the pump head flanges at high oscillations. Indeed, if the membrane deformation is particularly high (for instance at operating points with high stroke), the wave membrane may exit from the pump domain if no contact model is considered (see Fig. 8a). This results in a non-physical configuration that causes a sudden drop in the flow rate at the outlet. Instead, in presence of the contact model, the repulsive contact force v, defined in (2), is activated, prohibiting the membrane from exiting the fluid domain and penetrating into the wall (see Fig. 8b).
The effectiveness of the contact model depends on two parameters: the thickness of the contact layer e c and the penalty constant c c , see (2). If e c is too narrow or if c c is too low, the contact force enabled in the contact layer cannot push the membrane away, obtaining a result very similar to that reported in Fig. 8a. At the same time, too high values of e c and c c should be avoided as well, because the contact model may interfere with the natural wave propagation of the membrane deformation. In our simulations, we used e c 2 0:05; 0:2 ½ mm and c c 2 0:02; 0:2 ½ , depending on the magnitude of the oscillation parameters. In particular, for each set of oscillating parameters, the contact parameters were selected following a conservative approach, taking the smallest ones in these ranges that guaranteed non-penetration of the moving structure. Parameter values in the cited ranges allowed proper representation of the contact dynamics, with a variability in the mean flow at the pump outlet below 5%. Notice that values of contact parameters outside the given ranges may still be admissible, but have not been investigated in this work. Nevertheless, different choices of contact parameters are not expected to significantly affect the flow results, as long as non-physical phenomenon, such as mesh penetration or abrupt deviations from regular wave propagation are avoided.

Parametric Analysis
In this section, we report the results obtained for different operating points (OPs) of WMBPs, set by changing the head pressure (P-analysis), the frequency (F-analysis), and the stroke (S-analysis), see Table 5. Notice that the oscillation parameters shifted to lower frequencies and higher strokes than what was observed in ''Design Comparison'' section (see Table 3), because an increased hydraulic performance was found in this range for the J-shape design. We fixed the timestep Dt to 0.0004 s, so that we obtain 40 to 60 time points per cycle of oscillation for each tested frequency. The values of the membrane parameters q m , k m and G m are based on the material properties of the membrane. The specific values can be obtained from the membrane material suppliers, but can vary from lot to lot. * In Fig. 9 we reported the time evolution of the outlet flow rate Q in the P-analysis (left). As expected, Q decreases when DP increases, in accordance with the pump functioning. Specifically, the time-averaged outputs amounts to 3.45 L/min at DP ¼ 50 mmHg (P1) and to 2.66 L/min at DP ¼ 60 mmHg (P2). In the same figure, we also reported a comparison between the computed mean flow value and corresponding experimental in vitro data (right), obtained as described in ''Experimental Setup'' section, providing a preliminary result about the validation of our model.
Analogously, in Fig. 10, we reported Q for F-and S-analyses. These results show that the output increases when either one of the two parameters increases. Specifically, we have that the averaged-intime flow rate increases from 2.66 L/min for f ¼ 44 Hz (F1), to 3.09 L/min for f ¼ 50 Hz (F2) and further to 4.03 L/min for f ¼ 60 Hz (F3). Similarly, the averagedin-time flow rate amounts to 3.09 L/min for U ¼ 1:5 mm (S1), to 3.76 L/min for U ¼ 1:6 mm (S2) and to 4.39 L/min for U ¼ 1:7 mm (S3).
The membrane deformation was studied for different oscillation parameters by looking at the membrane envelopes, obtained by reducing the membrane crosssection to its centerline (see Fig. 3, right), and plotting its displacement during one cycle of oscillation. Membrane envelopes can provide important insight in terms of wave propagation and formation of fluid  We grouped them (possibly with overlaps) in three distinct sets: P where DP is changed; F where f is changed; S where U is changed.

*
The reader can contact the authors to ask for access to protected information on membrane dimensions and properties. pockets, by analyzing its distance from the pump flange. In Fig. 11, we reported the membrane envelopes for operating points F1, F3, S1, S3 (see Table 5). The black lines indicate the relative position of the pump head superior and inferior flanges, with respect to the superior and inferior edges of the membrane, respectively.
Notice that the distance of the membrane envelopes from the pump flanges decreases when f and U are higher (F3 and S3), especially with respect to the lower flange, indicating that the fluid pockets are better isolated during their propagation towards the outlet, thus limiting backflows. In addition, in all cases the amplitude of the membrane undulations increases while moving towards the center of the pump (decreasing n values), probably as a result of the decreasing thickness of the membrane in the radial direction. In particular, the maximum displacement is observed at the mem-brane tip, where contact occurs with the pump flange, without exiting from the pump domain (as discussed in ''Analysis of the Contact Model'' section).
Finally, we computed the solid Von Mises stress in the wave membrane to study the material resistance at different frequencies of oscillation. In particular, in Fig. 12 we showed this quantity for operating point F3 at time t ¼ 0:052 s and its maximum and volumetric mean values for each tested frequency. We can observe that the solid stress increased with the frequency of oscillation and that the maximum was achieved at the junction with the membrane holder (gray region) during the ascending phase of the membrane holder.

Hemocompatibility Analysis
Hemocompatibility-related adverse events in WMBPs, such as hemolysis or thrombosis, are closely  related to hemodynamic-generated stress exerted on blood cells. 23 In Fig. 13, left, we reported for F3 the fluid stress r scalar defined in Eq. (3) in the pump head, where the highest values of stresses are observed, especially on the pump head flanges, nearby the tip of the membrane, where the structure velocity is maxi-mum and where the membrane approaches closely to the wall.
In Fig. 13, right, we showed the time evolution of the volumetric mean of r scalar in the pump head during the last period of oscillation, for all operating points. Notice that the fluid stresses increase for larger values  of U and f, whereas, they decrease for increasing DP. The trend of r scalar is similar for all operating points, showing a peak at p, when membrane holder velocity is maximum.
In order to assess the potential of inducing insidepump hemolysis and thrombosis, we reported in Table 6 the statistics of r scalar and WSS for all operating points. In particular: the 'Max' column corresponds to the highest value of r scalar during the last period of oscillation; the 'Mean' column refers to the maximumin-time of the volumetric mean of r scalar in the pump head region; the area percentages represent the amount of wall surfaces exposed to WSS lower than 0.1 and 0.3 Pa; the volume percentages indicate the portions of pump volume with r scalar lowering 1 Pa, representative of low stress regions, or exceeding typical thresholds for hemocompatibility, i.e., 9 Pa for the degradation of the Von Willebrand factor, 24 and 50 Pa for the activation of platelets. 31 We remark that the area and volume percentages are computed at two different instants of the oscillation cycle: the WSS statistics refer to the time point t min that minimizes the volumetric mean of r scalar , when the oscillating structures have almost null velocity; while the volume percentages are calculated at the instant t max of maximum stress conditions, i.e., at maximum oscillation velocity.
Notice that all operating points show low stress conditions, with most of the pump volume (65-75%) with fluid stress smaller than 1 Pa and no occurrence of values exceeding 50 Pa. As a consequence, in our simulations WMBPs showed very low potential for shear-induced platelet activation or hemolysis, which are hypothesized to occur at shear stress conditions over 50 and 150 Pa, respectively. 24,33 Concerning the risk of thrombus deposition, the portions of the pump walls with low WSS (1.5-2.5% of the volume for WSS <0:1 Pa) are restricted to small regions below and above the magnet ring. Nonetheless, such regions are largely affected by the oscillation dynamics of the structures; indeed, in the same points, the WSS reaches physiologic values of 1-1.5 Pa 45 at time t max . Thus, being in absence of flow stagnation,  the potential for thrombogenicity is also expected to be low in WMBPs for the presented operating conditions.

Nominal Operating Conditions
The primary ambition of the WMBP system is to provide a physiological pulsatile flow to patients suffering of advanced heart failure. Therefore, the pump hydraulic and hemocompatibility performance should be studied over the entire heart cycle, with the head pressure DP varying in time. Here, for computational reasons, we limit the analysis to a fixed operating condition to initiate the process of prediction of pump performance using numerical models. In particular, we chose to study a typical operating condition (here referred to nominal operating point, NOP) generating physiologic mean flow at head pressure DP ¼ 80 mmHg, with oscillation parameters f ¼ 60 Hz and U ¼ 1:8 mm. This value of head pressure corresponds to the standard pressure difference between systemic and left ventricular pressure at diastole in a failed heart with LVAD support. 36 The pump system was simulated for T ¼ 0:05 s, corresponding to three periods of oscillation, using a fixed timestep Dt ¼ 0:0002 s, with contact parameters e c ¼ 0:2 mm and c c ¼ 0:2. Notice that the variation of the timestep and of the contact parameters with respect to default values in ''Results'' section was motivated by the larger membrane undulations observed during the simulations, probably due to the higher stroke parameter. The physical parameters are set as in ''Parametric Analysis'' section. Figure 14 shows the blood velocity and membrane displacement fields at t ¼ 0:0416 s (a), the trend in time of Q (b) and the membrane envelope (c) at NOP. The higher oscillation parameters of NOP resulted in a wider and more uniform membrane wave propagation, that generates a flow field with high blood velocities and no recirculation areas or stagnation points. The time-averaged LVAD output amounts to 6.25 L/min. In addition, maximum and volumetric mean of r scalar are plotted in time during the last period of oscillation (Fig. 15a). The maximum value in time and space is 56.4 Pa ,achieved at t 1 ¼ 0:035 s, while the volumetric mean in the pump head reaches a maximum of 5.3 Pa at time t 2 ¼ 0:0432 s, i.e., during the descending phase of the membrane holder. The histograms in Fig. 15b represent the volumetric distribution of the elemental scalar stress, i.e., the average value of r scalar in each element, for both time points. In particular, the highest elemental scalar stress (42.8 Pa at t ¼ t 1 ) is restricted to a volume of 0.003 mm 3 , close to the membrane tip. Since the mean local velocity in the element amounts to 114 cm/s, the residence time in the element is expected to be very low ($ 0:1 ms). Furthermore, for t ¼ t 2 , the stress values above 20 Pa are overall limited to less than 1 mm 3 in the pump volume.
In conclusion, we analyzed the relative contribution of normal and shear components present in the definition of stress r scalar , see Eq. (3). Figure 16 shows the ratio between extensional and shear stress components in the flow field at times t 1 and t 2 . The results highlight that the extensional stresses dominate in 45-55% of the pump volume, especially close to contact regions and in curved regions of the flow channel; on the other hand, shear is higher on most part of the pump walls and at the sides of the membrane holder and the magnet. When we restrict the analysis to regions with r scalar >9 Pa, shear stress was found to be predominant in 99:2% of the corresponding volume at t 1 and 71:4% at t 2 . Indeed, in the computation of r scalar the contribution of the normal stresses is penalized by a factor k ¼ 1=3, as used in most numerical studies. 6,16,20,33 However, other studies 18,21 proposed that extensional stresses actually deform red blood cells more than shear stresses and, hence, other values k>1 in Eq.
(3) could be tested in the future.

DISCUSSION
WMBPs, developed at CorWave SA, are a new family of cardiac assist devices based on an innovative technology where the propulsion of blood is obtained by the undulations of an immersed elastic membrane. This was introduced mainly to try to overcome the limitations of standard rotor-based LVADs, in particular to provide a pulsatile flow rate to the body and to reduce the risk of hemocompatibility-related adverse events.
We believe that the research on this new pump and the study of the influence of its design and operating points on its performance could be of interest for the bioengineering and medical community.
For such reasons, in this work we have presented a computational study of a new design of WMBP (Jshape membrane pump) and compared its performance with that of a previously studied pump (flat membrane pump). 34 .
Our numerical results, in good agreement with the experimental findings (see Fig. 9, right), showed a drastic increase in the hydraulic performance of the Jshape design with respect to the flat one. Such improvement could be ascribed mainly to the curved and elongated membrane geometry that allowed to separate the flow across the two sides of the membrane holder, reducing flow recirculations around the magnet, which causes hydraulic energy dissipation, differently from what was found for the flat design pump (see Fig. 6 in Reference 34). At identical operating conditions, the increase in the flow conditions led to higher shear rate conditions, which were however still significantly lower than the critical threshold of 42,000 s À1 : 4 The J-shape design produced the same hydraulic output of the flat one at lower oscillation frequencies (44-60 vs. 120 Hz) with a beneficial effect on the pump hemocompatibility (lower fluid shear stress) and an increase fatigue life of the membrane (lower solid Von Mises stress).
In our numerical results, we observed that the pump output increased of 0.35-0.45 L/min, for every increment of 5 Hz of oscillation frequency, while it increased by 0.6-0.7 L/min, when the stroke is increased by 0.1 mm. These data can then be used to further refine and optimize the membrane design and performance, in view of future in-human trials. In particular, we tested the nominal operating point, with higher membrane oscillation parameters (f ¼ 60 Hz, U ¼ 1:8 mm), that allowed to recover physiological output (' 5=6 L/min) at DP ¼ 80 mmHg, thanks to a better isolation of the fluid pockets during their transport in the pump head. Nonetheless, the operating point should be tested also at different head pres-sures to fully characterize the HQ trend and capture the global behavior of the pulsatile pump. With the aim of reducing the computational time and explore several scenarios, some studies considered lumped parameter models to rapidly provide qualitative information of the pump performance over the heart cycle for variable head pressure. 52 Another new contribution of the present work was the introduction of a contact model to handle at the numerical level the recurring impingement between the membrane and the pump flanges at high oscillations. Indeed, we considered here more realistic operating conditions of the membrane than what seen in Reference 34, with an increase in stroke of 40-70%, that allowed for a more effective transport of the fluid pockets towards the outlet without backflow, but lead to contact between membrane and flanges. We also identified suitable ranges for the contact parameters e c and c c , that allow to obtain a contact force that prohibits non-physical membrane penetration out from the pump boundary (see Fig. 8), without affecting significantly wave propagation and pump outflow.
Finally, we want to discuss the hemocompatibility of the new pump. In LVADs blood shear stresses are considered to be good indicators of the risk of hemocompatibility-related adverse events.
In particular, large fluid shear stress peak values are frequently markers for hemolysis potential, especially in case of long exposure time. In literature, the reference thresholds for hemolysis risk typically vary from 100 to 650 Pa, depending on the exposure time. 24,33,46 In our simulations, we observed reduced levels of blood shear stress, with most of the pump volume (65-75%) presenting values below 1 Pa, while never exceeding 60 Pa. Specifically, at nominal operating conditions, the maximum value (in time and space) of scalar stress is 56.4 Pa, significantly lower than the critical thresholds, and the pump volume exposed to stress higher than 40 Pa amounts to 0.0063 mm 3 . Moreover, the peak stress value is achieved on the pump head flange close to the membrane tip, in a region where the flow is very disturbed and the residence time is expected to be low.
Such results are particularly promising when compared with the hemocompatibility studies on modern LVADs working in similar flow-pressure conditions. For instance, in HeartMate II, a clinically-approved axial blood pump, 11 mm 3 of pump volume was exposed to stress levels higher than 150 Pa at 5.4 L/min, 13 but showed acceptable levels of plasma-free hemoglobin in experimental studies at similar operating conditions. 7 In the same experimental work, 7 the centrifugal blood pump HVAD showed lower hemolysis indices than HeartMate II, even though in Refer-ence 16 0:7% of its volume presented stresses values higher than 100 Pa at 5 L/min. Finally, the preclinical tests for HeartMate III, 13 one of the most recently developed LVADs, showed very good performance in both in vitro hemolysis testings and in in vivo bovine trials, even though in the numerical stimulations the scalar stress was found to be higher than 50 Pa in 52.7 mm 3 and than 150 Pa in 3.3 mm 3 (up to 350-400 Pa). Nevertheless, in future works, an hemolysis model will be developed in order to properly assess the risk of hemolysis and estimate the mean residence time in WMBP.
Similarly to hemolysis, there exist stress thresholds that are thought to be significant from the point of view of thrombogenesis: for instance, the unfolding of the Von Willebrand factor protein, that mediates the platelet adhesion, is triggered at fluid stress levels higher than 9 Pa, 24 while platelets are activated at shear stress higher than 50 Pa. 31 Both these thresholds are very rarely achieved in our pump, and only locally in space and time, see Table 6. Low fluid stress conditions should also be avoided , especially at the blood contacting surfaces, because they increase the likelihood of platelets deposition and aggregation, that can lead to thrombosis. 50 In our numerical results, the WSS statistics indicate that only 1.5-2.5% of the pump walls present very low WSS values (smaller than 0.1 Pa), especially when the frequency of oscillation is low. Nevertheless, the WSS in the same regions was found to be one order of magnitude higher at other time points, indicating that the low stress conditions are local in time and not accompanied by flow stagnation during the oscillation cycle.
Some limitations characterized this study. First, we omitted from the computational domain a set of thin suspension springs that are employed to maintain the magnet ring centered around the actuator. Their inclusion may deform the flow patterns near the magnet ring. Second, the membrane holder, that is actually a separate component embedded in the elastic membrane, is here modeled as part of the membrane assembly, as done in Reference 34. In principle, this simplification may affect the membrane displacement close to the junction of the membrane holder. Further, although we observed strains smaller than 10%, the assumption of linear elasticity for the membrane could be improved in future works by considering a finite elasticity model to better describe the membrane dynamics and its interaction with blood. Finally, we observe that we have used fluid-dynamics quantities (such as shear stresses) to obtain preliminary results about hemocompatibility. Although shear stresses could provide some indications in such a direction as highlighted by some works, 24,24,31,33,50 these are surrogate quantities that in general may or may not be definitively predictive of hemolysis generation in the device. More direct analyses provided by the introduction of an hemolysis model could be developed for future works to properly assess the risk of hemolysis and estimate the mean residence time in WMBP.

AUTHOR CONTRIBUTIONS
MM: conceiving of the presented idea, design of the computational framework, implementation of the code, numerical simulations, first draft and writing of the manuscript, discussion of the results. FC: writing of the manuscript, discussion of the results, experimental planning. CV: conceiving of the presented idea, design of the computational framework, writing of the manuscript, discussion of the results, project supervision.

DATA AVAILABILITY
The experimental data used for the validation of the numerical results in Fig. 9 are publically available on Zenodo (https://doi.org/10.5281/zenodo.4964021).

CODE AVAILABILITY
Not available.

CONSENT FOR PUBLICATION
Not applicable.