The Discrete-Continuous, Global Optimisation of an Axial Flow Blood Pump

This paper presents the results of the discrete-continuous optimisation of an axial flow blood pump. Differential evolution (DE) is used as a global optimisation method in order to localise the optimal solution in a relatively short time. The whole optimisation process is fully automated. This also applies to geometry modelling. Numerical simulations of the flow inside the pump are performed by means of the Reynolds-Average Navier-Stokes approach. All equations are discretised by means of the finite volume method, and the corresponding algebraic equation systems are solved by the open source software for CFD, namely OpenFOAM. Finally, the optimisation results are presented and discussed. The objective function to be maximised is simply pressure increase. The higher pressure increase the lower angular velocities required. This makes it possible to minimise the effect of haemolysis because it is mainly caused by high shear stresses which are related, among others, to angular velocities.


Introduction
Since advanced medical treatment is usually not enough to prevent the further decline of patients with heart failure, two treatments of patients with such a disease can be distinguished, namely heart transplantation and artificial heart blood pumps. The former is somewhat difficult because of the relatively high costs and, what is even more important, lack of donor organs [1], not to mention organ rejection and mortality rates. The latter approach to the heart failure problem, i.e. artificial heart pumps, have gained popularity due to constantly improved design and features.
Several designs of artificial hearts are known. The most complicated aim to replace the ailing heart whereas simpler pumps are designed in order to support it. The supporting devices fall into three groups: LVAD (the left ventricular assist device), RVAD (the right ventricular assist device) and BIVAD (the bi-ventricular assist device) [2]. Despite the fact that the human heart's nature is pulsatile, continuous flow pumps are the most popular solutions. This is because of the simplicity and size of the device. Furthermore, we can distinguish centrifugal and axial blood pumps. The latter being smaller in comparison to the former. However, axial blood pumps require significantly higher angular velocities in order to increase outlet pressure. This may lead to blood damage, i.e. thrombosis and haemolysis [3][4][5] in particular. Haemolysis is mainly caused by high shear stress [6]. In order to minimise the effect of haemolysis an optimisation process of the axial pump is undertaken. The optimisation results lead to improvement in the reduction of the wall shear stresses.
According to Behbahani et al. [6] the design optimisation process is a fully automated iterative improvement process, in which CFD simulation is coupled with an optimisation algorithm. What is more, this process is in contrast to a design approach. This is because the latter includes a human expert in the iteration loop.
Most published instances of CFD optimisation involve only regional optimisation, without parameterising the geometry of the whole pump [7]. What is more, in various published works, part of the optimisation remained a manual process [7]. For instance, Zhu et al. [8] optimised only the diffuser of an axial flow blood pump, with a fixed-shape straightener and a rotor by means of commercial CAD, CFD and optimisation solvers. Derakhshan et al. [9] investigated an optimisation of a centrifugal pump with only six geometry variables with the help of ANN and an Artificial Bee Colony algorithm. Zhang et al. [10] optimised, or in fact rationalised, several geometry parameters such as diameters, heights and blade shapes. However, no optimisation algorithm and method are reported, suggesting a trial and error approach rather than an optimisation process. Gouskov et al. [11] optimised a circulatory support pump with only four geometry variables by means of the little know LP-Tau sequence generator regarded as a global optimisation algorithm. The whole process was performed in three stages involving initial pump geometry, meaning that the part of the optimisation remains manual. Another method was developed by Hai et al. [7] for Archimedes screw rotary blood pumps with guide vanes by means of a commercial optimisation tool and a commercial CFD package. However, several simplifications were introduced such as a constant shaft diameter. Also, the head and tail of the shaft were missing. The optimisation problem involved only seven design variables. Frazier et al. [12] modified an existing axial-flow pump by increasing its inducer-impeller inlet angle hence increasing its pressure responsivity. Korakianitis et al. [13] tested the performance of sixty two axial pump impellers with varying outlet angles and number of blades, i.e. only two design variables. Furthermore, the gathered data allowed the estimation of the optimal axial impeller geometry for any desired operating condition. The optimisation (rationalisation) method appears to be generate-and-test or the so called exhaustive search.
In this paper a fully automated iterative improvement process is discussed, coupling CFD simulations with an efficient global optimisation algorithm. The optimisation process is carried out by means of a slightly modified Differential Evolution (DE) algorithm, which is regarded as one of the current state-of-the-art algorithms [14] and typically forms the basis of the best performing algorithms. A simple and effective method of geometry modelling is proposed parameterising the geometry of the whole pump by means of fourteen variables. There is no need for any intermediate CAD software since the geometry is created directly by means of freely available GNU Octave and Octave Geometry scripts [15]. The same concerns the optimisation algorithms [16] available for everyone. A steady state solution is also compared with transient solutions in terms of pressure increases for various time steps and two different methods.

Geometry Description
The whole geometry consists of a rotor and stator and is described by fourteen parameters (design variables) x = {x 1 , . . . , x 14 } listed in Table 1. First of all, the rotor blade shape is a non-linear helix given by where R stands for the radius of rotor blades. Further, the parameter t ∈ [0; 2πx 1 ] is related to the rotor blade pitch 2πx 1 . The non-linearity f (t) of the helix Eq. 1a is described by one parameter x 2 shown in Fig. 1. Two of fourteen design variables are discrete, i.e. the number of rotor and stator blades. The former is named x 3 and the latter x 8 . Next, the shape of the shaft is described by four points (x 4 , x 5 ), (x 6 , x 7 ) of the spline shown in Fig. 1. Stator blades are given by the so called camber line and the blade thickness. The former is represented by two variables (one point) (x 11 , x 12 ) of the spline whereas the latter by two variables (two points) (0, x 13 ), (1, x 14 ), see Fig. 2. Finally, the two remaining variables x 9 and x 10 are the stator blade upper and lower twist angles, respectively. The considered rotor radius R is 8 mm and the length of the shaft is 4.5 R.
The whole optimisation process is fully automated. This also applies to geometry modelling. Several random example geometries are shown in Fig. 3. What is important is that the proposed method is simple and effective. There is no need for any intermediate CAD software since the geometry is created directly by means of GNU Octave and Octave 'Geometry' script in a semi-discrete form (STL format).

Governing equations
Although blood is a suspension of blood cells in the plasma [17], the blood flow in large vessels as well as blood pumps can be regarded as a single-component and single-phase fluid, i.e. blood may be treated as a Newtonian fluid. This method is suitable for vessels larger than 0.1mm. If, however, the diameters are in the range of 0.1 to 1 mm, which is not the case here, then non-Newtonian constitutive equations are required in order to account for the non-Newtonian phenomena, such as shear thinning, yield stress and constant viscosity values at high shear rates. Numerical simulation of the flow inside the pump is performed by means of the Reynolds-Average Navier-Stokes approach. A closed system of equations [18] for incompressible fluid, involving a mass conservation equation, the Reynolds equation and two additional transport equations for the two-equation SST model, is solved. An additional equation for the eddy viscosity ν t is also necessary together with two blending functions F 1 and F 2 [19] ∇ ·ū = 0, In the above equations u is the velocity vector, p -pressure, ρ -density, ν -kinematic viscosity coefficient and D -strain rate tensor. The two additional transport quantities are the kinetic energy of velocity fluctuations k and the turbulence frequency ω. The shear stress transport (SST) model combines the k-ω model near the wall with the k-ε far from it. Constants marked with the subscript '3', namely σ k3 , σ ω3 , α 3 , β 3 are linear combinations of constants from the component models, i.e.

Equation discretisation
All the equations are discretised by means of the finite volume method, and the corresponding algebraic equation systems are solved by the open source software for CFD, namely OpenFOAM [20]. Since all the transport Eq. 2a-2e have common terms, the general transport equation for a quantity φ has the form of where the overall source term S φ should be linearised first, if necessary S φ = S u + S p φ.
The diffusivity for φ is denoted as Γ .
The integral version of the above transport equation over a control volume V P can be expressed as where |V P | is a measure of V P and S f is a surface normal vector pointing outward. The considered convex control volume V P around a centroid P consists of f planar surfaces S f . Divergence schemes include both convection terms ∇ · (φu) and other diffusive ∇ · (Γ ∇φ) terms and involve Gauss integration. The discretised convection term needs to be interpolated by means of cell centred values because the values φ f are located at the face centroids. Limited linear interpolation is used (the steady state case) or linear upwind (the transient case), both being second order accurate. Further, the discretised diffusive terms involve surface normal gradients (∇φ) f · S f , and are evaluated at a cell face that connects two cells. In order to maintain second order accuracy for non-orthogonal meshes, apart from orthogonal schemes, a non-orthogonal correction is considered.
The SIMPLE algorithm is used in order to solve pressure-velocity coupling, and the pressure equation is solved by means of the GAMG solver with the DIC smoother for the steady state version of the system Eq. 2a-2e. For the velocity fields and turbulent quantities standard solvers using a GS smoother are utilised. Under-relaxation factors are used in order to improve the stability of a solution. This is particularly important when solving steady-state flows. The assumed factors are 0.3 for pressure, 0.7 for velocity and 0.5 for the turbulent quantities k and ω. The transient version of the system is solved by means of the PISO [21] or PIMPLE (PISO + SIMPLE) algorithm, discussed further in Section 3.3.
The integrand dφ P dt of the left hand side of Eq. 4 is discretised by means of an implicit multi-level scheme What is more, this method is known to be second order accurate in time. Additionally, this approach is the so called three-level method because it requires the values of the unknown function φ P at three different time steps, namely φ n+1 P , φ n P and φ n−1 P .

Space and temporal discretisation
The flow domain is divided into three parts, see Fig. 4 (bottom). Apart from the rotating rotor, two additional steady pipes (cannulae) are considered. The cannula shape optimisation problem is a separate problem discussed elsewhere [14]. Furthermore, all three domains are discretised separately and merged by the so called arbitrary mesh interface (AMI). The total number of nodes is 1 689 534 and the total number of volumes is 1 495 877 where 1 381 669 of them are hexahedra. Formally, the mesh used may be classified as Cartesian, see Fig. 5 (left).
In order to make certain that the flow near the walls is properly resolved thin layers around the physical walls are generated. The quality of the mesh near the walls is inspected in terms of the maximal y + values. Maximal values of y + are below 1 for all the considered walls including the blades, see Fig. 5 (right). It has to be clarified, however, that the Open-Foam implementation of the k-ω family models (the SST among them) in the near wall region allows for a scalable wall function if 1 < y + < 300 or no wall function if y + < 6. Two options are then possible and two were inspected giving negligible differences in terms of pressure rises. However, the former appears to be more stable which is crucial when it comes to a fully automated optimisation process. Figure 6 demonstrates mesh convergence by showing the influence of the number of nodes on the pressure increase p ρ −1 . It is obvious that increasing the number of the mesh node above 2×10 6 has negligible effects on the pressure increase. This is because results are nearly constant above 2×10 6 . Given that the mesh size is crucial from the CFD calculations perspective and thus optimisation time, the corresponding number of nodes was chosen at the level of 1.7 × 10 6 .
A steady state solution (MRF) in terms of pressure increase is next compared with transient solutions for various time steps. Two methods are considered, namely PISO (Pressure-Implicit with Splitting of Operators) and PIMPLE.
The number of time steps corresponding to the CFD time step is N . For instance, if the time step corresponds to the 4 • angle of revolution then N = 90 and so on. A time step convergence for the considered pump for five revolutions and two different approaches (PISO and PIMPLE) is shown in Fig. 7. As for the PISO algorithm five different Fig. 7 Temporal convergence cases are presented namely 2 • , 1 • , 0.5 • , 0.2 • and 0.1 • time-step. Furthermore, the moving averages (the thick solid lines) of the pressure increase Eq. 7 are super-imposed. It is well visible that the PISO algorithm over-predicts the pressure rise in comparison to the steady solution (MRF -the solid black line). Decreasing the angle of the revolution below 0.2 • has no effects on the pressure increase. In the second case, i.e. PIMPLE, the situation is completely different. It is evident that decreasing the temporal resolution from 1 • to 6 • per time-step has a negligible effect on the results in terms of pressure increase. What is more, the PIMPLE algorithm gives results similar to a steady solution (MRF). Thus, it can be concluded that stationary solutions are representative. This is exceptionally important from the point of view of the optimisation algorithm. This is because such CFD calculations must be performed repeatedly. Finally, the transient simulations require typically 3 full revolutions to reach the pseudo-periodic state.

Boundary conditions
The main boundary conditions are: -Inlet. The specified constant volumetric flow rateV = 3 dm/min is directed perpendicularly to the inlet surface accompanied by the zero normal gradient pressure. Low turbulence intensity is also considered. This means that the turbulence intensity is at the level of 1% and viscosity ratio ν t /ν = 1 where ν = 3.3019 × 10 −6 m 2 s −1 and the reference density ρ = 1060 kg m −3 . -Outlet. The constant pressure distribution is assumed here together with zero gradient velocity for the flow out of the domain. This is because the outlet surface is located relatively far from the rotor. -Walls. The so called no-slip conditions is applied meaning that impermeability and adhesion requirements are forced. Rotating wall velocity n = 6000 rev/min is considered in the rotating frame of reference. The flow in the near wall region is modelled by means of the scalable wall function. -Interfaces. In order to allow for coupling between stationary pipes and rotating part of the pump the so called cyclic arbitrary mesh interfaces (AMI) are considered. A steady state approach is used rather than a full transient rotor-stator interaction as explained in Section 3.3. This is crucial for time consuming optimisation processes since CFD calculations need to be repeated hundreds of times. This approach is commonly known as a multiple reference frame (MRF) simulation.

General remarks
A discrete-continuous optimisation approach is needed to find the optimal shape. This is because two of fourteen design variables are discrete, namely the number of rotor and stator blades. The remaining variables such as the shaft, rotor and stator shapes as well as their thickness are continuous. Metaheuristic procedures or more precisely Differential Evolution (DE) are used as a global optimisation method to localise the optimal solution in a relatively short time [16,23]. Another commonly used term is nature-inspired metaheuristic. DE can be further classified as multi-point (population based), derivative free optimisation algorithms. Most importantly, no additional information about the objective function is required. This is because DE is not problem-specific, stochastic algorithms with randomisation and local search. Additionally, randomisation is introduced through the probability of crossover. This makes it possible to efficiently explore the design space and escape local minima.

Objective function
The optimisation problem is to find a maximal pressure increase p. Since most optimisation algorithms are designed for the minimisation of the objective function, the pressure increase is considered with a minus sign Furthermore, the argument of the global minimum value of the objective function is expressed as where D = 14 stands for the dimension of constraint space Ω or simply the so called optimisation domain. Thus, the objective function is subjected to box constraints listed in Table 1 Ω where L i and U i are lower and upper bounds, respectively. Box constraints are regarded as a special case of inequality constraints. This type of constraint is commonly met in optimisation problems and does not need any special treatment.

Algorithm
Differential Evolution [24] is a simple, fast and effective metaheuristic algorithm. Furthermore, DE is regarded as the next step in evolution of the Genetic Algorithms (GA). Crossover and mutation are utilised on floating-point vectors. Additionally, selection is also present in DE. Most importantly, an explicit update equation is provided in contrast with GA. Four main features of the algorithm can be distinguished, namely three different individuals selection, mutation, crossover and selection. Once three randomly individuals x a 1 , x a 2 , x a 3 are selected, a mutant vector v i is generated according to the so called DE/Rand/1/Bin The scale factor F , or the so called differential weight, is used in order to control the rate of population development and is assumed here to be F = 0.7. Furthermore, the trial vector y i is created via binomial crossover with probability C = 0.9 (line 10 in Fig. 8). The crossover probability regulates how much of the mutant vector is copied to the trial vector. It is possible to combine mutation and binomial crossover in a single vector equation by means of the Heaviside step (theta) function H. An auxiliary D-dimensional vector K consists of 0 and 1 (line 7 in Fig. 8) where U (0, 1) stands for the continuous uniform distribution realisations characterised by its minimum 0 and maximum value 1. This approach is somewhat different in comparison to the original algorithm [23,24]. What is more, it is now possible to combine mutation Fig. 8 The vectorised differential evolution pseudocode and crossover, into a single vector formula -line 10 in Fig. 8. Three different and randomly chosen individuals are indexed on the basis of a random permutation vector a (line 9 in Fig. 8). Additionally, line 8 corresponds to setting a random index of the vector K to 1 in order to guarantee that the y i = x n i . Here U {0, D − 1} represents the discrete uniform distribution realisations where 0 and D − 1 are the parameters of the distribution. The selection step is shown in line 12 in Fig. 8. Simply, the best solution of the trial vector y i and original individual x n i , in terms of the objective function value, is passed onto a next generation x n+1 i . What is interesting is that this step is fully deterministic in contrast with mutation and crossover. The algorithm terminates if a given stop criterion is satisfied, i.e. the maximum number of objective function evaluations. Finally, a uniform random distribution is generated within the search domain with a random seed based on the clock (line 2 in  Table 1 (the second column).
The presented algorithm [16] is almost identical in comparison to the original algorithm [24]. The differences include a vector representation, which is executed faster by mathematical packages such as GNU Octave. Furthermore, the presented version of the algorithm includes an option to check the range of variables (line 11 in Fig. 8), which allows to avoid non-physical configurations such as a negative number of blades.

Results
The  different population sizes N and iteration numbers n max are shown in Table 2. The best solutions were obtained for the 20 × 20 configuration of DE. The individual computing time for calculating one objective function, i.e. geometry modelling, discretisation and CFD calculation, is about 15 minutes on a i7-6850K 3.60 GHz processor (3 out of 6 cores involved). For example, the calculation for the entire optimization process with a population size of N = 20 and a number of iterations of n max = 20, i.e. 400 evaluations, takes about 4 days.
From Fig. 9 and Table 2 it arises that the best results are obtained for medium populations with the size N of 20 individuals. Furthermore, increasing the number of iterations n max above 20 does not bring further improvement. This applies to both population sizes, i.e. n max = 20 and n max = 15. Populations with larger sizes, i.e. N = 30, require more iterations to obtain results comparable to smaller iteration numbers n max . This is important due to the time-consuming CFD calculations. Out of all the optimisation processes, the best results are obtained for N = 20 and n max = 20. Due to the stochastic nature of the DE algorithm, it is obvious that restarting the algorithm may give slightly different results, which, however, are typically close to the optimal solution.
The optimal values of g are listed in Table 1. The upper part of Fig. 4 presents the optimal pump shape according to Table 1. It may be observed that the optimal geometry consists of 4 rotor and 8 stator blades being the upper accessible ranges. As for the remaining variables, most of them are located on the boundary of the considered box constraints which is a typical situation during constrained optimisation. It is worth mentioning that some constraints in Table 1 allow us to avoid non-physical configurations. Other constraints require preliminary optimisations or at least geometry creation and check.
The average shape evolution for the selected generations (iterations) for the 20 × 20 DE, reflecting convergence process, is presented in Fig. 10. It is interesting that the first shape (iteration 1) is simply an arithmetical average of purely random shapes. This is because the initial population is randomly generated. At the same time, the last shape (iteration 20) is similar to the optimal shape since the population here is nearly uniform. Finally, since these shapes are arithmetical averages of individual generations, it should be noted that none of these has been subject to any CFD calculations.
The upper part of Fig. 11 displays the best shape evolution resulting from the shape optimisation process (N = 20, n max = 20). The last shape (iteration 20) is the optimal solution shown also in Fig. 4. Interestingly, all shapes apart from the second (iteration 2) are similar in terms of rotor configurations. What distinguishes the second geometry from the remaining three is the shape of the shaft and stator blades in particular. The lower part of Fig. 11 presents the corresponding wall shear stresses (WSS) distributions. The higher the iteration number the smoother the distribution on the shaft and blades. Also, the maximum values of WSS are lower. One has to keep in mind, however, that high shear stresses are related primarily to angular velocities.

Conclusions
-The results of a discrete-continuous global optimisation process of an axial flow blood pump are presented. The optimal results are achieved by means of DE in a relatively short time.  Figure 11 presents the wall shear stress distribution which is responsible for haemolysis.
It may be observed that the highest values are localised, among other, on the rotor and stator blade tips. The higher pressure increase the lower angular velocity required. Thus the optimisation process leads to improvement in the reduction of the wall shear stresses and the effect of haemolysis. This is because haemolysis is mainly caused by high shear stresses which are related, among others, to angular velocities.

Limits of the current study
It should be mentioned that many ventricular assist devices designs are also equipped with guide vanes (flow straighteners) not only at the outlet, but also at the inlet of the device. This kind of flow straightener is missing in the current study. However, taking into account the presence of a straightener would not change the algorithm itself, but only extend the calculation time. The same concerns the gap sizes between rotor and housing. One has to also keep in mind that the current optimisations will, most likely, not be sufficient to develop an efficient and biocompatible blood pump. It is well known that the wall shear stress is responsible for the phenomenon of haemolysis. In addition to haemolysis, another phenomenon is exposure time. Considering the above criteria, if at all possible, could lead to other, less effective solutions than those presented in this paper. One possible solution is multi-objective optimisation that would take into account both pressure increases, haemolysis and exposure time. The first problem is the much longer multi-objective optimisation time in comparison with the single-objective problem. The main problem is the availability of a CFD model allowing an accurate prediction of haemolysis in such complex flows (stagnations, recirculations). This problem is constantly being analysed and works devoted to it are still being published [7,25,26]. According to Hai et al. [26] there is no clear guideline concerning the applicability or the limitations of the studied models. What is more, a greater problem is the lack of well documented experimental data, adequate for a validation of the discussed models [26]. It seems that one of the possible solutions to this problem is an experimental stand that is currently being prepared. The experiments carried out together with biologists will most likely answer the above questions. This issue will be addressed at length in a separate work.