Inducer and diffuser blade geometries as well as rotor blade tip-to-shroud gap were assumed fixed. The inducer was designed (SolidWorks, Dassault Systems, Walton, MA) with four blades parallel to the flow axis (i.e., inducer blade exit angle,
= 0), in order to prevent prerotation at rotor inlet. The diffuser was designed with five blades; the metal angles at inlet
and exit (
) were set at 66° and 0°, respectively, the latter in order to minimize whirl at pump exit.
The limitations in rotor magnet volume and losses in magnetic flux, both arising from flow area constraints, have been circumvented by transporting the primary flow path to the center of the rotor magnet cylinder. Blades were fixed to the centrally axial shaft (hub) and to the inside surface of the magnetic cylinder (tip). Smaller height blades were installed on the outside surface of the magnet cylinder to encourage forward flow in the secondary flow path between the magnet cylinder and the shroud (Figure 1). Rotor length and rotor magnet length
were taken as 40mm and 10mm, respectively, and rotor gap clearance was set at 0.2mm for all blade designs.
Design Variables and Constraints
Independent design variables were rotor speed
and rotor blade geometry. The parameters relating to rotor geometry were tip (rt) and hub (rh) radii (which, together, defined Ac and mean blade radius, rm), and rotor blade metal angles at inlet (βb1) and exit (βb2). Dependent design variables were shaft torque (T), pump flow (Q) and head (ΔH), overall hydraulic efficiency (ηh), area-averaged shroud shear stresses (SSsh), and backflow (BF) and forward flow (FF) expressed as percent of total forward flow.
Nominal pump operation point was set at Q = 5L/min and ΔH= 100 mmHg. Minimum acceptable value for ηh was 20%. Maximum allowable limit of SSsh was 550 Pa to prevent hemolysis3; and those for BF were 15% and 35% at rotor inlet and outlet, respectively, to minimize turbulence.
Torque was obtained from fluid momentum change across rotor blades4,5
where Δvur, the change in fluid tangential velocity across rotor blades, was obtained from velocity triangles (Figure 2A)
Tangential velocity at the median plane was defined as u = ωrm; and fluid exit angle from rotor blades was assumed 20 < βf2 < 30. Axial fluid speed, vx, was determined from flow rate and area
. Torsional-to-hydraulic power transfer was expressed using Euler work equation
Combining Equations 1 and 2
From force balance (Figure 2B), pressure rise (ΔP) across the blade (rotor or diffuser) was expressed as
where ρ is fluid density,σ is rotor blade solidity (chord–pitch ratio), and CD is drag coefficient. Mean fluid angle βm was defined as
where the subscripts i and o refer to blade inlet and outlet, respectively. Kinetic energy rise
across the rotor was calculated assuming constant Ac
Blade pressure losses caused by drag forces were obtained from Bernoulli equation
Pump efficiency was formulated as
Efficiency was expressed as a function of Q, and maximized with respect to Q by setting
= 0 (Appendix)
is the real root of Equation 6 within range [3 < Q(L/min) < 10], which maximized ηh. Diffuser-to-rotor solidity ratio
was obtained from blade geometry, and the 1×9 arrays Pr and Pd are functions of βf2 (Appendix). The rotor-to-diffuser drag coefficient ratio (
) was calculated after computational fluid dynamics (CFD) results became available (Section 2.6). Equation 4 was written separately for rotor and diffuser and the resulting equations were combined with Equation 5. After rearranging, rotor drag coefficient was expressed as
and plotted as a function of Q.
Motor battery voltage and armature resistance, R, were accepted as 12 Volts and 0.011 Ohm, respectively. Armature current range was selected as 0.29 <ia (Amp) < 0.37. Neglecting mechanical losses (ηm = 1), torque constant Kb was calculated by substituting6
in Equation 1 and rearranging
The number of stator coil turns, N, was estimated from
Where Bg, magnetic field intensity at fluid gap, was assumed to be a function of magnet material (Nd-Fe-B, N48-40 grade) and pump body material (titanium).
Beginning with ηh = 25% at the nominal pump operation point (Q = 5L/min and ΔH = 100 mmHg) and varying βf2, rhub, and rtip within the specified flow range, corresponding values of ω were calculated from Equation 3 using Matlab (The Mathworks, Natick, MA). Knowing and Q, fluid attack angle at rotor inlet (βf1) was calculated from velocity triangles (Figure 2A). Initially, rotor leading edge metal angle (βb1) was equated to βf1, assuming zero rotor incidence (ir = βf1 − βb1). Finally, Kb and N were solved from Equations 7 and 8, respectively.
The virtual prototype (SolidWorks) was fed in a commercial Computational Fluid Dynamics software (CFD, ANSYS Inc., Canonsburg, PA), which computed the nominal values (at Q = 5L/min, designated with the subscript “nom” from here on) of BF in primary and secondary flow channels at rotor inlet and exit, SSsh,nom, Hnom, and Tnom. The value of ω was iterated in CFD until experimental and nominal values of ΔH converged. At this point,ηnom was calculated from Equation 2. This provided the first ω -Q- ΔH- ηhcombination of the actual characteristic performance curve of the pump. For the range of 0 < Q(L/min) < 10, the ηh-Q and ΔH-Q plots for constant ω were obtained in Matlab, at increments of 0.5L/min. The flow rate, which maximized ηh, was identified as Q@ηmax and used for solving Equation 6 (Section 2.4). Mean value of ηh for the operational range was calculated from
The effect of rotor trailing (βb2) and leading (βb1) edge metal angles was tested in Models 1A–1F and 2A–2D, respectively. The effect of Ac on flow field and performance was tested in Model 3A. The respective values of βb2, βb1, and Ac, for which ηmax was highest, were selected as the optimum for that design variable.
Computational Fluid Dynamic Analysis
Fluid domain was discretized into an 8-million-cell mesh to solve momentum and Navier–Stoke equations. Assuming turbulent characteristics at rotor inlet and across diffuser blades (Re > 10,000), k-ε turbulence model was adopted; and in order to locate the nodes as close to the viscous sublayer as possible, enhance wall treatment was used.
Left ventricular assist device flow domain was divided into three zones representing inducer, rotor, and diffuser. Frame rotation was applied to rotor zone, and rotor blade boundaries were defined as rotationally moving walls with zero relative velocity. All other surface boundaries were defined as nonslip stationary walls. Pump inlet and outlet cross-sections were defined as mass-flow inlet and pressure outlet, respectively; and boundary intersections at rotor inlet and outlet were defined as interface. Shroud shear stresses were averaged over areas where shear exceeded 400 Pa.
Following the methodology described in Design Iterations, blade geometry optimization was conducted at the nominal design point.
Initial Shaft Speed and Blade Height
To approximate the starting point for blade geometry optimization by CFD, Equation 3 was manipulated using Matlab. When Ac was augmented 30% by increasing blade tip height from 7 to 8mm, ω changed from 9,923 to 8,627rpm (−13%; Table 3).This, however, shifted Q@ηmax, from 9.73 to 14.00L/min (+44%), i.e., farther away from the nominal flow. The location of Q@ηmax shifted right by ~20% with each 10% increase in Ac. Increasing theoretical fluid exit angle (βf2) from 20° to 30° (+50%) brought Q@ηmax left by ~8%, while increasing ω by only 1–2% (Table 3). To begin CFD iterations, rt and rh were selected as 4 and 7mm, respectively, and ω was set at 10,122rpm.
At this speed, shaft torque was 3.3 < T(mN-m) < 4.2 at the design point. With Kb = 11.4 mN-m/Amp and N = 9 turns/cm stator length, corresponding pump power range was 3.49 < W (Watts) < 4.44.
Fluid attack angle at rotor inlet (βf1) was determined as 81° from velocity triangles; and rotor blade attack angle (βb1) was initially set at this value.
Optimization of blade geometry at rotor exit
Remaining at the Qnom, rotor blade trailing angle was increased stepwise, starting from βb2 = 0° (Table 4). Shaft torque and pump head were observed to peak in Model 1E with βb2 = 14° (Tnom = 6.828 mN-m, ΔHnom = 107.2 mmHg), but ηnom was highest (= 19%) and backflow at rotor outlet was lowest (BFro, nom = 32.69%) in Model 1A with βb2= 37°. On the other hand,ηh -Q curve peaked in Model 1B with βb2 = 32°, while maximum attainable efficiency ηmax = 22.81% occurred at Q@ηmax = 8.20L/min. In Models 1A and 1E, highest ηmax coincided with 7.60 and 9.20L/min, respectively (Figure 3).
Optimization of blade geometry at rotor inlet
Continuing with Model 1B at the nominal flow rate, CFD iterations of blade attack angle, βb1, were initiated. Backflow at rotor inlet was observed to minimize (BFri,nom = 6.88%) as βb1 was lowered to 74° (Model 2A); but ηnom (20.74%),ΔHnom (108.7 mmHg), and Tnom (5.837 mN-m) maximized when βb1 was further lowered to 67°, 64.5°, and 58°, respectively (Table 5). The peak of the ηh -Q curve was highest in Model 2B (Figure 4), with ηmax = 24.33% (6.7% improvement over Model 1B) observed at Q@ηmax = 7.60L/min (0.20L/min closer to Qnomversus Model 1B).
Velocity streamlines exhibited no prerotation at rotor inlet, no whirl at diffuser exit, and minimal disturbance at magnet inlet and exit (Figure 5). Flow was measured to split smoothly between central (through the magnetic cylinder) and peripheral (through the magnet-coil gap) channels of Model 2B (Table 6). Approximately 10% [479.8/(479.8 + 4331.8)] of the total flow passed through the peripheral channel, where total backflow accounted for 9.59% (47.9/499.5) of forward flow. Maximum area-averaged shroud stress was 514 Pa.
Optimization of rotor blade height
In a last iteration, rh of Model 2B was increased from 4mm to 4.5mm, and Model 3A was created. Reduced Ac not only lowered BFri by 51% (from 13.4 to 6.56) and BFro by 11% (from 32.92 to 29.16), but also decreased ΔHnom by 12% (from 108.4 to 95.9 mmHg, Table 7). Increasing ω from 10,122 to 10,313rpm (+2%) recovered ΔHnom back to 100.9 mmHg with an associated increase in maximum mean shear from 514 Pa to 520 Pa (+1.2%). At this speed, ηnom was 21.59% (up by 4.1% versus Model 2B). Velocity streamlines improved versus Model 1B, particularly between the rotor blades (Figure 5). Maximum attainable efficiency of ηmax = 24.07 occurred at Q@ηmax = 6.90L/min. Within the 4 < Q(L/min) < 8 range, ηmean was 22.57% (6.23% below ηmax and 4.54% above ηnom) and rotor drag coefficient varied within the range 0.14 < CDr < 0.158.
The computational design process of a novel axial-flow LVAD is presented with particular emphasis on the optimization of rotor geometry under steady-state flow conditions. Bearing design and inducer-rotor and rotor-diffuser gap hemocompatibility effects were omitted from the present work, while only superficial attention was given to magnet design and electromagnetic power feasibility.
Primary design objective was the maximization of rotor hydraulic efficiency (ηnom) at the nominal design point. To this end, an attempt was made not only to increase the maximum attainable efficiency (ηmax), but also to shift the entire efficiency curve along the flow axis so that the flow which maximized efficiency (Q@ηmax) coincides with nominal flow. It was also desired to have as flat a ηh -Q curve as possible so that a relatively wide range of flow values yield efficiencies close to ηmax. The latter condition particularly mattered since pump flow would be fluctuated for clinical purposes.7 Indeed, in Model 3A, ηmean was only 6.23% below ηmax, which occurred only 1.9L/min to the right of Qnom (Table 7).
Increasing blade trailing angle at rotor exit, βb2 an important parameter affecting pump performance in terms of head, torque, and flow, reduces the tangential component of exiting fluid velocity, hence the fluid attack angle (βf3) at diffuser inlet (Figure 2A). CFD results showed that increasing βb2 effectively decreased backflow and steadily increased ηmax (Table 4). Model 1B with βb2 = 32° was selected as its ηmax = 22.81% was highest among other geometries (Figure 3). Associated losses in ηnom (19.00 to 18.38%), head (107.2 to 101.8 mmHg), and torque (6.828 to 5.593 mN-m) were accepted in hopes of subsequent recovery by manipulating Ac.
At the inducer–rotor interface, visualization of fluid streamlines by CFD illustrated the relationship between rotor incidence, ir, and backflow (Figure 6). Velocity vectors (arrows) show forward (right) moving fluid around the inducer (horizontal on the left-hand side of each frame) toward the leading edge of the rotor blade (diagonal on the upper-mid region of each frame). As βb1 is reduced, rotor blade rotated clockwise and incidence increased. The boundary layer progressively separated from the inducer surface and percent backflow between two adjacent rotor blades increased (arrows pointing left, against the flow). Backflow was accentuated in the radial direction, from mid-plane (left column) to shroud (right column) as the blade radius, hence tangential velocity increased. Clearly, flow patterns were smoother in Model 2A (Figure 6, top row; I = 7°) than in Model 2B (Figure 6, second row; I = 14°), as confirmed by the 94% higher backflow (6.88 versus 13.4) in the latter model (Table 5). However, since ηmax and ηnom were highest with βb1 = 67°, Model 2B was selected, again, in anticipation of reduced backflow through subsequent optimization of Ac.
As expected, increasing rhub in Model 3A reduced BFri,nom and BFro,nom (Figure 5; Table 7) possibly because the hub-to-tip pressure gradient was lowered, better aligning velocity streamlines. But the associated reduction in Ac caused vx to increase, therefore βf1 to decrease (Figure 2A), which led to a slightly smaller rotor incidence. Considering data presented in Table 4, a reduction in rotor incidence angle may entail an increase in ηmax within that range (moving from Model 2B to 2A), although data are insufficient to make that interpolation. On the other hand, reducing Ac is always associated with a significant drop in ΔHnom (108.4 mmHg to 95.9 mmHg with 0.5mm increase of rh, Table 7). A less costly alternative would be to increase βb2, which forces Q@ηmax left, toward the design point (Equation 6); and the accompanying drop in ΔP would not be as much (Table 4) as it was when rh was increased (Table 7). Further, respective pressure losses can be compensated by a smaller increase in shaft speed, inducing a smaller rise in shroud shear.
It should also be remembered that overall efficiency can be optimized by modifying diffuser blade leading edge angle (i.e., manipulating diffuser incidence), which may prove a more practical design strategy for improving efficiency than modifying rotor deviation, an issue which will be addressed in subsequent studies.
This is the first of two reports on the virtual design and testing of a new Turkish axial-flow LVAD. Iterative optimization led to the geometry of Model 3A, which satisfied design criteria associated with flow patterns (hemolysis) and hemodynamics (physiology) within the size (anatomy) and power (electromagnetism) constraints. The focus of the subsequent report will be stator geometry optimization (including axial and radial gaps), while motor and bearing design will continue in parallel.
Results are presented with the understanding that all assessments should be repeated under transient conditions in order to have clinical bearing as the pump speed will be undulated. Also, pump geometry cannot be considered final until present computational findings are validated in physical (mock circuit, PIV) and biological (animal experiments and clinical studies) performance tests on physical prototypes.
From trigonometric relations (Figure 2A)
For rotor and diffuser blades, respectively
. Substituting into the efficiency formula
Taking the derivative of ηp with respect to Q
Executing the derivatives, rearranging and equating to zero, an 8th order polynomial in Q@ηmax is obtained
1. Frazier OH. Ventricular assist devices and total artificial hearts—A historical perspective. Cardiol Clin. 2003;21:1–13
2. Untaroiu A, Wood HG, Allaire PE, et al. Computational design and experimental testing of a novel axial flow LVAD. ASAIO J. 2005;51:702–710
3. Giersiepen M, Wurzinger LJ, Opitz R, Reul H:. Estimation of shear stress-related blood damage in heart valve prostheses—in vitro comparison of 25 aortic valves. Int J Artif Organs. 1990;13:300–306
4. Dixon SI, Hall CA. Two-dimensional cascades in, Fluid Mechanics and Thermodynamics of Turbomachinery, 6th ed. 2010 Burlington, MA: Elsevier
5. Turton RK. Fundamental principles in, Principles of Turbomachinery, 2nd ed. 1995 London, UK: Chapman & Hall
6. Hanselman DL Brushless Permanent Magnet Motor Design. 20062nd ed Lebanon, OH Magna Physics Publishing
7. Bourque K, Dague C, Farrar D, et al. In vivo
assessment of a rotary left ventricular assist device
-induced artificial pulse in the proximal and distal aorta. Artif Organs. 2006;30:638–642
Keywords:Copyright © 2013 by the American Society for Artificial Internal Organs
left ventricular assist device; axial-flow heart pump; computational fluid dynamics; velocity triangle; turbomachinery electro-magnetic efficiency; hydraulic efficiency; drag coefficient; incidence; secondary flow pathway; magnet cylinder