Sebastian Hallier
email: sebastian.hallier@uni-rostock.de

Benjamin Torner

Jitendra Kumar

Frank-Hendrik Wurm

Miniaturization of an implantable pump for heart support

Keywords:

The use of implantable pumps for the heart support has proven to be a promising option for the treatment of advanced heart failure. Avoiding hemolysis and achieving high efficiency rates represent two of the main challenges in the miniaturization process. For the miniaturized pump configuration, the outlet guide vane is replaced by a volute. A method, based on the concept of constant angular momentum, is developed to facilitate the design process. The index for global prediction of the hemolysis MIH is modified for a local application, which enables to locate and understand the sources of hemolysis. A reduction of the pump's dimensions by more than 60 % is achieved, while improving the hydraulic efficiency and reducing hemolysis rates. An optimization of the impeller airfoils is still in progress.

INTRODUCTION

Heart failure is a form of cardiac disease in which the heart's pumping power is weaker than normal. This chronic disease is spread over twenty million people worldwide, out of that one million people need heart transplantation because of their advanced heart failure. In contrast, only 3000 people are donating hearts each year, which urges us to develop a technical solution as a substitute. Technical solutions can be distinguished as, (a) total artificial heart in place of natural heart, (b) heart support system called as Ventricular Assist Devices (VADs). Both solutions can be realized as intracorporeal (inside the body) or extracorporeal (outside the body) solutions. The use of intracorporeal solutions brings much higher life quality for the patient. The patient can leave the hospital and enjoy a "normal" life with work and leisure activities. The application of an artificial heart or a VAD could be following: The development status of total artificial hearts is very low. Only few applications are reported and the remaining life expectancy of the patients is relatively small. In the majority of heart failures, the left ventricle is damaged primarily. Therefore the main focus of research is dedicated on VADs for the left ventricle. The current research paper deals with an intracorporeal VAD for the application b), c) and d). The following requirements for the applications of heart support system can be defined:

• flowrate Q = 2,5 -10 l/min ; head H = 80 mmHg • high efficiency (very low heat and temperature rise tolerance for organism) • lowest damage of erythrocytes and thrombocytes in the blood inside the support system (erythrocyte damage has the higher importance) • low noise level (patients and their families feel very disturbed by noise emission) • life time longer than 2 years at least with low MTBFvalues • small dimensions • easy to implant (safe and easy attachment to the heart) • costs as low as possible A promising concept is the realization of a VAD as an axial flow pump. It has the advantage that the drive can be integrated into the pump easily and the most affected left ventricle can be supported by this concept optimally.The axial flow pump consists of the hydraulic part of the pump and the drive where the rotor is placed within a magnetic bearing. These pumps are available today and have already been implanted in a couple of thousands patients. It is usually placed below the heart as shown in Figure 1. The biggest disadvantage of the available solution is the relatively large dimension which restrict to use it. These VADs cannot be implanted in children and the use in women is limited because of the required space inside the body. The main target of our research work is the drastic miniaturization of the hydraulic part with respect to all of the above mentioned requirements. A simple way to reduce the dimensions of a pump is to increase the rotational speed. Unfortunately, this leads to non-acceptable damage of erythrocytes. The damaging mechanisms of the erythrocytes are complex. Very high shear rates for a short exposure time or lower shear rates for longer time, both can rupture the blood cells membrane, which is called hemolysis. Therefore, one of the leading parameter in the development of a VAD is to reduce the hemolysis inside the pump. Miniaturization of the pump and reduction of hemolysis are focused as main objectives in the current research work.

METHODS

Design of the pump

The requirements on the pump and the available space lead to an unusual design of the pump to some extent. The hydraulic part of the existing pump consists of an inlet guide vane, an impeller and an outlet guide vane. Miniaturization requires the modification of the pump, which leads to the idea to replace the outlet guide vanes by a volute with a radial discharge line. This configuration allows the attachment of the VAD in the direct vicinity of the heart. The suction line can be implanted directly into the heart and the discharge line is located parallel to the apex of the heart. The miniaturization gives the possibility to implant by using minimal-invasive methods.

For the impeller design two different methods are used, a semi-empirical method with consideration of cascade effects [START_REF] Wolf | Ein einfaches Berechnungsverfahren für Verdichtergitter[END_REF] and a method of singularities. Both methods are tested, in order to find the best method for the required pump data and small dimensions which lead to an unusually high …… ratio between the gap (impeller to impeller casing) and impeller diameter. For the design of the casing a method based ontthe principle of constant angular momentum is developed and used. One of the main challenges is to redirect the axial impeller outflow into a tangential outflow graft. Spiral volutes can be designed with various different cross-sectional shapes. For this specific application the shape is defined by splines in order to facilitate the most efficient use of available space and to leave many degrees of freedom to the designer. The result is an asymmetric volute that expands in radial and axial direction.

Fluid mechanical optimization

The investigation and optimization of the flow field is done by using numerical methods. Three-dimensional numerical models of the entire pump are created in order to perform URANSE simulations at different load conditions. For the investigation and optimization the commercial software package ANSYS CFX is used as a solver. To perform transient simulations, the k-ω-SST model is selected as it combines the advantages of the k-ε and the k-ω model. Curvature correction and the Gamma Transition model are used in order to account for the strongly curved streamlines and the low Reynolds numbers inside the pump respectively. Blood has non-Newtonian properties up to a shear strain rate of approximately 100 s -1 . At higher shear strain rates Newtonian properties can be assumed [START_REF] Oertel Jr | Bioströmungsmechanik. Vieweg+Teubner[END_REF]. The shear strain rates in all components of the pump are much higher than 100 s -1 so Newtonian properties with a constant dynamic viscosity of 0,0035 mPas is assumed for the blood. All simulations are performed on structured meshes with maximum values for y + aimed to be around 1 in order to fully resolve the boundary layer.

Optimization regarding hemolysis

As mentioned above, one of the main assessment criterion of the design of VADs is a limited hemolysis generation. Hemolysis is defined as the destruction of red blood cells (erythrocytes) under the influence of high shear rates. As a result of the cell disruption (lysis), the hemoglobin is released into the blood plasma irreversibly, whereby the erythrocytes lose their ability to bind oxygen [START_REF] Oertel Jr | Bioströmungsmechanik. Vieweg+Teubner[END_REF].

The standard way for the numerical hemolysis prediction is the use of a global criterion. In this project, the standard way is extended by use of a global and a local hemolysis prediction model. Both models are based on an empirical power-law approach, which takes into account the parameters involved in the lysis, shear stress and time of exposure to these stresses [START_REF] Oertel Jr | Bioströmungsmechanik. Vieweg+Teubner[END_REF]. The power-law correlation of Giersiepen et al. [START_REF] Giersiepen | Estimation of shear stress-related blood damage in heart valve protheses-in vitro comparison of 23 aortic valves[END_REF] is used, which rests upon the measurement of the ratio of plasma-free hemoglobin (∆ ) to the total hemoglobin concentration ( ) in the blood (Equation ( 1)).

= ∆ = 3.62 ⋅ 10 ⋅ . ⋅ . (1)
For the global hemolysis prediction model, the quantity Modified Index of Hemolysis MIH is derived from the damage fraction by means of = , ⋅ 10 . The MIH is an often used quantity for evaluating blood damage in medical devices. Following the idea of Garon and Farinas [START_REF] Garon | Fast Three-dimensional Numerical Hemolysis Approximation[END_REF], a time averaged, global MIH value is defined through Equation ( 2).

! = " 1 # $ % &'( ) * . ⋅ 10 ( 2
)
where % is the rate of hemolysis production per unit time. This quantity can be depicted with respect to the power law by Equation ( 3).

% = 3.62 ⋅ 10 . ⁄ ⋅ . . ⁄ (3)
The correctness of the constants and exponents in Equation ( 1) and ( 3) has been called into question by research workers. This is due to the fact that the power law together with these exponents overrates the hemolysis to a significant degree. The reason for this error was additional damage to the blood cells due to the seal in the experimental device [START_REF] Taskin | Computational Characterization of Flow and Hemolytic Performance of the UltraMag Blood Pump for Circulatory Support[END_REF]. In spite of this fact, these constants and exponents have been widely used by many researchers.

In complex medical applications it is common to replace the shear stress tensor by a scalar representation [START_REF] Pauli | Transient Stress-and Strain-Based Hemolysis Estimation in a Simplified Blood Pump[END_REF]. In the present paper the scalar shear stress is calculated using the eddy viscosity, dynamic viscosity and the shear strain rate (Equation ( 4)).

= , + , . ⋅ /0 = , + , . ⋅ 12 23 4 25 6 7 46

Through Equation ( 2), ( 3) and ( 4) and the time averaged quantities of a numerical simulation, it is possible to determine a single value, which can be used to estimate and classify the overall hemolysis behavior of a blood pump. However, local hemolysis information cannot be derived from ! .

The second, modified way of hemolysis assessment allows to predict hemolysis locally. This opens the possibility to pinpoint the location of high lysis inside the computational domain of the medical device [START_REF] Garon | Fast Three-dimensional Numerical Hemolysis Approximation[END_REF]. The model based upon solving a scalar transport equation in the computational domain. The solved scalar quantity is the linear damage fraction 8 . This parameter can be derived from the relation 8 = √ :.;<= . A hyperbolic partial differential Equation ( 5) is obtained through this definition. This equation is implemented in ANSYS CFX by use of the option Additional Variable [START_REF]CFX Documentation[END_REF].

' ' 8 = > 2 2 + 3 ⋅ ?@ 8 = % (5)
The source term % of this equation is the same quantity as in Equation ( 3). In post treatment, the linear damage fraction is transformed into local MIH values A = 8 .

⋅ 10 for the hemolysis estimation in the medical device.

A connection between the local and global model can be found through Equation [START_REF] Pauli | Transient Stress-and Strain-Based Hemolysis Estimation in a Simplified Blood Pump[END_REF]. A verification of the local hemolysis criterion is feasible using this equation.

! = 10 ⋅ " 1 # $ B3 ⋅ CD 8 '7 E * . = 10 ⋅ " 1 # $ % &'( ) * . (6)
As mentioned above, the hemolysis in blood pumps is overestimated when ! is calculated using Equation ( 2) or [START_REF] Pauli | Transient Stress-and Strain-Based Hemolysis Estimation in a Simplified Blood Pump[END_REF]. Because of the overevaluation, its not reasonable to compare the numerically determined, overall MIH value with tolerable hemolysis indexes for blood pumps from experiments (e.g. in [START_REF] Franco | Advanced Therapy in Cardiac Surgery[END_REF]). However the numerically computed, global MIH value can be used to qualitatively compare the overall hemolysis behavior of different blood pumps. So it provides a means to classify medical assistant devices and is therefore a valuable engineering tool [START_REF] Garon | Fast Three-dimensional Numerical Hemolysis Approximation[END_REF]. Since the local hemolysis model also uses the Equation (3) as a source term, hemolysis is overestimated with this approach as well. Nonetheless, it is possible to identifiy and highlight regions with high lysis qualitatively.

Use of the modified method for hemolysis prediction

For the prediction of local hemolysis, additional transient numerical simulations are performed for the new impeller at nominal load. An additional transport equation is solved where the calculated scalar quantity is the linear damage fraction 8 . As initial condition, the damage fraction is assumed to be zero in the computational domain at time = 0 and similar to Taskin et al. [START_REF] Taskin | Computational Characterization of Flow and Hemolytic Performance of the UltraMag Blood Pump for Circulatory Support[END_REF], the transported quantity at the inlet and the scalar flux through the walls are set to be zero as boundary conditions. The calculation is executed over 40 impeller revolutions with time averaging for the last three revolutions. This long time period is due to the initial condition, because the blood damage by hemolysis must develop from zero to a statistically steady state. During the last revolutions, attention is paid to that the term F B3 ⋅ CD 8 '7

E

shows a periodical behavior at the outlet of the rotational domain to guarantee reasonable time averaged values.

Optimization of the airfoil-shape regarding hemolysis at varying inflow conditions

Airfoils produce the lowest flow losses when the vector of the velocity meets the airfoil tangential to the camber line. The losses, and in our application the hemolysis, can increase significantly at varying inflow angles. The operating condition of VADs vary permanently with multiple factors like the remaining activity of the failed heart or the physical stress of the patient. In order to ensure low rates of hemolysis in the broadest possible range of operating conditions, an optimization process is in progress. A pulsating flow is realized as inflow condition for the optimization.

A full cardiac cycle takes around 0.8 seconds and the impeller of an axial VAD is rotated by several thousand revolutions per minute. Therefore more than one hundred revolutions of the impeller about its axis are simulated for a single cardiac cycle. In order to test different optimization algorithms the optimization problem is transferred into a 2D setup with an airfoil in a relative frame of reference as the first step.

The geometrical model consists of an airfoil section at a given radius of the impeller. Cascade effects are accounted for by the angular spacing between translational periodicities (Equation ( 7)). Assuming swirl-free flow at the inlet and constant revolutions per minute, the velocity components are derived from the velocity triangle. For the relative frame of reference the circumferential component G H is the reversed circumferential velocity and is constant over time. Since the flow rate varies over the cardiac cycle, the meridional velocity is time dependent as shown in Equation ( 8). This is realized using a discrete Fourier transformation (DFT) with a length of series of three (Equation ( 9)). The result is a pulsating meridional velocity, corresponding to a flow rate that varies from 3.2 to 6.1 l/min over a period of 0.79 s. The resulting inflow condition is shown in Figure 2. Radial velocity components are neglected.

G M = N M = # O (8)
P = Q + R S 4 NTI >2JU V @ + W 4 IUC >2JU V @ X 4Y (9)
The airfoil geometry is defined by a constant stagger angle and camber line which is superposed by a variable but symmetric thickness distribution. For the parametrization the bi-superellipse-profile (BSP) method by Bross [START_REF] Bross | Entwicklung neuer Schaufelgitter aus Profilen variable Geometrien zum Einsatz in Leiträdern drallgeregelter Turbomaschinen -Teil II[END_REF] is used. With this method the shape of an airfoil can be defined by a set of parameters just like NACA airfoils are described by their digits. The algorithm for describing the airfoil is based on the superellipse-equation (10).

Z [ W \ M + Z 5 S \ ] = 1 (10)
The camber line is described by a circular arc. The upper and lower surface of the airfoil are each defined by three separate super-ellipse curves for the front, middle and back part. This way the thickness distribution of the upper and lower surface could be described independently from one another but doing lso would alter the blade angles at the leading and trailing edge. This would mean that the circular camber arc does not coincide with the aerodynamic camber. In order to maintain the same pressure increase at BEP the blade angles must not be varied. For this reason the upper and lower surface of the airfoils in this study is always modified symmetrically to the camber line. A visualization of the BSP geometry definition is shown in Figure 3. Variable parameters for the optimization are the maximum thickness dmax, the location of maximum thickness xd and three more parameters n1, n2 and m2 which define the slope of the airfoil. All other parameters are kept constant.

A genetic algorithm from the MATLAB Global Optimization Toolbox is used in the first to create the populations and evaluate the fitness of the individuals because of its ability to find a global optimum. The maximum value for the wall shear stress that occurs on the airfoil over a full cardiac cycle is defined as the fitness function since high shear stresses are a main cause for hemolysis. In addition, low wall shear stresses indicate a smooth guidance around the airfoil in off-design conditions.

The transient simulations have been carried out with a time step size of 5•10 -5 s so that a full cardiac cycle is completed within 1580 time steps. In addition 200 initial time steps are given for each simulation to settle. These are ………... ………………………………………………..

Q [l/min] t [s]
disregarded for the evaluation of the fitness function. The optimization by using BSB-profiles is still in progress and results will be shown in future studies.

RESULTS AND DISCUSSION

The optimization process of the new components lead to a miniaturized pump with dimensions which are 60% less than those of the existing product, while being designed for the same flow rate. This will allow a much broader use of this new VAD. The placement of the new VAD compared to the existing one is shown in Figure 4. It can be seen, that the required space within the body is drastically reduced. The VAD could be implanted by using minimal invasive methods.

Despite the miniaturization, the new pump design outperforms the existing design in terms of hydraulic efficiency as shown in Figure 5. Hydraulic efficiency is very important in VADs since the dissipated losses lead to heating of the fluid and structures which may result in clotting of the blood. One of the reasons of the improved efficiency is the shortened cannulas that result from the intraventricular placement of the pump. Here and in all following diagrams, the operating points for both designs are the same and cover the normal range of an average adult's cardiac output. This excludes the upper of required flow rates, since these rates usually occur for very short periods of time.

N ^= _ . .,`M& &&&&&&&&& a _ . . &&&&& b 2 3 (11)
The plot of a normalized pressure loss coefficient, defined as in Equation (11), is shown in Figure 7 for both designs at nominal load. Even though the flow from the impeller approaches the guide vane blades without incidence, wide regions of pressure losses can be observed within the first part of the vane passage. This is due to a non-uniform flow field. Furthermore, there is a separation of the flow nearby the trailing edge, which causes mixing losses in the wake region. Inside the volute, there is a small zone of high losses right behind the volute tongue. It is caused by a minor separation zone. The losses increase over the volute's circumference and into the diffuser but remain on a much lower level ………… In Figure 6, the pressure head curves for the existing and new impeller are shown at the operating condotions where both pumps deliver the required normalized pressure head of 1 at Q/Qopt=1. As a result of the higher efficiency and the lower losses inside the volute, the impeller needs to provide a lower pressure head. The flatter progression of the pressure head curve is required by medical reasons and is realized in the impeller design. In Figure 8, the curves for the global hemolysis index ! of both designs are shown. While the values are nearly the same at BEP, the new pump shows significantly lower rates of hemolysis at both part load and overload conditions, despite the miniaturization of the pump. The increase in efficiency and decrease of hemolysis are caused by lower pressure losses inside the volute. Another explanation for the lower rates of hemolysis is the lower and flatter pressure-head curve of the new impeller.

………………. …………

The theoretical specific work c .d developed by the impeller is defined in Equation ( 12) where 3 N H disappears for swirlfree inlet flow. In that case the generated pressure head is the product of the impellers rotational velocity u2 and the fluids circumferential velocity component cu. High tangential velocities induce high shear stresses inside the clearance gap while the circumferential velocity leads to increased shear stress behind the impeller due to higher fluid velocities. This correlation is supported by the increasing gap in hemolysis rates and in the pressure-head curves of both designs towards part load and overload conditions. The local assessment of hemolysis allows the understanding of the main reasons for hemolysis. Local hemolysis is analyzed for the impeller. In a first step, the local hemolysis prediction model is verified. The averaged global MIH value is calculated on the basis of Equation ( 6) and the term B3 ⋅ CD 8 in post treatment. It shows a relative deviation of less than two percent from the MIH value determined with Equation (2). So, the local model is able to reproduce a global MIH value, which is comparable to the value of the global prediction model. Another source of hemolysis in that region is the leakage flow path between the rotor blade tip area and the impeller casing. To analyze this, the time averaged source term % of Equation ( 5) is shown on the impeller surface in Figure 10. This term describes the average rate of hemolysis production. High values between % ≥ 0.3 … 0.5 can be observed at the whole blade tip area. From this follows that the plasma-free hemoglobin concentration is increased (the erythrocytes are damaged because of shear stresses and release their hemointo the plasma) when the blood flows from the blade pressure side to the suction side through the gap. Because the plasma-free hemoglobin concentration and the local MIH value are linked through the Equation ( 1) and the relation A = ⋅ 10 , high MIH values are observable in that section. Furthermore, Figure 10 shows that there is a region with a high rate of hemolysis production at the leading edge of the blade. Thus, this region appears to be a source of lysis, too. However, only a small area with high MIH values can be identified on the blade suction side in channel 1. That is because the blood with a high plasma-free hemoglobin concentration at the leading edge is transported into the direction of the impeller casing by centrifugal forces. Together with the damaged blood from the leakage flow it is accumulated in the leakage vortex. This vortex structure can be clearly seen in Figure 9

CONCLUSION

A numerical analysis and optimization of a transplantable heart pump was presented. It was possible to reduce the pump dimension by 60 % without negative changes in the pump characteristics. There are even improvements in the new design, for example a higher hydraulic efficieny and a reduced potential for blood damage. Furthermore, regions of pressure losses were analyzed and visualized by using CFD methods to find the main reasons for the better characteristics of the new, miniaturized turbomachine.

In order to further improve future pump designs, two new concepts were tested in the paper.

The first one was the analysis of local blood damage by using A &&&&&&&& . Sources of hemolysis were found at the blade's tip and leading edge and the impeller casing region. The idea of finding local blood damage sources could be used to further optimize impeller parameter like the blade's leading edge angle or the clearance gap.

The second concept was the optimization of the airfoil shape regarding hemolysis at pulsating inflow conditions. Unfortunately, the optimization process is still in progress. These results will be shown in future studies.

NOMENCLATURE

  (a) Bridge to bridge: use until another therapy can be started (often used for a couple of hours only). (b) Bridge to recovery: use until a heart recovery happened (it happens in few cases only). (c) Bridge to transplant: use until a heart transplant is available (often used for a couple of months or years). (d) Destination therapy: use as a final solution until the end of life.

Figure 1 .

 1 Figure 1. Placement of the existing VAD-design.

Figure 2 .

 2 Figure 2. Volume flow rate equivalent to the meridional velocity for the simulated cardiac cycle.

Figure 3 .

 3 Figure 3. Bi-super-ellipse-profile definitions for a) camber line b) thickness distribution and c) the superposition of camber and thickness [9].

Figure 4 .

 4 Figure 4. Placement in the body of the new VAD compared to the existing one.

Figure 5 .

 5 Figure 5. Hydraulic efficiency of both pump designs.

  to the guide vane. As a result, the integral values for the pressure loss inside the guide vanes (7157 Pa) are almost twice as high as inside the volute (3637 Pa).

Figure 6 .

 6 Figure 6. Normalized pressure head curve of both impellers.

Figure 8 .

 8 Figure 8. Normalized global MIH of both impellers.

Figure 7 .

 7 Figure 7. Comparison of normalized pressure loss coefficient within the outlet guide vane (left side) and the volute (right side). N ^ definition bases on equation (11).

  The time averaged, local MIH valuesA &&&&&&&& are displayed in a cut plane through the impeller domain in Figure9. The values are normalized with the global MIH value calculated with Eqation (6). High magnitudes above be seen in the clearance gap region nearby the impeller casing. A reason for the high MIH values can be found in high shear stresses, which are present in the outer boundary layer owing to the no-slip condition of the fluid.

Figure 9 .

 9 Figure 9. Contour plot of the normalized, time averaged local MIH values in a cut plane through the rotational domain.

  as a large area with high local MIH values above A &&&&&&&& ! &&&&&&&&&& i ≥ 4 nearby the impeller casing and the inlet of channel 1. The plasma-free hemoglobin concentration is distributed within the blade channel while the convective movement of the fluid. As a consequence, large areas with ………. noticeable at the lower cut plane in channel 1. A Figure 10. Contour plot of the time averaged source term %̅ on the impeller surface in an isometric view, the cut plane shows the location of the contour plot in Figure 9.

Abbreviations

BEP:

Best Efficiency Point

MTBF:

Mean Time Between Failures

SST: Shear Stress Transport

VAD:

Ventricular Assistent Device

MIH:

Modified Index of Hemolysis


Shuai Yang

Dazhuan Wu
email: wudazhuan@zju.edu.cn

Peng Wu

Leqin Wang

Investigation on the transient characteristics of the pump system using MOC-CFD coupled method

Keywords:

The dynamic characteristics of pump response to transient events were investigated by combining the Method of Characteristic (MOC) and Computational Fluid Dynamics (CFD) together. In a typical pump-pipeline-valve system, similar to the reactor system, the pump is treated as three-dimensional CFD model using Fluent code, whereas the rest is represented by one-dimensional components using MOC code. Meanwhile, a description of the coupling theory and procedure ensuring proper communication within the two codes is given. Then the pump transient characteristics caused by downstream valve closing and opening were analyzed and compared with the steady characteristic to study the effect of fluid inertia on the pump transient performance. It was found that the H-Q curve in transient operation evidently deviates from the steady-state value and shows two distinct patterns of deviation, and the cause of the deviation was further explained by the comparison of pump internal and external characteristics. All the results

showed that MOC-CFD is an efficient and promising way for simulating the interaction between pump model and piping system.

INTRODUCTION

The security and stability of the pump and pipe system have received extensive attention in hydraulic engineering, especially in nuclear system [START_REF] Gao | Analysis of reactor coolant pump transient performance in primary coolant system during start-up period[END_REF][START_REF] Zhang | Simulation of hydraulic transients in hydropower systems using the 1-D-3-D coupling approach[END_REF]. Transient operations, such as the rapid closure of a valve, a sudden switch over and start procedure of pumps and pump failure will cause sudden flow variation in the pipe system. This phenomenon, known as water hammer, will produces severe impact force and cause transient effects on the pump during the hydraulic transition process [START_REF] Li | Numerical simulation of the transient flow in a centrifugal pump during starting period[END_REF][START_REF] Tanaka | Transient behavior of a cavitating centrifugal pump at rapid change in operating conditions-Part 1: Transient phenomena at opening/closure of discharge valve[END_REF]. The dynamic interaction between transient pipe flow and pump are closely related with the safety and the operating stability of the whole system. So it is essential to study the transient interaction between them.

Although dynamic interaction between water hammer and pump had been validated by Ismaier and Schlücker [START_REF] Ismaier | Fluid dynamic interaction between water hammer and centrifugal pumps[END_REF] through an experiment, numerical calculation was not conducted due to the lack of effective solution methods. Up to now, one dimensional MOC (Method of Characteristic) and three dimensional CFD (Computational Fluid Dynamics) analyses have been two relatively independent ''modeling cultures''. MOC is commonly used in solving the hydraulic transients in pipe and giving important information at system levels because of its feasibility and advantages in modeling complex systems [START_REF] Afshar | Water hammer simulation by implicit method of characteristic[END_REF][START_REF] Deng | Hydraulic transients induced by pigging operation in pipeline with a long slope[END_REF][START_REF] Tian | Numerical simulation and optimization on valve-induced water hammer characteristics for parallel pump feedwater system[END_REF]. CFD is increasingly being used to capture complex local 3D features inside the fluid machinery [START_REF] Ding | Demonstration and validation of a 3D CFD simulation tool predicting pump performance and cavitation for industrial applications[END_REF]. In general it is impractical to simulate a whole system with a pure 3D CFD model considering computational resource and time, and consequently separate analysis of individual components within the system, such as a pump or a pipe, is considered in isolation. However, the 1D MOC is not mature enough to be applied in the calculation of the hydrodynamic characteristics of a pump under transient operating conditions for two reasons. On one hand, it cannot provide the detailed information, like flow field structure and pressure distribution, inside the fluid equipment. On the other hand, the pump dynamic characteristics, simulated solely based on single MOC, is a quasi-steady result due to replacing the pump model with steady-state performance curve [START_REF] Al-Khomairi | Use of steady-state pump headdischarge curve for unsteady pipe flow applications[END_REF]. Therefore, the two simulations methods are to be seen as complementary, and a multi-scale modeling approach, namely MOC-CFD, can be proposed by combining both models together to utilize the strengths of both approaches [START_REF] Wu | MOC-CFD Coupled Approach for the Analysis of the Fluid Dynamic Interaction between Water Hammer and Pump[END_REF].

In this paper, MOC-CFD coupled method was used to study the pump dynamic characteristics. In a typical pump-pipeline-valve system, the downstream valve is closed and opened rapidly to achieve a rapid change in operating conditions. The transient characteristics of the pump were obtained and compared with and steady characteristics. Then the deviation of the dynamic H-Q curve from the steady-state curve was further explained through the external characteristic and the detailed flow field structure inside the pump. Meanwhile, the effectiveness and accuracy of coupled simulation method were discussed during the research.

DESCRIPTION OF HYDRAULIC SYSTEM

As demonstrated in Fig. 1, the system is an open loop, which consists of a large reservoir with a free surface at ambient pressure, a centrifugal pump, a ball valve, a turbine flowmeter and stainless steel pipes. The pump is an in-line pump with inlet and outlet diameters of 0.032 m and rotates at a constant speed of 2320 rpm during the transient events. The valve inlet and outlet diameters are the same with the pipe and the values are 0.032 m. The flowmeter is used to measure the steady-state flow velocity. The stainless steel pipes with 0.032 m internal diameter, 0.002 m wall thickness and 0.02 wall friction factor are used to connect these hydraulic parts. The length of upstream pipe is 1.2 m and the length of downstream pipe is 6 m. The vertical distance between free surface of reservoir and pipe outlet is equal to 0 m. The measured wave speed is about 1000 m/s. Fig. 1 The hydraulic system schematic diagram

NUMERICAL METHOD

In this paper, MOC-CFD coupled method was proposed and developed to analyze the hydraulic system under transient operating conditions. It includes two parts, one-dimensional calculation using the Method of Characteristic (MOC) and threedimensional simulation using the Computational Fluid Dynamics (CFD) [START_REF] Yang | Dynamic analysis of the pump system based on MOC-CFD coupled method[END_REF].

The MOC is a numerical solution of equations that govern unsteady fluid flow in pipelines. At any interior grid intersection point (Q, H) in pipeline, the two compatibility equations are given [START_REF] Wylie | Fluid Transients[END_REF]:

: P C H C BQ    (1)
:

M C H C BQ    (2
) where parameters CP and CM are functions of the flowrate and head at a previous time step. By introducing appropriate boundary conditions, the flow rate and pressure at the next time step can be easily calculated by the values in the current time step during the MOC calculation. In current study, the MOC code written using Visual Basic was developed.

CFD uses numerical methods to solve the governing equations of fluid flow. The continuity equation and momentum equations are the fundamental governing equations, which can be formulated as follows [START_REF] Ashgriz | An introduction to computational fluid dynamics[END_REF].

         V  0 t (3)        ij V F τ  D p Dt (4)
in which, ρ is the fluid density, t is time, V is the fluid velocity vector, p is pressure, τij is the stress tensor, and F is the external body force. In this paper, Fluent software was employed and the data interface is accomplished through the user-defined functions (UDF).

Based on the experimental hydraulic system, the coupled model is established, as shown in Fig. 2. The pipe and valve are modeled in MOC code and the pump is modelled in Fluent code. Meanwhile, two data interfaces were developed to realize the exchange and synchronize the numerical values between both simulation codes.

Fig.2 The schematic diagram of coupled model

In each time step, C + and C -equations determine the transient system characteristics, and Fluent solver provides pressure datum to the MOC code every certain number of iterations and then updates inlet velocity boundary to complete the coupling process once. The maximum iterations can be set to a lager value (e.g. the multiples of coupled numbers) to ensure that the co-simulation achieves convergence after several coupled data exchanges before entering the next time step. When the flow rate and head obtained from Fluent simulation also satisfy the C + and C -equation at the interface of boundaries, convergence criterion is reached. The MOC code then advances to the next time step, and updates the C + and C -equations of MOC. The Fluent also runs into the next time step to couple the MOC code. At the beginning of the next time step, the flow field inherits the results of the current time step. In addition, the time steps used in the CFD simulation are in synchronization with that used in the MOC code during the transient coupling. It is clear that the information, including pressure and volume flow, is kept the same at the shared interfaces. But in the process of data exchange, cross-section averaged quantities are transformed, either at the interface from the 1D to the 3D code or at the interface from the 3D to the 1D code. After finishing the whole simulation, transient coupling ends.

The dimensions of impeller and volute are shown and provided in Fig. 3. In order to validate the independence of the grid, mesh with different number cells in different parts had been tested. Results show that the observed variation of head was under 1% when the cell number is more than 2 million, as shown in Fig. 4. Velocity in axial direction is specified as inlet boundary, while the average static pressure field is defined as the outlet boundary. And the used model is a realizable k-ε turbulence model, which is suitable to simulate the rotating flow inside the pump. The time dependent term scheme was first order, implicit. The pressure-velocity coupling was calculated through the SIMPLE algorithm. The second order upwind scheme was applied for the spatial discretization of the momentum, turbulent kinetic energy and turbulent dissipation rate.

In addition, in the coupled simulation, the pump unit is equivalent to a pipeline node of MOC calculation, and the flows of pump inlet and outlet are approximately equal to each other, which make pump model like a data transmission carrier, so an incompressible model can be used for the simulation of pump in Fluent code. Moreover, this treatment will achieve a good numerical stability.

Fig.3 3D pump model and mesh

In the coupled model, the valve is modeled in MOC code, and its transient operations are used to generate a rapid change of flow state. For the initial steadystate flow, the valve is described as follow [START_REF] Wylie | Fluid Transients[END_REF].

 0 0 0 ( ) 2 d G Q C A gH (5)
in which Q0 and H0, respectively, are flow-rate and hydraulic loss of valve in fully open state, and (CdAG)0 the area of valve opening times the flow-rate coefficient. Here, (CdAG)0 is 0.018. For the transient flow, the relationship between the flow-rate QV and valve hydraulic loss HV are as follows [START_REF] Wylie | Fluid Transients[END_REF][START_REF] Azoury | Effect of valveclosure schedule on water hammer[END_REF].

  0 0 V V Q Q H H (6)
in which, τ is the valve opening, and the transient process of valve closure will be expressed using the equations of valve opening τ over time. Valve closing:

1.5 1 ( 0 ) 0 ( ) c c c t t t t t t                  (7)
in which tc is the valve closing time. Here, tc equals to 0.042 s. Because the duration of the unsteady flow is very short, the pump is assumed to operate with constant speed.

For the transient simulation of valve closure, analyses were started (time step t0) using initial values of the flow field obtained from the above steady cosimulation. The time step in both codes was set to 0.0012 s, which is obtained from the ratio of pipe grid distance to wave speed in pipeline and is a compromise between calculation accuracy and computational efficiency. In the transient coupled process, the Fluent solver provides pressure datum to the MOC code every 30 iterations and then updates inlet velocity boundary. At this point, one coupling is completed. The maximum iterations were set to 450 corresponding to coupling of 15 times in each time step, and this number has been proved to be enough for the co-simulation achieving convergence before entering the next time step. The valve position was then altered to go into the next time step and to update the C+ and C-equations of the MOC code. The Fluent also runs into the next time step to couple with the MOC code. In the whole process, the volume flow was again chosen as monitoring parameter, as shown in Fig. 6.

In Fig. 6, the flow-rate Q in every time step has achieved convergence before entering the next time step. In the whole process, the flow gradually decreases with the closing of the discharge valve. After the valve is fully closed, the flow fluctuation is occurred around zero and continues for a longer time under the action of water hammer.

The pump curves of inlet and outlet pressure are also plotted versus the time in seconds and compared with the experimental data, as shown in Fig. 7. In the figure, downstream valve closing will cause a positive pressure wave moving from the downstream to the pump, whose propagation leads to the fluctuations. The fluctuations also continue for a longer time until water hammer is dissipated. In addition, a conclusion can be drawn from the comparative research. The coupled numerical results had shown a good match with experimental data, which proves that the MOC-CFD coupled method is reliable

Valve opening process

In the piping system, the discharge valve is needed to adjust when the pump is completely started up. However, the downstream valve rapidly opening will change the piping system resistance, which will cause sudden flow variations in a pipe and result in water hammer [START_REF] Wu | The transient flow in a centrifugal pump during the discharge valve rapid opening process[END_REF], so this section is mainly aimed at studying the pump transient characteristics when the discharge valve rapidly opens.

The reliability of MOC-CFD coupled method has been proved through the valve-closure experiment, so this section will directly analyze the valve opening process using this method. Analysis model is the same as that used in section 3.1. However, the coupled analysis of pump opening has a certain difference compared with the above valve closing. Pump is first under the shut-off operating condition and the flow in the pipe system equals to 0 m3/h. Therefore, there is no need to carry out the steady coupling to find the initial condition and just directly set the pump inlet velocity as 0 m/s. The initial detail flow field inside the pump was obtained after 2000 iterations in the Fluent code. Then the pump operating data including the inlet and outlet pressure was transferred to MOC code to obtain the initial operating condition of the whole system.

The transient process of valve opening can also be expressed with the equations of valve opening τ over time t.

Valve opening:

                1.5 (0 ) 1 ( ) o o o t t t t t t (8)
in which to is the valve opening time. Here, to equals to 0.042 s. In this transient process, the pump is also assumed to operate with constant speed.

The transient coupled analysis for valve opening was also started from the initial steady coupled results. The parameters setting and the coupled control scheme were same with those of valve closing. In the process, the volume flow was again chosen as monitoring parameter. The history curves of pump flow-rate and pump inlet and outlet pressures were obtained and shown in Figs. 8 and Fig. 9, respectively. Fig. 8 The history curve of pump flow As shown in Fig. 8, the history curve of flow-rate Q starts from 0 m3/s and gradually increases with the opening of the discharge valve. However, the curve represents a slight fluctuation. After the valve is fully opened, the flow reached the maximum value 4.78 m 3 /h. Fig. 9 The history curve of pump inlet and outlet pressure

During the whole process, the inlet and outlet pressure of the pump synchronizing with the flow curves represent a considerable difference, as shown in Fig. 9. The fluctuation of pressure curves is much stronger than that of flow-rate curve. The valve rapidly opening will produce a low pressure wave, which moves from the downstream to the pump at the wave speed. When arrived at pump outlet position, the wave lowered the outlet pressure. The biggest pressure drop for outlet is about 50 kPa. Then the wave passed through the pump body and decreased the pump inlet pressure. Moreover, this whole process is repeated under the action of elasticity of fluid, which results in the strong fluctuations. When the valve is fully opened, the action of fluid friction and wall friction damps out the vibration and eventually causes the fluid to flow smoothly.

Comparative analysis

The transient process of valve closure and opening will cause the instability of piping system, so the pump response to the change of system resistance will present a certain transient characteristic. Therefore, the pump transient characteristics under two kinds of water hammer events were arranged and compared with the pump steady characteristic to study the effect of fluid inertia on the pump transient performance.

As shown in Fig. 10, the H-Q curves are also obtained to relate the pump flow-rate to its head. However, the H-Q curve in transient operation evidently deviates from the steady-state value and shows two distinct patterns of deviation, one is above the steady-state curve and the other is below the curve. The one above the steady-state pump H-Q curve is associated with flow reduction and is attributed to the positive pressure wave that results from downstream valve closing. The transient H-Q finally rotates around the pump shut-off point with fluctuations and the rotating radius decreases with the effect of water hammer fading away. However, downstream valve opening will generate a negative pressure wave, which moves from the valve to the pump and results in lower pressures than a pump would usually provide. The transient H-Q curve with fluctuations moves from the pump shut-off point to the mass flow points. Moreover, the wave volatility is weakened but the frequency is increased in the whole process. The wave volatility for the valve opening event is much lower than that of the valve closing event. Fig. 11 shows the relative velocity vectors distributions on two surfaces, one is defined by the inlet and outlet pipe axis intersection and the other is the middle stream surfaces S1m of the impeller. The shade of the velocity vectors indicates the vector magnitude. However, the flow field structure inside the pump under different operating state is different. For the steady case, several big primary vortexes exist near the volute tongue and rotate in the opposite rotational direction of the impeller. For the transient case of valve rapidly closing, vortexes on the middle stream surfaces S1m are squeezed, stretched, and moves to outer diameter of the impeller, and those on the surface defined by axis intersection are generated and developed when compared with those in same zone under the steady condition. However, for the transient case of valve rapidly opening, vortexes on the middle stream surfaces S1m are expanded and become much clearer, and those on the surface defined by axis intersection develop into larger vortexes and begin to block the passage when compared with the steady condition. Therefore, the positive pressure wave will strongly damage the vortex structure inside the pump, leading to the transient H-Q curve above the stead-state curve, while the negative pressure wave will promote the generation and development of vortex structure, leading to the transient H-Q curve below the stead-state curve.

CONCLUSION

In this paper, a typical pump-pipeline-valve system, similar to the reactor system, was established to study the transient characteristics of pump response to transient events. And a coupled method, namely MOC-CFD, was used to predict this dynamic interaction. The coupled method is fully flexible and capable of providing detailed information only where needed while providing system level information in the rest of the domain. What's more, this multi-scale method avoids the need of imposing approximated boundary conditions to the pump model which would badly affect the reliability of the simulation itself.

In the transient analysis for the rapid process of valve closure, the coupled numerical results had shown a good match with experimental data, which proves that the coupled method is reliable. Then the transient analysis for the rapid process of valve opening was also carried out using MOC-CFD coupled method, and the pump transient characteristics under two kinds of water hammer events were arranged and compared with the pump steady characteristic. The H-Q curve under transient operating condition evidently deviates from the steady-state value and shows two distinct patterns of deviation, one is above the steadystate curve and the other is below the curve. The one above the steady-state pump H-Q curve is associated with flow reduction and is attributed to the positive pressure wave that results from downstream valve closing. However downstream valve opening will generate a negative pressure wave, which moves from the valve to the pump and results in lower pressures than a pump would usually provide. In addition, the effect of fluid inertia caused by transient events was further analyzed from the perspective of internal flow. The results show that the positive pressure wave will strongly damage the vortex structure inside the pump, while the negative pressure wave will promote the generation and development of vortex structure, and the actions are account for deviation of the dynamic H-Q curve from the steadystate curve.

Fig. 4 3 . 1 Fig. 5

 4315 Fig.4 the pump heads with different mesh numbers

Fig. 6

 6 Fig. 6 Transient coupled process

Fig. 10

 10 Fig.10 Comparisons of the pump transient and steady characteristics

  Fig.11 comparisons of relative velocity vectors at the same flow value (2.3 m 3 /h)

CKNOWLEDGMENT

This study was carried out as a part of the National Natural Science Foundation of China (Project No. 51276213). The support is gratefully acknowledged.


Hiroki Kobayashi

Ryosuke Hagiwara

Satoshi Kawasaki

Masaharu Uchiumi

Kazuyuki Yada

Yuka Iga

Numerical Analysis of Suppression Effect of Asymmetric Slit on Cavitation Instabilities

Keywords:

In the present study, numerical analysis carried out around the cyclic flat-plate cascade with a asymmetric slit, so as to examine the suppressing or controlling effect of the slit for cavitation instabilities such as a cavitation oscillation which is similar to cavitation surge, and rotating cavitation. These instabilities cause various problem for the turbomachinery, for example, rotating cavitation causes an asynchronous axial vibration, and cavitation oscillation causes a pulsation of working fluid by resonance phenomenon of the system. Especially, in liquid propellant rocket engine, suppression process for these instabilities bring increase in cost of the launch. Therefore, it is thought that to find effective suppression technique is significant for turbomachinery.

In this paper, two types of the flat-plate three blades cascade which have an asymmetric slit on each blade were analyzed, and compare with the result of cascade without slit. As a result, the cavitation oscillation is perfectly suppressed in both of two types cascade with an asymmetric slit. Also, in the one case of them, sub-synchronous rotating cavitation and rotating stall cavitation are suppressed. Furthermore, super-synchronous rotating cavitation which is not observed in single-stage cascade, is detected by arranging an asymmetric slit. These results indicate the possibility of suppressing cavitation instabilities or controlling the type of the cavitation instabilities by the arrangement of the slit. Moreover, the head performance at the almost same cavitation number is equal or slightly increasing by arranging the asymmetric slit.

INTRODUCTION

Cavitation is a phenomenon which liquid evaporates and vapor bubbles occur in the low pressure region such as a pressure of liquid is roughly less than vapor pressure. Usually, this phenomenon is observed in the high-speed fluid machinery, for example, pumps and hydraulic machinery. Also, the cavitation exerts a bad influence such as vibration, noise, and performance decrement on these hydraulic machineries. Moreover, the cavitation brings about the oscillations which are called cavitation instabilities in the turbomachinery on the occasion. The cavitation instabilities are generally classified into two type, which is cavitation oscillation and rotating cavitation. The cavitation oscillation [START_REF] Tsujimoto | Observations of Oscillating Cavitation of an Inducer[END_REF] leads to flow rate and pressure fluctuation, also, it causes an oscillation for entire of pump systems by a resonance phenomenon. The rotating cavitation [START_REF] Hashimoto | Experimental Study on Rotating Cavitation of Rocket Propellant Pump Inducer[END_REF] causes asynchronous axial vibration in turbomachinery. As the example that the accident caused by this phenomenon, launch failure in the Japanese liquid propellant rocket happened due to the rotating cavitation in 1999. Therefore, it is important to find the suppression and control method of these instabilities for designing of the highly performance and reliable turbomachinery. For the above reasons, the methods which aim to suppress the cavitation instabilities have been developed and adopted in the liquid rocket engine such as the using pogo-suppresser [START_REF] Morino | Vision and Description) POGO Fluctuation in Liquid Rocket[END_REF] and extension of the casing diameter at inducer inlet [START_REF] Shimagaki | Effect of the Casing Configurations on the Internal Flow and unsteady Pressure Fluctuation in Rocket Pump Inducer[END_REF]. However, these methods are not perfect owing to the complex characteristics of cavitation.

In these days, we can obtain commercial software and open source codes easily, and numerical simulation by the software is able to utilize in designing of the fluid machineries. But numerical simulation has not been utilized in designing of fluid machineries with cavitation and theoretical analysis and experiment are still predominant. It is because the accuracy of the CFD for cavitation is not enough. In benchmark simulation by 6 research groups and 4 software venders by Turbomachinery society of Japan,it was reported that even time averaged lift and drag can't be predicted by the all commercial software and the in-house codes [START_REF] Kato | Industry-University Collaborative Project on "Numerical Prediction of Cavitating Flows in Hydraulic Machinery -Part1: Benchmark Test on Cavitating Hydrofoils[END_REF] . In such a situation, our group developed numerical method which is suitable for unsteady cavitation by using homogeneous model and compressible scheme, and has several experiences of the numerical simulation for cavitation instabilities, which is although arising in 2 dimensional cascade. Hence, it is considered that the numerical simulation for development of suppression technique of cavitation instabilities is possible as application of our numerical method.

In this paper, the suppression effect for cavitation instabilities by using slit on blade is investigated. A Numerical Analysis of Suppression Effect of Asymmetric Slit on Cavitation Instabilities -2 numerical simulation of unsteady cavitation in the cascade with an asymmetric slit is performed, and results are compared with cascade without slit. The occurrence region of cavitation instabilities, cavity volume, mass flow rate, and head performance of the asymmetric slit are discussed and compared with results of cascade without slit. These results indicate the suppression and control effect of an asymmetric slit for cavitation instabilities in cascade without decline of performance.

NUMERICAL METHOD 2.1 Cavitation model

In this study, a locally homogeneous model of a compressible gas-liquid two-phase medium [START_REF] Iga | Numerical Analysis of Cavitation Instabilities Arising in The Three-Blade Cascade[END_REF] was used for numerical simulation of cavitation. This model allows to consider of a gas liquid two-phase field as a pseudo-single phase medium. Therefore, the Navier-Stokes equations for continuum can be applied to a cavitating flow field in which there is discontinuity between the gas and liquid phases. As a governing equation, the compressible gas-liquid Navier-Stokes equations are used for the gas-liquid medium.

 

j j t x         EE Q S (1)
,

i u Y          Q j j i j ij j u uu p uY                 E 0 , 0 j ij         E 0 0        S ,
where ρ, p, and u are the density, static pressure, and velocity of the mixture phase respectably. The mass fraction of the gas phase is denoted Y. u and μ , which is necessary for estimating the stress tensor τ in Eq. ( 1), is determined from the following equation:

                  i j j i ij ij x u x u u div 2 1 2     , (2)
   1 1 2.5 l g          , ( 3
)
where 𝛼 is the volume fraction of the gas phase which means the void fraction, and subscript g and l mean the gas and liquid phase. The government equations in Eq. ( 1) are closed by an equation of the state for the locally homogeneous compressible gas-liquid two-phase medium [START_REF] Iga | Numerical Analysis of Cavitation Instabilities Arising in The Three-Blade Cascade[END_REF]. In this equation, the local equilibrium between gas and liquid phase at a pressure and temperature, and a linear combination of the masses of the gas and liquid phases are assumed. Additionally, the liquid is assumed to be compressible, and the gas is assumed to obey the idealgas low:

        1 c l c g c p p p K Y p T T R Y p p T        , (4)
where 𝐾 𝑙 and 𝑅 𝑔 are the liquid constant and gas constant. 𝑝 𝑐 and 𝑇 𝑐 are the pressure and temperature constant of the liquid, respectively, and are assumed to be constant. The speed of sound in the two-phase medium is determined from Eq. ( 4) is found to be in good agreement with that obtained in an experiment on variations in the void fraction [START_REF] Iga | Numerical Analysis of Cavitation Instabilities Arising in The Three-Blade Cascade[END_REF]. Thus, this numerical method is able to reproduce pressure wave propagation in a gas-liquid flow field. This is important for numerical simulation of the mutual interface between cavitation and fluid machinery systems.

Γ in the source term of S denotes instantaneous equilibrium evaporation model [START_REF] Iga | Numerical Analysis of Cavitation Instabilities Arising in The Three-Blade Cascade[END_REF]. In the evaporation model, an empirical constant is not included. Therefore, it is robust to rapid pressure jumps during pressure wave propagation. For the reason above, this model can be applied to flow field such as the cavitation oscillation and other phenomenon.

In this locally homogeneous model, the cavity surface which is expressed as the void gradient, can be treated as the contact discontinuity problem in a compressible fluid. Therefore, the thickness of the cavity surface depends on the resolution of the computational grids. However, it can reproduce various types of unsteady cavitation.

Numerical scheme

In this study, the discretization of the governing equations Eq. ( 1) is based on finite difference method. It is necessary to simulate stably the cavity flows, which has a large density jump between gas and liquid phase. Accordingly, the total variation diminishing (TVD) scheme is applied to preserve the monotonicity of the solution. Specifically, the explicit TVD-Mac Cormack scheme [START_REF] Yee | Upwind and Symmetric Shock-Capturing Schemes[END_REF] with second-order accuracy in time and space is used. In this way, no turbulence model is applied because reliable turbulence model has not been found for simulate for unsteady cavitating flow.

Validation

In order to develop a suppression technique for cavitation instabilities by using numerical simulation, numerical method which can reproduce the cavitation instabilities is needed. Our numerical method by homogeneous model and compressible scheme has been validated for time averaged and unsteady characteristics in several cavitating flow field. The present numerical method was validated [START_REF] Iga | Numerical Study of Sheet Cavitation Break-off Phenomenon on a Cascade Hydrofoil[END_REF] through comparisons with experimental data about timeaveraged lift and drag in non-cavitation and cavitation conditions around a Clark Y 11.7% single hydrofoil and the time-averaged pressure distribution on a Clark Y 11.7% cascade hydrofoil in non-cavitation conditions in several cascade arrangements. Additionally, the availability of an existing RANS turbulent model of single-phase flow for the CFD of a cavitating flow was discussed [START_REF] Iga | Mechanism of Propagation Direction of Rotating Cavitations in a Cascade[END_REF]. Although, as shown in Figure 1. which is for time averaged pressure distribution on a hydrofoil, the RANS simulation is effective in non-cavitation condition but ineffective in cavitation condition. Additionally, in break-off frequency of sheet cavitation on a hydrofoil, the frequency of developed sheet cavitation can be predicted by simulation without RANS model as shown in Figure 2., which is closer to experimental frequency than that with RANS model. Consequently, we decided not to use a turbulence model in the present study in order to reduce computational costs. However, the macroscopic disturbance in the flow field, which is caused by cavity volume variation, can be reproduced in the present numerical simulation.

Additionally, the authors have researched the cavitation instabilities arising in three-blade cyclic cascades by the present numerical method in the past. Seven kinds of cavitation instabilities, caused by several different mechanisms, have been reproduced without adding any additional models or boundary conditions [START_REF] Iga | Mechanism of Propagation Direction of Rotating Cavitations in a Cascade[END_REF], the occurrence map of cavitation instabilities is shown in Figure 3. The propagation velocity ratio of rotating cavitations and frequency of cavitation oscillation can be reproduced as shown in Figure 4. and



Experiment by Kawanami [11] : U= 5. which corresponds well-known unsteady characteristics of cavitation instabilities in inducer of liquid propellant rocket engine. Furthermore, two-kinds of break-off mechanism of cavitation in cascade was reproduced by the present numerical method, one is reentrant jet which is in the case of low cascade interference and the other is cavity surface instability acconpanied by pressure wave propagation inside the cavity which is in a case of high cascade interference [START_REF] Iga | Numerical Study of Sheet Cavitation Break-off Phenomenon on a Cascade Hydrofoil[END_REF]. Therefore, the present numerical method is available for numerical simulation of unsteady cavitation and cavitation instabilites, hence, it is considered be able to utilize for development of supression techeneque of caviation insatabilities with some extent of reliability. In the Cascade 2, the asymmetric slit is arranged in each blade. The position of the asymmetric slit in blade 1, 2 and 3 are 0.5C, 0.6C, and 0.7C, in the case of Cascade 3, the position of the asymmetric slit are 0.3C, 0.5C, and 0.55C from leading edge respectably. The width of slit is 0.05C each cascades. In present cascade arrangement, the cascade throat is located at the point of 0.48C. Consequently, the slit of blade1 is located slightly downstream of the throat in Cascade 2. On the other hand, the slit of blade 1 is located upstream of the throat in Cascade 3. As a boundary condition, constant total pressure, flow angle, and void fraction are applied at inlet. At the outlet boundary, constant static pressure is applied. Under this condition, the static pressure at the inlet is extrapolated and the inlet velocity is calculated by the inlet total pressure. Besides, the density and velocity at the outlet is extrapolated. Although the flow angle is constant in the boundary condition, the flow rate can change according to the change of prewhirl flow. Also, a nonslip condition is imposed at the wall boundary of the blade. In addition, the cyclic boundary condition is imposed on the third flow channel in the cascade in order to reproduce the circumferential instabilities in the cascade, such as rotating cavitation.

The inlet flow angle is αin fixed at a value of 9.5°and the time averaged inflow velocity Uin is approximately 13.5 m/s. The calculation is carried out with changing cavitation number σ, which is a dimensionless quantity that indicate the ease of cavitation. σ can be changed by changing the outlet static pressure. A mesh per blade-to-blade channel is 361 71 points. The distance between the inlet boundary and leading edge is two chord lengths, and the distance between the trailing edge and outlet boundary is three chord lengths. The cavitation number σ and static pressure coefficient ψ are estimated using the following equations:

Suppression of cavitation instabilities by asymmetric slit

The head performance of present cascades are drown in Figure 7. At high σ value region, the head performance of the Cascade 2 and 3 is better than that observed in the Cascade 1. Around head drop region, the head performance of Cascade 2 and 3 are almost same or slightly better than that of Cascade 1 although the head drop point is almost same. As the reason why the head increases in the cascade with slit is as following: The local angle of attack of rear blade takes the positive value by the pressure gradient between the suction side and pressure side at the slit. Therefore, the flow decelerates and the pressure rises twice at front and rear blades in the cascade, then head performance is improved in Cascade 2 and 3.

Figure 7 is also colored according to predicted cavitation instabilities. The occurrence of cavitation instabilities is judged by the visualization of aspects of the cavity and waveforms of variations in cavity volume, upstream pressure and upstream flow rate. In Figure 8 and 9, typical aspect of cavitation oscillation (C.O) and super synchronous rotating cavitation (Super-S R.C) are shown respectively [START_REF] Yuka | Numerical Analysis of Three Types of Cavitation Surge in Cascade[END_REF] In C.O, three sheet cavities develop and shed, and three cloud cavities collapse in Numerical Analysis of Suppression Effect of Asymmetric Slit on Cavitation Instabilities -7 Sub-S R.C and R-Stall C is not detected in all σ region. In Cascade 2, the Super-S R.C, which is not detected in Cascade 1, is observed. These results indicate that the cavitation instabilities are suppressed or the type is controlled at certain σ region by arranging the asymmetric slit on the blade, and changing the slit position.

Flow field around the asymmetric slit

Figure 10. shows the flow field around the slit by void fraction and velocity vector distribution (upper) and pressure distribution (lower) of the Cascade 2 (Blade 2). As shown in Figure 10, the sheet cavity length is limited by a jet flow which flows from pressure side to suction side through the slit. This slit jet is caused by pressure gradient between the slit, which is increased by collapse of the cloud cavity in the channel of pressure side and also increased by the cavity growth in the channel of the suction side. When a re-entrant jet flows into the sheet cavity from upper side of the slit jet, a vortex cavity is generated in the vicinity of trailing edge of the front blade. Also, the slit jet induces increasing of the local angle of attack of the rear blade, which promotes the growth of a small sheet cavity at leading edge of the rear blade. Then the cavity on the blade with a slit will be composed by three vortex structures as shown in upper side of Figure 3. After that, these three cavities separate from the blade, shed to downstream and collapse with interfering with each other.

Time averaged characteristics

Table . 1 shows the time average of the cavity volume in each blade and total cavity volume in cascade 2 and cascade 3 with an asymmetric slit at almost same σ condition about σ = 0.10. The cavitation instabilities are suppressed in the both cases. From Table 1, the total cavity volume values of Cascade 3 are less than the case of Cascade 2 in all blades. Also, the cavity volume of each blade is getting larger with increasing front blade length in each cascade. As mentioned in previous section, the slit has an effect to limit the cavity length by jet flow.

Therefore, in case of front blade length is shorter, the sheet cavity becomes shorter.

Table . 2 shows the time average of mass flow rate of each slit in cascade 2 and cascade 3 in same case in Table 1. In the Cascade 2, the time averaged mass flow rate increase incrementally with front blade length getting longer. However, the mass flow rate of Blade 1 in Cascade 3, is larger than the Blade 2, 3 of Cascade 3 and Blade 1 of Cascade 2. In case of front blade length is long such as Blade 2 and 3 in the Cascade 2, the effect of pressure gradient increasing is large due to the high pressure rising which is caused by cavity collapse at downstream. Also, the slit which arranged on the Blade 1 of the Cascade 3, is covered with sheet cavity long time than other slit. Then the pressure gradient is comparatively large. Because of the above reason, time averaged flow rate in each slit also varies by the random arrangement of the slit. When the intensity of the slit jet is different, each volume of the three cavities which was shown in Figure 3. Is considered to change. Then, the cavity volume is decided not only by the slit position but also by the flow rate of slit jet as shown in Table 1 in which cavity volume in Blade 1 in Cascade 2 is different from that in Blade 2 in Cascade 3 in spite of the slit positions are same in the blades. Figure 11. shows the waveform of total cavity volume (upper) and that of upstream pressure(lower) in Cascade 1 and Cascade 3 at the same cavitation number (σ = 0.2). In the Cascade 1, the total cavity volume oscillates with regular cycle and the pressure oscillates as pulsation which is the characteristic of C.O [START_REF] Yamamoto | Instabilities in an Cavitating Centrifugal Pump[END_REF]. In Cascade 3, the regular cycle is not observed in total cavity volume and upstream pressure. It means C.O is suppressed and the break-off cycle of sheet cavity in each blade is not in the same phase in Cascade 3. By the way, it is known that Strouhal number of cavity break-off cycle based on cavity length is known to remain constant [START_REF] Pham | Investigation of Unsteady Sheet Cavitation and Cloud Cavitation Mechanisms[END_REF], it means break-off cycle depends on the cavity length. Then it can be considered that the cavity break-off cycle in each blade is different each other in a cascade with asymmetric slit because the cavity volume is different as shown in Table 1. Therefore, C.O, in which cavity oscillation with constant cycle repeats in same phase in all blade, is completely suppressed in Cascade 2 and 3.

Unsteady characteristics

The waveform of cavity volume in the each blade is shown in upper figure in

Conclusion

In the present study, the numerical simulation of the cavitating flow around the cascade with an asymmetric slit is carried out, and examine the suppression and control effect of the slit for cavitation instabilities in the cascade, from the view point of head performance, occurrence region of cavitation instabilities, mass flow rate of slit and cavity volume. The result are summarized as follows.

-The all cavitation instabilities can be completely suppressed in a cascade in which one slit locates upstream of the cascade throat. On the other hand, in the case of a cascade in which all slits locate inside the cascade throat, cavitation Table 1. Time average of cavity volume in each blade and total volume in cascade 2 at σ = 0.10 and in cascade 3 at σ = 0.11 Table 2. Time average of mass flow rate of each slit in cascade 2 at σ = 0.10 and in cascade 3 at σ = 0.11 Numerical Analysis of Suppression Effect of Asymmetric Slit on Cavitation Instabilities -9 oscillation is perfectly suppressed, but some rotating cavitation is observed.

-By arranging the slit as random in each blade, the cavity volume and the slit jet fluctuates nonperiodically in each cycle in each blade. It is because that pressure gradient at the slit is different in each cycle in each blade because of the relative situation between pressure side and suction side at the slit is different in each cycle in the asymmetric slit cascade -The head performance of the cascade is not declined and slightly improved in some case by arranging the asymmetric slit in each blade.

Figure 1 .Figure 2 .

 12 Figure 1. Grid convergence characteristics of the single wing surface mean pressure distribution, Comparison between RANS turbulence model and laminar flow analysis (NACA0015, αin = 8deg)

Figure 3 .Figure. 4

 34 Figure 3. Occurrence map of cavitation instabilities of the present three-blade cyclic cascade [9]

Figure 5 . 3 . 3 Figure 6 .

 5336 Figure 5. Comparison with empirical frequency and actualfrequencies of cavitation oscillation[START_REF] Yuka | Influence of Pipe Length on Cavitation Surge Frequency in a Cascade[END_REF]

  pv is the saturated vapor pressure. σ and ψ are calculated from the time-average flow field of the numerical result. The flow rate coefficient the flow angle αin. In this study,  is corresponding to 0.0963, where Ua and Ut are the axial inlet flow velocity and the circumference inlet flow velocity respectably.

Figure 7 . 6 Figure 9 .Figure 8 .

 7698 Figure 7. The head performance and occurrence region of cavitation instabilities

Figure 11 .

 11 in cascade 3 in σ = 0.2. The wave form indicates that the non-periodic waveform of total cavity volume in Cascade 3 shown in Figure 11 (upper) is yielded not by a difference of a cycle in each blade but by non-periodicity of each cycle in each blade. Also the waveform of mass flow rate of each slit jet has irregularity as shown in lower figure in Figure5. The result indicates that pressure gradient at the slit is different in each cycle in each blade because

Figure 10 .

 10 Figure 10. The void and vector distribution (upper), and pressure distribution (lower) of Blade 2 in Cascade 2 (σ =0.13)

Figure 12 .

 12 Therefore, strong non-periodicity of cavity oscillation can be realized and C.O is completely suppressed by arranging asymmetric slit on cascade. Besides, as shown in Figure7, cavitation instabilities are completely suppressed in Cascade 3, on the other hand, rotating cavitations are observed in Cascade 2 although C.O is completely suppressed. Judging from the result, the non-periodicity of cavity oscillation in Cascade 3 is considered to be stronger than that in Cascade 2, which relate to the slit position.

Figure 11 .Figure 12 .

 1112 Figure 11. Comparison between Cascade 1 and 3 about waveform of total cavity volume (upper) and upstream pressure (lower) (σ = 0.2)