Flow Characteristics and Coherent Structures in a Centrifugal Blood Pump

Blood clot formation can be initiated by local flow conditions where regions of high shear and long residence time regions, such as flow separation and stagnation, have been identified as risk factors. This study highlights coherent structures, some of which not yet considered in the literature that may contribute to blood clot formation in the ECMO (Extra Corporeal Membrane Oxygenator) circuit. The centrifugal ECMO pump investigated in this study is compact and delivers adequate volume of blood with relatively high pressure in order to compensate for the large pressure drop in the membrane oxygenator. These requirements lead to regions with high shear in several different parts of the pump. In the narrow gap between the pump house and the impeller body (the magnet) a Taylor-Couette-like flow is observed with azimuthally aligned wavy vortices, which are also pushed towards the bottom of the pump-house by the flow generated by the blades. At the bottom gap between the impeller house and the pump house one finds spiraling flow structures, due to the rotation of the former structure. Separation bubbles are found near the tongue of the pump and at the lee sides of the blades. Such vortical structures have in literature been identified as regions where platelets may be activated whereby clots may develop.


Introduction
In certain severe medical situations, the function of critical organs may temporarily be lost or insufficient to maintain life. When no other treatment option exists, a heart and lung machine can provide temporary support to allow for organ recovery. This mechanical system consists of artificial replacements that manage the blood pumping function of the heart and/or the oxygenating function of the lungs. ECMO (Extra Corporeal Membrane Oxygenation) therapy is such a life-saving treatment shown to be efficient in treatment of the patients suffering from acute respiratory failure with or without cardiac failure. Introduced Niclas Berg niber@mech.kth.se 1 Linné FLOW Center & BioMEx, KTH Mechanics, SE-100 44 Stockholm, Sweden in 1976 by Bartlett et al. [1], ECMO was initially intended for infants. However, the H1N1 influenza pandemic in 2009 showed that ECMO treatment is effective in adults as well, supported by the results presented in two controlled randomized studies reporting improved survival rates in adults [2] as well as newborns [3].
Due to the non-physiological blood flow conditions in these systems, the use of mechanical circulatory support, such as ECMO, is not free of complications. The complications are associated with blood damage including increased risk of hemolysis, activation of platelets, thrombosis and emboli. Thrombus have been observed in tubing connectors, cannulae, at the entrance of the oxygenator as well as in the blood pumps [4][5][6][7]. Regarding blood pumps, thrombus formation and hemolysis have been reported to be induced independent of the type of the pump (i.e axial, roller or centrifugal) [4,8,9]. In the ECMO setting, centrifugal pumps have been reported to be less hemolytic, but more thrombus prone as compared to roller pumps [4,7]. Both complications are highly linked to the local flow environment where regions of low shear rates (< 100 s −1 ) and pathologically high shear rates (> 10 000 s −1 ) are both problematic [7,10,11]. Considering blood clot formation, initiation of platelet activation has been reported to occur in regions of high shear and stagnation, but also in recirculation zones where the time during which the platelets are exposed to shear becomes an important factor [12]. Another flow characteristic identified as a potential platelet activator is vortical flow structures that allows blood clots to develop [13]. Risk of hemolysis on the other hand is increased during events when the pump is exposed to negative inlet pressures, elevating the risk of cavitation [14,15]. Due to the design of the pumps combined with them operating at high rotational speeds, the above mentioned flow characteristics will appear, especially regions of high shear.
Since the mid 1990s, CFD has been used to investigate the flow field in blood pumps [16][17][18]. Commonly, Reynolds Averaged Navier-Stokes (RANS) solvers combined with different turbulence models have been adopted, for which a detailed review is provided by Fraser et al. [19], stressing the importance of performing detailed flow field characterization using Large Eddy Simulations (LES) or Direct Numerical Simulations (DNS). Moreover, although blood is known to be non-Newtonian, most studies assume a Newtonian fluid behavior due to the large shear rates appearing in these pumps [5,[20][21][22]. Recent publications, focused on the flow features and its association with hemolysis and pump thrombosis, have been presented [5,21,22] where Wu et al. [21] combined numerical simulations with clinical observation for HeartMateII (axial blood pump). However, studies of the details of flow structures formed in the pump has yet not been reported although studies of centrifugal pumps used in other engineering applications emphasizing that LES is favorable in order to capture the complexity of the flow field found in these pumps [23][24][25].
In this study, a centrifugal blood pump is investigated, whose design is characterized by a magnetically levitated motor without bearings on which the impeller blades are mounted. In the lower part of the pump where the rotating cylindrical magnet is enclosed by the pump housing, two regions appear resembling two, for the fluid dynamic research community, classical configurations; one characterized by a rotating cylinder placed in a non-moving larger cylinder volume and the other a rotating circular plate placed at a small distance above a non-moving wall. These configurations ressemble the Taylor-Couette flow and the appearance of Ekman and Bödewadt layers, respectively. In 1890, Couette observed that for low inner cylinder rotation rates, an azimuthal flow will develop in the small gap between two concentric cylinders with either both cylinders rotating or only one of the cylinders rotating while keeping the other still. With increasing rotation rate, instabilities occur, leading to the formation of toroidal vortical structures superimposed on the azimuthal flow, described both experimentally and theoretically in the 1920s by Taylor [26,27]. This set the basis for a research area that thereafter has attained great attention for a wide range of parameters (Reynolds and Taylor numbers) [28][29][30][31][32][33][34][35][36].
Another intricate flow feature is expected to occur in the lowermost part of the pump where the rotating pump house drives the flow in the gap to the stationary pump enclosure. Simplified versions of such flows have been investigated thoroughly in the last century; first by von Karman [37] who calculated the flow for an infinite disc placed in a fluid at rest, and later extended by Batchelor [38] for the flow between two rotating discs. For finite configurations with the rotating plate/plates enclosed by a cylindrical wall, the basic flow [38] and its stability have been extensively studied [39][40][41][42][43][44]. Three limiting cases are commonly referred to as: the von Karman layer found in the flow of a rotating disc in a quiescent fluid, the Bödewadt [45] layer observed for a rotating fluid over a stationary plate and, finally, the Ekman layer [46] found when the fluid is rotating at an angular velocity similar to that of the rotating plate.
The flow in the centrifugal pump is far from ideal due to the interactions between the different regions of the pump. For example, the flow generated by the impeller blades introduces flow into the gap between the magnet and the housing. This flow has in turn an impact on the flow structures and thus also on the shear stresses appearing in the flow field. As mentioned above, the level of shear appearing in the flow affects platelet activation and red blood cell damage inducing thrombus formation and hemolysis. The work presented here is inspired by our previous work suggesting a non-invasive method for blood clot detection in a centrifugal ECMO pump from the acoustic spectra of the flow [47]. In that study, swirling flow motion at inlet and outlet was observed associated with a low frequency (2-10 Hz) that was amplified in the presence of a blood clot in the circuit. However, the study did not reveal the location of the clot nor its initiation. Thus, the aim with the current study is to gain insight in the details of the flow structures forming in the centrifugal blood pump, in particular focusing on the vortical structures that may act as platelet activators, increasing the time platelets are exposed to shear, normal as well as elevated levels.

Method
The flow in the centrifugal ECMO pump is assumed to be incompressible, turbulent and Newtonian. The Large Eddy Simulation equations are given by the following: where ρ = 998 kg/m 3 is the fluid density, ν = 0.89 · 10 −6 m 2 s −1 the kinematic viscosity (water), p and u are the space filtered (i.e. resolved) pressure and velocity, respectively, and τ SGS represents the subgrid stress. Although whole blood shows non-Newtonian behavior, the viscosity in this study is assumed to be constant. The motivation adopting this assumption lies in the dependency of viscosity on the local blood cell concentration and shear rate diminishing with increasing shear rate. At shear rates exceeding approximatelyγ = 100 s −1 , the fluid exhibits a Newtonian behavior [48]. The shear rates found in the pump clearly exceed this threshold. Thus, the choice of assuming that the fluid is Newtonian is well motivated.
In this study, the commercial software Star CCM+ was used. Equations 1 and 2 were discretized using blended central differences combining second order upwind with a central difference scheme for the convective terms. The PISO algorithm was used for pressure correction and the solution was evolved in time through an implicit second order time stepping scheme. 1440 time steps per impeller rotation period were used throughout the simulations, obtaining a Courant number less than 1 for all cells in the core region of the computational domain. We use this restriction not for stability reasons, but rather in order to maintain the accuracy of the small scales of the flow and the computation of the shear rates. The subgrid stress was modeled implicitly through numerical dissipation and no additional model is included (τ SGS = 0), since the spatial resolution was found to be fine enough to resolve about one order of magnitude of the inertial subrange.

Geometry and boundary conditions
The pump geometry is shown in Fig. 1a. The outer wall of the pump forms a closed surface, and the impeller blades are mounted on a levitating cylindrical magnet. The geometry was generated in a CAD software, and is based on the design of the Thoratec CentriMAG TM pump and will be referred to as the Baseline geometry. A second geometry was generated with the distance between the outer radius of the casing and the magnet housing increased from 0.5 mm (Baseline geometry) to 1 mm (Widened gap geometry), similar values as for the centrifugal pumps investigated by Fraser et al [5].
The inlet and outlet pipes were extended (see Fig. 1a) to allow the ingoing flow to develop and to reduce the influence of the outlet boundary condition. Flow rates and pressures obtained at inlet and outlet in experiments were applied as boundary conditions. Given the relatively high Reynolds number of the incoming flow (Re ∼ 10 4 ), a turbulent flow is expected. Hence, synthetic turbulence with 5% intensity and a length scale of 0.5 mm was added as described in [49]. A constant pressure condition was employed at the outlet. At the walls, no-slip condition was applied.
The geometry was discretized with an unstructured mesh of polyhedral cells. Local refinement was applied near the blades and all sharp edges of the geometry. In order to  mesh (b, c). Arrows indicate flow and rotation direction. The insets (i1) and (i2) show zoom-ins of the mesh for the widened and the baseline geometry, respectively resolve the boundary layers, 5 to 10 prism layers were added near the walls. Three different meshes were used, obtained by successive global refinement. The total number of cells were 3.8M (coarse), 7.7M (medium), and 15.8M (fine), respectively. Cutplanes of the medium mesh is shown in Fig. 1b and c. To account for the impeller rotation, a sliding mesh approach was applied. The moving reference frame (MRF) approach was also considered and tested. However, although more computationally efficient, the methodology was found to be unable to capture the rich dynamics in the impeller-tongue interaction.
For the sliding mesh approach, the computational domain is split into three parts: the inlet pipe, the volute along with the outlet pipe and the region surrounding the impeller. The two former were kept stationary throughout the simulation whereas the mesh surrounding the impeller was continuously rotated, corresponding to a rigid body rotation with a constant rotation rate. Based on our previous work [47], pump rotation rates of 1500, 2000, 2500, 3000 and 3500 rpm were considered. A summary of the considered cases is shown in Table 1.

Phase and time averaging
For all simulations, velocity and pressure data were sampled 240 times per revolution. Each case ran for 20 rotational cycles, requiring 3h on 128 processors per cycle, i.e. each simulation allocated a total of 7680 CPUh. To obtain statistics of the flow, phase averaging was applied. For a timeseries u(t 0 ), u(t 1 ), ..., u(t N ) the mean,ū i , and the rms, u rms i , at phase i (time = t i ) are respectively introduced as, and, where N p is the number of samples per period, and N s is the total number of samples for the considered phase. For the current simulations, the phase averaging was performed at the passing frequency of the four larger impeller blades, i.e. N p = 60. At least N s = 80 samples were used for each phase average.

Mesh sensitivity
The sensitivity of the results to the mesh resolution was investigated by considering the total pressure gain over the pump, as well as the velocity sampled over lines. The average increase in pressure was evaluated by computing the spatial average of pressure over planes at the inlet and outlet of the pump (see the inset in Fig. 2), and then taking the time average of the difference over 10 rotation cycles. The comparison of the results shown in Fig. 2 showed no significant difference between the medium and fine meshes. The velocity profiles were evaluated over a line in the magnet center (A-A') and in the lowermost gap (B-B'). These locations were chosen due to the coherent structures developing here, of importance to resolve in detail for the aim of this study. The velocity was linearly interpolated to 100 equidistant points, and then time averaged over 10 rotation cycles. The profiles for the case with a rotation rate of 3500 rpm are shown in Fig. 3. A maximal relative difference of around 4% is found for the line A-A' and 2% for line B-B' between the medium and fine mesh when normalized by the inlet velocity. A larger difference is found between the coarse and the fine mesh.
Given the difference in results between the fine and medium mesh, the results will henceforth be reported for the medium mesh.

General flow description
The phase averaged velocity for two different planes for a rotation rate of 2500 rpm is shown in Fig. 4. Considering the incoming flow, the majority of the flow is deflected radially At the outlet of the pump (location V), a clear separation region can be noted. This is more clearly shown in Fig. 5 where the root mean square of the velocity magnitude is shown for different rotation rates. Elevated fluctuation levels can be found in the impeller wake and downstream of the tongue region near the outlet. Furthermore, it should be noted that only a small part of the impeller induced radial flow leaves the pump through the outlet. Instead, the majority of the flow re-circulates back into the pump along its upper part. This description is similar for both geometries investigated here.
Due to the magnet-mounted impeller, two narrow regions are introduced (locations III and IV). The impeller passage forces the fluid down into region III, creating a downward directed flow. This flow enters IV where the rotation of the magnet induces a secondary flow. These two latter regions will be discussed in more detail below.

Coherent structures in the gap between the magnet and pump housing
The narrow gap between the outer radius of the magnet and the pump presents a geometry resembling the well-known Taylor-Couette flow. For Taylor-Couette flows the most  important (non-dimensional) parameter has been identified to be the Taylor number, T a = ωdR i /ν, where ω is the angular velocity, d is the gap of 0.5 mm, the inner cylinder radius R i = 16 mm, and ν is the kinematic viscosity. Depending on the rotation rate of the two cylinders as well as whether the cylinders are co-rotating or counter-rotating, different flow regimes are distinguished [28,31]. In this study, the values of the Taylor number are ranging from 1300 to 3000. Thus, considering the work by Andereck [28], the onset of a wavy vortex flow, modulated waves and turbulent Taylor vortices are expected to be observed in our flow setup.
At low Taylor numbers, viscous forces dominate over centrifugal forces. The flow is aligned with the azimuthal direction and the velocity is increasing linearly with radius from the outer to the inner cylinder. As the centrifugal forces increase, this base flow becomes unstable and counter-rotating flow cells appear. Figure 6 shows snapshots of the azimuthal vorticity along the centerplane for 1500, 2500 and 3500 rpm, where elongated vortex streaks are found in all cases. As shown, the number of vortices remains constant with increasing rotation rates. Furthermore, increasing the gap width from 0.5 mm (Baseline geometry) to 1 mm (Widened gap geometry), leads to an increase in vertical wavelength of the vortical structures, or a decrease in the number of observed azimuthal vortices for all pump rotation rates.
The stability of the vortical structures can be greatly influenced by superposing a coflow on the Taylor-Couette flow, as pointed out by Tsameret and Steinberg [36]. While the flow may become unstable, the instability (at moderate Taylor number) shows a convective behavior, i.e. the instability grows in the reference frame of the flow. Therefore, perturbations can be convected out of the system before having grown substantially. In this study, a co-flow is generated by the passage of the impeller blades forcing the fluid into the narrow gap. The time-averaged vertical velocity (z-direction) is shown in Fig. 7. On average, the flow displays a downward motion in most parts of the domain with a varying velocity. A minor upward flow can be observed around the area located below This flow behavior is enhanced with increasing rotation rates. Moreover, increasing the gap width, the region closest to the impeller passing is characterized by a stronger downward motion as compared to the Baseline geometry. However, further down (decreasing z), this downward motion decreases to become less strong as compared to the Baseline geometry. The residence time of platelets and red blood cells in the Taylor-Couette region will be determined mainly by the flow rate, which, as shown in Fig. 8a, increases with the rotation rate. Only a small difference can be observed between the two different gap widths. An estimate of the flow through time through the region can be obtained as t ∼ h/ū, where h = 11 mm is the height of the region, andū = Q/(π((R i + d) 2 − R 2 i )) is the mean vertical velocity. The obtained results are presented in Fig. 8b.
The timescale of pump rotation period, is between t rot = 0.04 s (1500 rpm) and 0.017 s (3500 rpm). Comparing this to the convective timescale shows that fluid particles will typically complete 3-4 revolutions in the Taylor-Couette region for the Baseline geometry, and 6-8 with the Widened gap.

Bödewadt layer instability
The vortical structures advected out of the Taylor-Couette region enter the lowermost region of the pump, a region comprised of an upper rotating plate and a lower stationary plate separated by a distance of 1 mm. This configuration resembles that associated with the development of both a Bödewadt layer [45] and an Ekman layer [46] as mentioned in the introduction of this article. The radial and azimuthal velocity distribution in the lowermost region is shown in Fig. 9 for the Baseline geometry. In the center of the gap, the centrifugal forces are balanced by the pressure gradient, yielding no net acceleration. Approaching the lower wall, the flow decelerates due to the no-slip condition at the wall, while the pressure gradient remains unchanged. The result is a net acceleration directed radially inwards. Along the upper wall, opposite to the above, the velocity gradually increases, providing an outward going flow. Both the upper and lower Ekman layers are separated. The lower layer constitutes a Bödewadt layer reported to become unstable at lower rotation rates as compared to the upper layer, the Ekman layer [43]. For all cases studied here, spiraling structures are found within this layer, shown by the iso-contours of the λ 2 criterion in Fig. 10. As indicated by Fig. 10, there was little difference between the two geometries investigated in terms of the spiraling motion developing in the lowermost region of the pump.

Conclusions
The numerical simulation of the flow in the centrifugal pump reveal interesting regions that could not be observed in steady state RANS computations or in experiments (due to the lack of optical access). Our initial URANS studies (not published) did not reveal any indication of the vortical structures in the gaps. Among the most important flow structures, undulating Taylor-Couette vortices was observed and captured with LES. Moreover, these vortices have a slow vertical flow component leading to long residence times. The flow in the core of the magnet house is intricate, characterized by unsteady fluid motion along the axis of the magnet with low speed at both directions. In the narrow gap between the rotating magnet house and the pump house itself, a complex set of vortices are formed. Regions of stagnation and high shear, as well as vortical structures have previously been identified as potential sites where blood clots can develop. In this study, high shear regions as well as elongated vortical structures have been identified. Such regions with prolonged residence time have been found in the narrow regions surrounding the levitating magnet. These structures show similarity with either the well-studied Taylor-Couette or the Ekman layer flows. The significance of the different structures for the blood pump is the impact on the thrombogenicity of the ECMO system. As small particles (such as platelets) tend to accumulate at the center of vortical structures, it is suspected that these structures may contribute to the blood clot formation that has been observed in such centrifugal pumps. This point is further investigated in Fuchs et al. [50]. The LES results are highly beneficial in order to identify the two main components contributing to hemolysis and platelet activation, namely the strength and location of the shear stress and the residence time of the fluid transporting platelets and red blood cells.