Jaeyeol Cho

Nahmkeon Hur
email: nhur@sogang.ac.kr

Jongrak Choi

Jiwon Yoon

Numerical simulation of oil and air two-phase flow in a planetary gear system using the overset mesh technique

Keywords:

Planetary gear systems are widely used in differentials, transmissions and brake systems of many industrial machineries. Often, the components inside such a system are fully or partially submerged in oil, which serves for simultaneous lubrication and cooling. Hence, this oil recirculation is a key element for the durability and performance. The focus of this paper is two-phase flow analysis of a planetary gear system at different initial oil level and, also, its application in a transaxle brake assembly. Usually, due to the complex multi-phase flow phenomena inside the transaxle chambers, their simulation is arduous, and hence, little technical information exist about the details of their multi-phase phenomena. In order to imitate a real case situation, the numerical simulations are solved in an unsteady framework using the Realizable k-ε turbulence model together with the volume of fluid (VOF) method. Overset mesh technique is used to simulate the rotating planetary gear and the brake pad, which allows more realistic simulation of gears engagement. The aims are to predict flow characteristics, quality of lubrication and also heat transfer inside the brake during its operation. The results are presented and discussed in terms of lubrication and cooling efficiency. The outcomes can be used to improve the design of existing brake oil lubrication.

.

INTRODUCTION

Planetary gears are usually used for combination of power shafts with different types of gears which might be used in differentials, transmissions and transaxle brake systems in many machineries, such as automobiles and trains. A planetary gear is composed of a sun gear at the center, planet gears and a ring gear. Often oil is used for the lubrication and cooling of the mechanical parts engaged in the rotation, including gears, brake pads and shafts. Since, the mechanical part might be partially submerged in the oil, the oil and air form two-phase phenomena, which are rather impossible to be fully predicted and understood without numerical simulations. In this case CFD method can be an appropriate solution.

Tadashi Yamada [1] simulated thermal characteristics of the the coolant inside of the powertrain. He predicted oil and air multiphase flow characteristics and conjugate heat transfer using a mesh regeneration method for rotating gears and VOF (Volume of Fluid) method. This numerical method was validated by comparing with experimental data, and the results were used for improving powertrain performance and efficiency. Yazdani and Soteriou [2] studied oil jet effect at the simple 2D gear system using a mesh regeneration technique and VOF method. They compared cooling efficiency of oil according to velocity of oil jet. Stosic [3] investigated oil and air multiphase flow and temperature characteristic inside of compressor using screw gears. Similar to the previous works, mesh regeneration method and VOF method were applied. Nirvesh Mehtha et al. [4] simulated oil flow inside a gearbox containing three axes and ten gears using Multiple Reference Frame (MRF) method for 3D rotating gears. Because of limitation in the MRF method, their simulations did not include the overlapped regions of the rotating gears. Li Li et al. [5] investigated relation between oil churning loss and oil level by gear simulation using mesh regeneration method.

Based on the difficulties and limitations of the existing numerical techniques, simulation of a compact rotating gear system in the presence of circulating oil and considering heat transfer problem in a single framework needs a very high computational resource. In addition, these kinds of simulations have to be carried out using unsteady multiphase flow assumptions. Furthermore, an appropriate mesh technique with sufficient mesh resolution is required, which usually results in a larger number of computational cells, which has been a drawback for many scholars. Moreover, the MRF method which is usually used for the simulations is not inherently appropriate for these simulations, which is due to the fact that MRF requires some simplifications in the overlapped regions, which causes considerable discrepancies between the simulations and the experiments. Hence, a few comprehensive sources and reports are available, in which all of the previous important simulation parameters are appropriately considered.

This research aims at investigation of the effects of initial oil level in planetary gear and unsteady heat transfer and oil circulation inside a transaxle brake. To do so, first, the required mesh resolution for overset mesh technique and other necessary CFD factors such as timestep are inverstigated during the planetary gear simulation. Later, the transaxle is simulated taking into account the results of the planetary simulations, and the predictions about the two-phase interaction of the oil and air, in the presence of heat transfer are investigated. The study is carried out using a commercial CFD solver, STAR-CCM+ V 10.02 on the computer cluster of CFD-ERC of Sogang University, Korea. The typical solution times for the planetary gear and transaxle are 90 hours and 180 hour using 80 cores (Intel Xeon E5-2660 v2, 2.20 GHz CPUs).

GOVERNING EQUATIONS

The numerical simulations consider the mass and momentum conservation equations for incompressible flows. In addition, for the pressure and velocity coupling, the SIMPLE scheme is adopted. This scheme consider relation between velocity and pressure correlations to enforce mass conservation and to obtain pressure distribution. The mass conservation equations can be expressed as:

() 0 j j u x     (1)
The momentum conservation equation is :

() i ji i j i ji u p uu s tx x           (2)
The VOF method is used for oil and air multiphase flow simulation. This model considers the volume fraction for each fluid and tracks free surface of the immiscible fluid regions. The location of interface between oil and air can be determined from the distribution of volume fraction in domain. Fluid properties such as density and viscosity are calculated based on volume fractional average. The equations physical properties are calculated as :

ii i     ( 3 -a ) ii i     ( 3 -b ) () pi i p ii c c      ( 3 -c )
where i  is the volume fraction of phase I; and evidently:

1 i    (4)
The conservation equation for the volume fraction of one phase is written by:

() () i Vi Si g ii i dV v v da t D sd V Dt            (5)
Because of fast rotating gears, the flow regime is turbulent, (Re=10,000,000). To take into account the effects of turbulent flows inside, the Realizable k- turbulence model is used. This model is based on transport equation for turbulent kinetic energy and dissipation rate and has been reported to have certain advantages in the modeling rotating and swirling flows. The transport equations of kinetic energy and dissipation rate can be presented as:

, () [ ( ) ] 1 ( ) 2 ( ) 3 i j jk j ii ij jh t i ii t ii k ku k tx x ug S xx uu k xx                        (6) 1 2 32 4 , () [ ( ) ] 2 [ ( ) ] 3 1 ( ) t j jj ii tt i ii t ht i i u tx x uu CP k kx g u CC C kx kx                                         (7)
where k1 234 0.009, 1.0, 1.22, 0.9, 0.9, 1.44,

C 1.92, C 1.44, C 0.419 hm CC              
Overset mesh, also known as Chimera, is used to discretize a computational domain with multiple different meshes that overlap each other in an arbitrary manner. Overset mesh does not need any mesh modification after generating the initial mesh. An overset mesh simulation needs background cells enclosing the entire solution domain and multiple smaller regions containing the bodies with in the domain. In overset mesh, cells are consist of active, inactive and acceptor types. At the active cells, discretized governing equations are solved. At the inactive cells, no equation is solved. The active cells and inactive cells change interactively by motion of bodies. In addition, acceptor cells separate active and inactive cells between background region and overset region and these cells are used to couple solutions between two overlapping grids. Solution values at donor cell at background region and overset region transferred to acceptor cell through interpolation. The donor cells are the active cells from other mesh that are nearest the acceptor cell. All the solution between donor cell and acceptor cell are calculated simultaneously. This means all the meshes are implicitly coupled.

GEOMETRY AND MESH RESOLUTION

Planetary gears are one sub-assembly part in the transaxles. Hence, prior to the simulations of the transaxle, a series of simulations for a panletary gear has been conducted in order to provide a technical CFD foundation for the transaxle. Planetary gear consists of a sun gear, planet gears and a ring gear. The dimensions of the simulated planetary gear are sun gear diameter D2=105 mm, planet gear diameter To make proper overset mesh, the smallest hexahedral cell size is chosen as 1mm. Therefore, two cells are located in narrowest gear gap. To make acceptor cells stable, cells have to be overlapped as many as possible. For sufficient number of mesh for gear gap, multiple prism layers are used at the wall surfaces of gears. The numerical simulations have been tested for different number of prism layers, and it has been found that at least 4 prism layers within the 2 mm thickness are required on each gear wall to avoid the numerical divergence. The CFD mesh domain contains four rotating and on stationary sub-domains. Finally, considering all necessary mesh parameters, the total number of hexahedral cells become approximately four million for the plenetary gear simulations. A cutsection of the mesh resolution is shown in Fig. 3. In addition, the transaxle brake consists of brake pad, planetary gear system and oil reservoir. Brake pad is engaged with the sun gear of planetary gear system. Brake pad has four holes which are designed for oil exchange between reservoir and planetary gear room. The same geometry and mesh resolution with the planetary gear simulations has been chosen resulting in a total eight million hexahedral cells for this model

RESULTS AND DISCUSSIONS

For the CFD simulations, the oil density, specific heat, conductivity and viscosity are assumed as 0.8461 Kg/m3, 0.4683 kcal/kg•K and 0.1141 kcal/m•hr•K and 46.35 cSt, respectively. To study the effects of ilnitila oil level on the heat transfer coefficient, the planetary assembly with different initial oil levels 30%, 50% and 70% is investigated firstly. After this part, the transaxle is simulated with the oil level of 50%. The rotation rate of planetary gear system in transaxle brake is same with planetary gear simulation. This operating condition for the transaxle is defined to be very close to the highest speed condition of an industrial transaxle brake.

Planetary gear simulation

The planetary gear simulations are continued until the steady state condition is achieved (the steady state condition judged by instanteneous volume fraction of oil distribution, the heat transfer coefficient of planetary gear surface and churing energy loss of oil in whole domain and is observed after approximate operational time of 2 sec.). The values of steady state volume of fraction distribution, velocity magnitude and also heat coefficient versus the initial oil level are showen in Figures 678, respectively. As it can be seen, for the cases with higher oil level, more oil is dispersed on the gears. In the Figure 6, volume fraction of oil distribution according to initial oil level are shown. In case of 30%, oil does not considerably reach to the central and upper regions of planetary gear, and merely trivially scattered oil reachs to middle and upper regions. In case of 50%, oil reaches higher, but still only relatively low amount of oil reachs to the upper region. However, when the initial level of oil is 70%, all regions are noticeably effected by lump of oil. In the Figure 7, velocity magnitude distribution after the operational time of 2 second are shown. Near the sun gear teeth, velocity magnitude is highest because of the effects of sun gear rotation. In the Figure 8, average heat transfer coefficient of gear surface is shown according to initial oil level. In case of 30%,50% and 70 % the heat transfer coefficients are claculated as 7000 W/m 2 •K, 11000 W/m 2 •K and 17000 W/m 2 •K, respectively. These results show that initial oil level is a crucial factor for improving the cooling effects in a planetary gear assembly. Hence, to enhance the cooling effect, increasing the amount of oil can be considered. In the Figure 9, churning energy loss at steady state is shown according to initial oil level. In case of 30% churning energy loss is calculated about 0.0002 J/kg, about 0.0005 J/kg for 50% and about 0.0012 J/kg for 70% initial oil level. These results show that the higher oil level cause more energy loss of rotating gear. When initial oil level increases 30% to 70%, energy loss increases more then six times higher. This energy loss is related to the oil resistance against the rotation of the gears. A very similar trend has been reported by other scholars [1,5], in which the increasing the oil level causes increasing churning energy loss.

. Transaxle brake simulation

For the transaxle brake simulations, oil properties are same with the planetary gear system, the and initial oil level is assumed to be 50%. Rotation speed of the planetary gear of brake is same with previous planetary gear simulation. Brake pad is engaged with planet gear cover and it rotates at the same rotating speed. When brake is operating, the heat is generated on brake pad and gears surfaces. The heat generation during the operation is calculated and assumed to be 1500 W for brake pad and 600W for gear surface (these values are assumed to be close to the real case situation). As shown in Figure 3, there exists a cold oil reservoir at the left side and planetary gears at the right side of the brake pad. In addition, four holes are designed on the brake pad for oil transfer between the left and right sides for the sake of cooling. Also these holes are made to increase surface attached with oil and air for effective cooling whole brake pad. However, the problem of overheating still exists, which is caused by ineffective oil flow between the cold oil reservoir and the planetary chamber. This issue is related to the churning of the oil around the brake pad and oil despersion apart form the brake pad as a result of centrifugal forces. Similar to the previous part, the simulation is continued until the results of VOF distribution become steady state (after approximately 2 sec.). Hence, as it can be seen in Figure 10(d), when the flow reaches its steady state condtion, the oil flow between the two sides is severly decreased. As shown, there exist only the air phase around the brake pad. So, there is no oil recirculation between cold reservoir region and planetary gear chamber. By heat generation of brake pad, left side of the brake is heated first. At first, because of the rotation of planetary gear, heat of planetary gear dissipates well. However, as the rotation continues, the temperature inside the planetary chamber increases. This problem is also addressed by Yanzhong Wang et al.

[6] as "spin flow effect". They simulated one pair of rotating gears with stationary oil and oil jet and concluded that rotation of the gears decreases the quality of lubrication and oil circulation inside the chamber. Accordingly, in the present research the "spin flow effect" increases with the rotational speed of the brake pad. As a inevitable consequence of this "spin flow effect", heated oil remains in the the planetary chamber of the brake and the amount of heat trnasfer rate is very lower than the required amount of heat transfer for cooling planetary chamber of the brake. As time increases, this phenomenon causes the oil temperature rise inside the planetary chamber (as depicted in Figure 11).

To over come this problem, several methods are suggested.

According to planetary gear result, using more oil and filling the oil reservoir to a higher level is one solution to increase cooling effect. It is also helpful to devise a by-pass pipeline for oil exchange between the oil reservoir and planetary chamber, in order to compensate the spin flow effect.

CONCLUDING REMARKS

In the present research, the unsteady turbulent multiphase flow analysis of a planetary gear system and a transaxle brake are simulated in the presence of heat transfer using the overset mesh technique. Prior to run the transaxle simulation, the planetary gear is simulated using similar fluid properties at different initila levels of oil in order to achieve the appropriate resolution for the overset mesh technique and other CFD settings. The planetary gear simulations are carried out for oil levels of 30%,50% and 70%. The heat transfer coefficient versus oil level has been calculated. In addition, the oil churning energy loss versus oil level is plotted. Using the same mesh resolutin, the transaxle brake simulation has been conducted with 50% oil level, and the phase distribution and heat characteristics inside the brake are obtained. Furthermore, the reason of overheating problem of transaxles is found to be the blockage of oil circulation inside the chamber between two sides of the brake pad (spin flow effect). This phenomenon is mainly caused by the brake pad fast rotation. According to the results, it is suggested to add a by-pass oil pipe to connect two oil chambers on both sides of the brake pad in order to compensate the spin flow effect.

Figure 1 .

 1 Figure 1. Geometry of planetary gear

Figure 2 .Figure 3 .

 23 Figure 2. Mesh arrangement of planetary gear

Figure 4 .

 4 Figure 4. Mesh arrangement of transaxle brake

Figure 5 .

 5 Figure 5. Brake heat generation region

Figure 6 .

 6 Figure 6. Steady state oil distribution in planetary gear according to initial oil level (30%, 50%, 70%)

Figure 7 .Figure 8 .Figure 9 .

 789 Figure 7. Steady state velocity distribution in planetary gear according to oil level (30%, 50%, 70%)

Figure 10 .Figure 11 .

 1011 Figure 10. Volume fraction of oil inside of brake

Table 1 .

 1 Rotation rate of gears

	Sun gear	758 rpm
	Planet gear	Rotation Revolution	363 rpm 136 rpm
	Planet gear cover	136 rpm
	Ring gear	stationary

Mass

N. Stosic, On heat transfer in screw compressors, International Journal of Heat and Mass Transfer, 72:517-530, 2015.


Junghyun Lee

Nahmkeon Hur
email: nhur@sogang.ac.kr

Insik Yoon

Numerical investigation of cavitation in centrifugal pump impeller for prediction of NPSH

Keywords:

Reducing cavitation problem in large hydraulic turbo machineries such as pumps has been a hot topic for many years. This problem is understood to be directly related to the inlet pressure and flow rate and also a parameter known as Net Positive Suction Head (NPSH). Cavitation negatively affects the pump performance and even can damage the impellers and blades. In this paper, the cavitating flow in a centrifugal pump is numerically studied in order to calculate the NPSH value. Because of low NPSH, the fluid on the pressure side of the impeller evaporates. Therefore, cavitation occurs over the pressure side of impeller, starting from leading edge towards the trailing edge. For the CFD simulations, the simplified Rayleigh equation and the VOF (Volume of Fluid) method are adopted to simulate the cavitation phenomenon. The two-phase flow simulations have been carried out under the various operational conditions. The results are presented and compared in terms of the effects of hydraulic head and NPSH on cavitation.

INTRODUCTION

Cavitation is one of the hazardous fluid dynamic phenomenon which is frequently encountered in hydraulic turbo machineries. It is a complex process that is influenced by many factors such as fluid static pressure and saturation temperature, and it may cause serious damage to blades and casing. Hence, many experiments and numerical simulations have been carried out on this problem from various aspects. Among them, computational fluid dynamics (CFD) provides a wide range of possibilities to predict the flow parameters and to observe the phenomena that may occur under different conditions.

Medvitz et al. [START_REF] Medvitz | Performance analysis of cavitating flow in centrifugal pumps using multiphase CFD[END_REF] captured the characteristic performance trends associated with off-design flow as well as the rapid drop in head coefficient at low cavitation numbers. Coutier-Delgosha et al. [START_REF] Coutier-Delgosha | Stoffel Experimental and numerical studies in a centrifugal pump with two dimensional curved blades in cavitating condition[END_REF] investigated two-dimensional curvature blade geometry in cavitating and non cavitating conditions using different experimental techniques. They also simulated a three-dimensional numerical model in order to study cavitating flows. Steinmann et al. [START_REF] Steinmann | Numerical and experimental investigations of the unsteady cavitating flow in a vortex pump[END_REF] investigated the unsteady cavitating flow in a vortex pump. Their CFD results turned out to be appropriate to increase the knowledge about vortex pumps but the calculated values were not accurate. Zhao et al. [START_REF] Zhao | Numerical simulation of cavitation flow under high pressure and temperature[END_REF] simulated a cavitation flow on a 2D NACA0015 hydrofoil under high pressure and temperature. They concluded that the heat absorption in the cavitation area reduces the local temperature and the saturated vapor pressure. Grist [START_REF] Pouffary | Numerical simulation of 3D cavitating flows: analysis of cavitation head drop in turbomachinery[END_REF] proposed a new concept for the Net Positive Suction Head (NPSH) performance indicator. The validity of the new indicator over a wide range of pumps operating conditions has been studied. The numerical analyses were also used to study cavitation flow field in the turbomachinery due to the impact of turbulence models [6~9]. Among various approaches to simulate the cavitating turbulent flows, the Reynolds-Averaged Navier-Stokes (RANS) equations are widely used to deal with the turbulence to obtain suitable computational resources. Lei et al. [START_REF] Lei | Numerical simulation of cavitation in a centrifugal pump at low flow rate[END_REF] predicted the critical point of cavitation in centrifugal pump. For considering the impact of the compressibility of the cavitation flow, a computational model was developed based on the full cavitation model and the RNG k-ε turbulence model. Zang and Chen [START_REF] Zhang | CHEN Numerical analysis of cavitation within slanted axial-flow pump[END_REF] expressed that the pump performance drop is relevant to the instability of cavitating flow on the blade suction surface. Liu et al. [START_REF] Liu | Hui Experimental investigation and numerical analysis of unsteady attached sheet cavitating flows in a centrifugal pump[END_REF] reported that as the cavitation number decreases, the length of the sheet cavities grows substantially. He et al. [START_REF] He | Analysis of cavitation behavior in a centrifugal pump[END_REF] analyzed the cavitating flow field in a centrifugal pump and concluded that two reasons of rapid head drop are the reduction of NPSH and growth of the energy dissipation in the vortices.

Following the above investigations, in the present study numerical investigation of cavitating flow in centrifugal pump is conducted in order to investigate the NPSH values with development of the attached sheet at the entrance of the impeller. The reduction of total dynamic head in low NPSH condition is discussed. The numerical simulations are solved by computer cluster at Multi-phenomena CFD Engineering Research Center at Sogang University. The cluster has AMD Opteron 64 bit processor (2.66 GHz). Typical computational time for each simulation of the cavitation is 10 days using 36 parallel cores.

GOVERNING EQUATIONS

In order to poccisely predict the flow inside the impeller, unsteady turbulent simulations are adopted. The simulations have been conducted using both multiple reference frame (MRF) for the non-cavitating flows, and sliding mesh for the cavitating flows. MRF is a relatively simple, robust, and efficient steady state to simulate rotating machinery. MRF assumes a weak interaction between different domains. However, sliding mesh technique can deal with strong interactions. The sliding mesh model is the most accurate method for simulating flows in multiple moving reference frames. Therefore, it is more desirable for cavitating simulations. The sliding mesh is computationally more expensive than MRF. The governing equations for the mass and momentum conservation become:

() 0 k k u tx      (1) () () i ji i j i ji uP uu S tx x          ( 2
)
where ρ is the density of the fluid, u represents the velocity vector, τ symbolizes the shear stress tensor, P denotes the pressure, and S is source term.

The non-cavitating cases directly use the above equations for single-phase liquid. However, the VOF method, which is used for cavitating simulations, assumes that all immiscible fluid phases in a control volume share the same velocity, pressure, and temperature fields. The equations are solved for an equivalent fluid. The flow properties inside the domain are linearly dependent upon the volume fraction phases. The equation employed for the VOF can be summarized as:

ii i     (3) ii i     (4) () Pi i p i i c c      (5) i i V V   (6)
where αi, ρi, μi and (cp)i are the volume fraction, density, molecular viscosity and specific heat of the i th phase, respectively. The conservation equation that describes the transport of volume fractions αi is:

() i ii ii g i VS V dD dV v v da s dV dt Dt             (7)
where i S  is the source or sink of the i th phase, v is fluid velocity, and Dρi/Dt is the material or Lagrangian derivative of the phase densities ρi.

The fluid undergoes the cavitation (phase change) when the static pressure at a particular location within the liquid falls below the saturation vapor pressure of the liquid. The cavitation option is available under VOF simulations in STAR-CCM+, and it was the simplified Rayleigh equation as:

2 2 3 sat l DRP P Dt           ( 8
)
where R is bubble radius, ρl is the liquid density, and Psat simbolizes the saturation pressure corresponding the temperature at the bubble surface, and P∞ is the pressure of the surrounding liquid in domain.

For the cavitating cases, different NPSH values are obtained by vapor pressure Pv involved in the definition of net positive suction head while the upstream pressure is kept constant. The NPSH is defined as:

in v l PP NPSH g    ( 9
)
where Pin is the total pressure at the pump inlet, Pv is the vapor pressure depends on the liquid temperature (25°C).

CFD SIMULATIONS

Domain sizes and boundary conditions

The schematic views of the centrifugal pump model including a shrouded impeller and a volute are shown in Figure 1. The pump inlet and outlet diameters are 50 mm. The impeller has five single curved blades with inlet and outlet diameters of 65 mm and 138 mm, respectively. A volute is used to adjust pressure conditions around the impeller. Impeller is considered to rotate at a speed of 3550 rpm, and volute and casing are considered as stationary domains.

For the non-cavitating cases, different flow rate values are applied to the outlet while the inlet relative pressure is assumed to be zero. In the cavitation cases, the inlet boundary condition with different pressure values is chosen (Table 1). Prior to the unsteady simulations, the models are solved at the same rpm and flow rate using steady state simulations, and the results are used as the initial values for the unsteady simulations.

Mesh sensitivity analysis

The models are meshed and solved using the commercial CFD package STAR-CCM+ Ver. 906. Polyhedral mesh is used, and five prism layers are used at the vicinity of the walls (Figure 2). A mesh sensitivity analysis based on the pressure head of single-phase is performed prior to the main simulations, in order to validate the mesh resolutions. As it can be seen from Figure 3, when the number of cells is above 2,800,000, the head values stay constant. Hence, in the present simulations, the case with 2,800,000 cells is chosen for all of the non-cavitating and cavitating simulations.

RESULTS AND DISCUSSIONS

Non-cavitating simulation (single-phase)

In order to check the ideal case and also the accuracy of the results, the single-phase steady state simulations have been performed at various flow rates from 5 m 3 /hr to 35 m 3 /hr. The rotational speed differ from 2400 rpm to 3550 rpm. These non-cavitating simulations are done by Standard k-ε turbulence model, and the MRF method is chosen. Figure 4 shows the non-cavitating H-Q curves obtained from experiments and the CFD results. The numerical simulations results show a good agreement with experimental data, especially at low rpm values.

Cavitating simulation (multi-phase)

Having single-phase simulations carried out, cavitating analyses are performed using both Standard k-ε and SST_ k-ω turbulence models. The simulations for cavitating Figure 5 shows the total dynamic head drop curves with different NPSH values obtained from the CFD results. The results are extracted after six revolutions have been passed(average over one revolution). Here, numerical simulation of the cavitation phenomenon are done by changing the inlet pressure at design flow rate (30 m 3 /hr).

When the pump pressure becomes less than a critical value(the stagnation pressure), the liquid in the centrifugal pump will change to vapor. In order to analyze the distribution and shape of the cavitation bubble in crucial conditions, Figure 6 shows the cavity inside impeller channel for different NPSH values. With the decrease in NPSH, a large area of cavity emerges in the impeller channels and expands over the pressure side, towards the blade trailing edge and also the opposite suction side of the channel. Hence, the fluid flow is affected, and the pump head drops. In general, it can be clearly seen that as the cavitation grows, the centrifugal pump performance decreases, which is beacuse of the swelling of the bubble. For NPSH=2.11 m, the impeller channels are seriously blocked by the cavities, which result in an abrupt pump head drop as shown in Figure 5.

In addition, Figure 7 presents the total pressure distributions along the blade near the inlet. With the reduction of NPSH, the pressure at the impeller inlet decreases. Pressure distribution at the impeller demonstrates that the low pressure region indicate the cavitation phenomenon.

Figure 8 shows the variation of volume fraction with different NPSH, inside the impeller at design flow rate (Q=30 m 3 /h). The cavitation first appears on the pressure side near the blade leading edge. Similar to Figure 6, with the decrease in NPSH, the cavity expands inside the impeller channels and occupies a considerable volume. Hence, the pump performance is affected (as shown in Figure 5). At this condition, the centrifugal pump cannot work normally, and the impeller blades may be damaged due to the air bubble collapse.

Finally, Figure 9 shows the velocity vector plot for different NPSH. With the decrease in NPSH, the velocity vectors become very turbulent at the blade pressure side near the leading edge and generate a sever vortex. The low velocity area expends gradually. In case of NPSH=3.77 m, the fluid flow near the pressure side of blades 2 and 3 becomes faster than other regions. This phenomenon could result in slight increase in NPSH curve, as it is reported by some researchers [START_REF] Grist | The NPSH performance indicator, NPSH (π)[END_REF][START_REF] Liu | Hui Experimental investigation and numerical analysis of unsteady attached sheet cavitating flows in a centrifugal pump[END_REF]. However, this effect has not been observed in the present research.

CONCLUDING REMARKS

In this paper, the numerical simulations of the cavitation problem have been carried out for a centrifugal pump. For the validation with experimental data, the single-phase steady state simulations are performed at various flow rates from 5 m 3 /hr to 35 m 3 /hr and rotational speed from 2400 rpm to 3550 rpm. The chacteristic curve shows a good agreement with the experimental data.

After steady state simulation, the cavitating flow are investigated at design flow rate (30 m 3 /hr) with two different turbulence models, Standard k-ε and SST_k-ω. The cavitation models are solved based on the VOF approach and simplified Rayleigh equation. Simulation results are presented and analyzed in various NPSH values. With the decrease of NPSH, the cavitation bubbles are generated starting from near the impeller leading edge. When NPSH=2.11 m, the impeller channels are seriously blocked by the cavities. Therefore, the fluid flow is affected by a large bubble sheet, and the pump head drop. Therefore, in accordance with the cavitation growth, the hydraulic head sharply decrease beacuse of blockage effect. The present paper shows that numerical simulation can provide useful information for the operating condition in centrifugal pumps. Furthermore, this study can also be used as the basis for the further research of attached sheet cavitation in centrifugal pumps

ACKNOWLEDGMENTS

Figure 1 .

 1 Figure 1. The schematic views of the pump models with impeller geometry.

Figure 2 .

 2 Figure 2. Computational grid.

Figure 3 .

 3 Figure 3. Grid sensitivity analysis.

Figure 4 .

 4 Figure 4. Comparison between simulation and experimental results at different rpm.

Figure 5 .

 5 Figure 5. The NPSH curve in different turbulence conditions.

Table 1 .

 1 Different boundary condition with NPSH values Inlet (MPa) 0.1 0.2 0.25 0.27 0.29 0.3 0.4 NPSH (m) 19.1 29.4 34.5 36.5 37.4 37.4 37.6 cases have been carried out at design flow rate (30 m 3 /hr) at a fixed rotational speed (3550 rpm). The pump cavitation performance is shown by calculating the NPSH and the total dynamic head at different inlet pressure conditions.

Figure 6 .Figure 7 .Figure 8 .Figure 7 .

 6787 Figure 6. The volume fraction contours at different NPSH values. (a: 2.84 m, b: 2.14 m, c: 2.11 m)


Shixi Ma

Dengji Zhou

Huisheng Zhang

SA-PSO Hybrid Algorithm for Gas Path Diagnostics of Gas Turbine-1 SA-PSO Hybrid Algorithm for Gas Path Diagnostics of Gas Turbine

Keywords:

With an ever-increasing power demand in the world and also growing share of power generation by gas turbines, continuous and risk-free operation of these devices is of high significance. All types of gas turbines are susceptible to performance deterioration because of the working conditions and polluting environment etc. As a result, health monitoring and performance diagnosis are two of the most important priorities of gas turbine manufacturers and users. Non-linear model based diagnostic method, with higher accuracy, is suitable for GPA(Gas Path Analysis) diagnostics. One way of non-linear model based diagnostic method is to transfer the Fault Diagnosis into an optimization problem. Genetic algorithms have been applied to solve this high-dimensional optimization problem. Although it has overcome the shortcomings of other traditional methods for easily being fallen into the local minimal, it is criticized for the slow convergence speed and sensitive to algorithm parameters. In addition, it may eventually result in diagnostic delay and wrongly diagnosis. In this paper, a new diagnostic method is used for the gas path fault diagnosis of gas turbine. The method used particle generator which was designed by means of fast SA (Simulated Annealing) to produce an optimized initial swarm of particles for PSO (Particle Swarm Optimization). Combining the global search ability of SA and the high efficiency of PSO, the new method has a good performance in diagnosis speed and accuracy..

INTRODUCTION

The national energy structure adjustment has accelerated the application of gas turbines. Now it is playing a gradually important role in many areas due to its low toxic emissions, great reliability and broad fuel adaptability, on the other hand, Gas turbine experience degradations over time that cause great concern to gas turbine users on engine reliability, availability and operating costs [START_REF] Li | Gas turbine performance prognostic for condition-based maintenance[END_REF]. In order to improve reliability and reduce maintenance costs, a further improvement named CBM (Condition Based Maintenance), based on current conditions of the gas turbine rather than historical operating data has drawn more and more attention [START_REF] Li | Gas turbine performance prognostic for condition-based maintenance[END_REF][START_REF] Fadlun | Measurement error influence on gas turbine operability for condition-based maintenance and reliabilityavailability improvement[END_REF]. This will call for the development of reliable diagnostic techniques.

There are many approaches for gas turbine fault diagnostics, such as gas path analysis, oil analysis, visual inspections, X-ray checks, vibration monitoring, noise monitoring, turbine exit spread monitoring, and so on [START_REF] Zedda | Gas turbine engine and sensor fault diagnosis using optimization techniques[END_REF][START_REF] Li | Performance-analysis-based gas turbine diagnostics: a review[END_REF]. Gas Path Analysis (GPA) is one of the most powerful tool among them [START_REF] Li | Performance-analysis-based gas turbine diagnostics: a review[END_REF], it can quickly isolate failures of engine components. and provide supports for engine maintenance.

In decades, a variety of techniques have been proposed to develop the GPA. A linear model based method was introduced to solve this problem firstly in 1967 by Urban [START_REF] Urban | Gas path analysis applied to turbine engine condition monitoring[END_REF]. And then in order to take into account the non-linearity of gas turbine behavior, a non-linear model based method with conventional optimization was firstly introduced in 1990 by Stamatis et al [START_REF] Stamatis | Gas turbine component fault identification by means of adaptive performance modeling[END_REF]. Unfortunately, conventional optimization may stop at a local minimum leading to a wrong diagnosis result. Especially when the number of components within turbine is large, the smearing effect may so strong that degraded components may not be recognized. In recent years, different methods have been used to overcome this disadvantage. Genetic Algorithm (GA) was used by Zedda and Singh in 1990 to estimate the engine performance through optimizing an objective function by a real coded GA [START_REF] Zedda | Gas turbine engine and sensor fault diagnosis using optimization techniques[END_REF][START_REF] Zedda | Gas turbine engine and sensor diagnostics. 13 th ISABE[END_REF], and the result showed a high level of accuracy. Particle Swarm Optimization (PSO) algorithm was used in a gas turbine health status estimation by Ying in 2015 [START_REF] Ying | Nonlinear steady-state model based gas turbine health status estimation approach with improved particle swarm optimization algorithm[END_REF]. The method was tested for a three-shaft marine engine and the case studies showed that the approach can accurately search and isolate the degraded components and quantify the degradation for major gas path components.

The accuracy of fault diagnosis is significantly improved by using these new intelligent optimization methods. However, the solving speed and diagnostic accuracy of these methods have an inverse trend. Faster and more reliable method is essential to diagnostic. In this paper, a Simulated Annealing (SA)-PSO hybrid algorithm is applied in gas turbine fault diagnosis. This proposed approach has been tested in six cases. The result shows that it has an excellent performance in diagnosis speed as well as the accuracy. The impact of different PSO models is investigated. Analysis and conclusions are made accordingly, compared with the typical GPA methods, such as GA, PSO, and traditional method.

G a s T u r b i n e F a u l t D i a g n o s i s U s i n g O p t i m i z a t i o n m e t h o d

When the gas turbine encountered a physical fault, it will result in the changes of characteristics. In the last step of the process, the changes of measurement parameters will be collected by sensors. The diagnosis process is to estimate the degradation extent of engine components based on sensor information, the process is shown in Fig.

1 [START_REF] Urban | Parameter selection for multiple fault diagnostics of gas turbine engines[END_REF].

During the steady state operation of the gas turbine, the relationship between measurement parameter z and the state parameter x can be represented by a nonlinear equation:

() z F x v   1
Where, ν is the measurement noise. The principle of diagnosis is shown in Fig. 2. The health parameters x of the gas turbine determines the state of the gas turbine, the change of different states are represented by the change of the measurement parameter z. Assuming an initial health parameter 𝑥̂, the prediction of the measurement parameter 𝑧̂ can be simulated by the gas turbine model. The optimization process is to minimize the objective function by adjusting the health parameters.

The objective function is the deviation between actual measurement parameters z and the predicting parameters 𝑧̂ . A minimization of the objective function is carried out iteratively until the best predicted engine health parameters 𝑥̂ for real x are obtained.

( )

i i objectionFun z z      2
The objective function is divided into different types, such as the absolute value type, square type, etc.

2 3 | ( ) 4
| i i i i objectionFun z z objectionFun z z        
Marco Zedda has compared the performance of classification in optimization process between these two types of objective functions [START_REF] Zedda | Gas turbine engine and sensor fault diagnosis using optimization techniques[END_REF], and the result showed that equation 3 was a better criterion.Thus equation 3 is adopted in this paper as the objective function.

P r o p o s e d A l g o r i t h m

There are several adjustment variables in gas turbine diagnosis, so we need to search for the objective function in the high-dimensional space. Certainly, the search process should be quick and accurate for diagnostic requirements. The searching methods can be broadly divided into 2 categories: traditional optimization methods and intelligent optimization methods.

The basic idea of traditional nonlinear programming algorithms, such as the improved newton method, gaussnewton method, gold segmentation.etc. [START_REF] Nelder | A simplex method for function minimization[END_REF][START_REF] Santa | Diagnostics of gas turbine engines based on thermo-dynamic parameters[END_REF], is to obtain the gradient direction of the objective function, and searching along the rising or falling direction to get the minimum value of the objective function. These methods have obvious advantages, such as fast calculating speed, high efficiency, reliability, and stability. Although the traditional gradient based algorithms have many advantages, the objective function in practice is generally not "flat" with many local minimum points in the feasible region. Therefore, the gradient based method may obtain the local minimum point leading to the wrong diagnostic result ultimately [START_REF] Zedda | Gas turbine engine and sensor fault diagnosis using optimization techniques[END_REF].

The intelligent optimization algorithms, that differ from the traditional methods, have the global searching ability which can avoid being plunged into the local minimal effectivly. In addition, they usually have a better performance [START_REF] Zedda | Gas turbine engine and sensor fault diagnosis using optimization techniques[END_REF][START_REF] Ying | Nonlinear steady-state model based gas turbine health status estimation approach with improved particle swarm optimization algorithm[END_REF]. These intelligent optimization algorithms includes: the Simulated Annealing (SA) algorithm, GA, PSO algorithm, etc.

Simulated Annealing Algorithm

Simulated annealing (SA) algorithm was put forward by the Metropolis [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF], based on iterative strategy for solving stochastic optimization problems. Its inspiration comes from the similarity between solids annealing process and optimization problems.

Firstly, the SA algorithm will randomly select an initial solution x and the initial temperature T in the solution space. At each step, it considers a neighbour state x' from the current state x, and probabilistically decides between moving the system to state x' or staying in state x. These probabilities ultimately lead the system to move to states that can fulfil the diagnostic rquirement or until a given computation budget has been exhausted. In fact, the probability can be determined according to Boltzmann-Gibbs:

' 1, if ( ) ( ) () exp( / ( )), otherwise f x f x Pf f KT       (5)
Where T is the temperature in each cycle, K is the Boltzmann constant, Δ = 𝑥 -𝑥 ′ . In most situations, the temperature decreases according to the following formula:

' ( ) ( ), 0 1 T t aT t a    (6)
Where a is an exponential cooling coefficient; a new state is updated as follows [START_REF] Alexandridis | A medical diagnostic tool based on radial basis function classifiers and evolutionary simulated annealing[END_REF]:

|2 1| 1 [(1 ) 1] rand x T T         (7)

Particle Swarm Optimization Algorithm

Particle Swarm Optimization (PSO) algorithm was proposed by Kennedy and Eberhart (1995) [START_REF] Kennedy | Particle swarm optimization[END_REF]. For details, a standard text can be referred Kennedyet al., 2001 and Poli et al 2007) [START_REF] Kennedy | Particle swarm optimization[END_REF][START_REF] Poli | Particle swarm optimization an overview[END_REF]. In this section a brief introduction of PSO algorithm is presented. Each solution of the optimization problem is called a particle. Particles fly in the n-dimensional search space with a certain speed, using the fitness function to evaluate the merits of the particles, and the particles dynamically adjust their position and velocity according to the flight experience and the best particle position. Standard PSO algorithm is described as follows: Supposing the searching space is d-dimensional, the N particles are grouped as the swarm. The position and velocity of the i th particle at the k th generation is represented as a d-dimensional vector:

  , , , 12 ,
, ,..

i k i k i k d ik x x x x  (8)   , , , 12 ,
, ,..

i k i k i k d ik v v v v  (9)
The i th particle's velocity and position could be update by the equation ( 8) and ( 9)

, 1 , 1 , , 2 , ( ) ( )
i k i k i k i k k i k v w v c rand Pbest x c rand Gbest x            (10) , 1 , 1 , i k i k i k x v x     ( 11
)
Where Pbest is the best current positions of all particles, Gbest is the historical best position of the whole swarm (global or neighborhood), w is the inertia coefficient, c 1 and c 2 are the acceleration factor, the values of w, c 1 , c 2 could be adjusted depending on the specific problems. In the original algorithm, w = 0.7, c 1 = c 2 = 2.

The original PSO algorithm has also some disadvantages. All individuals which learn from the same global best particle will make the algorithm converge faster. However it is inclined to access the local minimum point. Thus an improved algorithm called optimal particle disturbance strategy PSO (OPDS-PSO) is chosen [START_REF] Zhao | A perturbed particle swarm algorithm for numerical optimization[END_REF]. Compared to the original PSO, OPDS-PSO adds a small disturbance on the global best particle based on an adjustable perturbing normal distribution. The other particles will be updated by the speed updating function:

( , ) k k Gbest N Gbest    (12) , 1 , 1 , , 2 , ( ) ( )
i k i k i k i k k i k v w v c rand Pbest x c rand Gbest x             (13)
𝐺𝑏 𝑘 ′ is the global optimal particle after adding a disturbance, which is produced by the normal distribution. σ is the deviation of normal distribution function. Perturbation algorithm will let the particles move to the neighborhood of the global optimal position, so this algorithm has better searching ability to avoid premature.

Because the significant drawback of original PSO, in this paper OPDS-PSO is chosen as our test method, and σ = 0.005. As for the parameters w, c 1 , c 2 , we design w = 0.7, c 1 = c 2 = 2, as PSO Model I, Trelea [START_REF] Trelea | The particle swarm optimization algorithm[END_REF] investigated massive dynamic processes, and recommended two sets of parameters: w = 0.729, c 1 = 1.494, c 2 = 1.494 and w = 0.6, c 1 = 1.7, c 2 = 1.7. In this paper, these two sets are named as PSO models II and III. In the referenced papers, better results were obtained by using a decreasing inertia weight w from 0.9 to 0.4 linearly during the whole process [START_REF] Ying | Nonlinear steady-state model based gas turbine health status estimation approach with improved particle swarm optimization algorithm[END_REF], so it will also be tested in this paper as PSO Model IV.

Diagnostic Procedure

SA-PSO hybrid algorithm is used as an effective optimization tool to predict health parameters based on a self-adaptive gas turbine model. The main process of this diagnostic method includes two parts: particles generator based on fast SA to produce an optimized initial swarm of particles for PSO, and PSO algorithm is responsible for the precise searching in diagnosis. The main process of diagnosis is shown in Figure 3.

(1) SA initialization of particle position Firstly, a global search is implemented based on SA in feasible area, this will lock the approximate location of the global optimum. The result will be the initial point of particle swarm optimization.

(2) Objective function The objective function should reflect the difference between simulation result and the measurement parameters. (3) PSO searching process In each cycle, the speed and position of the particles is updated in PSO algorithm in accordance with the above equation for different PSO model I-IV. So the algorithm can achieve an accurate and fast seaching in the particular region to get the global optimal solution.

(4) Termination criteria Termination condition is the maximum number of iterations or minimum difference in every generation.

A p p l i c a t i o n a n d A n a l y s i s

In this paper the test object is a GE LM2500 gas turbine which is a two-shaft gas turbine engine used for driving. The basic specifications of this engine are as follows:

total air flow rate : 32 kg/s, total pressure ratio :21 power output: 29.8 MW thermal efficiency: 38%. The diagnostic process is tested by a PC with 3.3GHz i5 processor and 4GB RAM. The modelling of the gas turbine and optimization methods are carried out using Matlab. Six performance parameters to be estimated by using 11 measurement parameters in this paper is described in Table 1.

The optimization algorithm needs to search in the sixdimensional space in which the six parameter variations are from -0.2 to 0.1. -0.2 means that efficiency or flow rate is reduced by 20%. The other parameters for the method are shown in Table 2

Model Validation

Fault implanted here is a typical compressor fouling fault and turbine erosion. According to statistical analysis, more than 70% of engine performance degradation are related with compressor fouling [20]. It was assumed that compressor fouling fault causes a 3% of flow rate decrease and 1% efficiency degradation [20,21], and the turbine erosion causes 2% of flow rate increase and 1% efficiency degradation. Six cases are shown in Table 3. Case1 to case 3 are three single component faults(compressor fouling fault, high pressure turbine erosion fault and power turbine erosion fault) and case 4 to case 6 is multiple component faults(combination of these three kinds of faults). SA-PSO Model I is chosen to test the diagnosis result, the first 80 cycles are seaching by SA algorithm and the last 120 cycles is taken charge by PSO Model I. The related SA-PSO seaching results are shown in Table 4. It can be seen that the degraded components are sucessfully isolated and the predicted result is accurate due to the fact that the principle of the new method is a global optimun seaching method. It also can be seen that the multiple component faults diagnosis(case 4~case 6) is much harder comparing to the single component fault diagnosis. The error level of case 1 to case 3 is 10e-3 comparing to case 4~case 6 with a 10e-2 error level.

Result analysis

Comparison between original PSO and SA-PSO

Case 6 is chosen as the test case. It can represent the searching process due to the complexity with multiple component faults. Two diagnosis processes, with the original PSO algorithm and proposed SA-PSO method, are shown in Figure4. It can be seen that these two are both global optimum searching methods. After a number of generations the two methods both get the accurate diagnostic result but with different speeds. For the SA-PSO method, 80 cycles is made by the SA algorithm for a global seaching, and about 120 cycle is made by the PSO algorithm for refining search.

The original PSO method is much slower than the SA-PSO method. For original PSO algorithm, about 200 more cycles is needed to reach the same accuracy compared to SA-PSO method. It means that about half diagnostic time will be saved.

Table 3. Implanted degradation of different cases

Component Symbols

Implanted degradation (%)

Case1 Case2 Case3 Case4 Case5 Case6 CP DGC -3 0 0 -3 0 -3 DEC -1 0 0 -1 0 -1 HT DGT 0 2 0 2 2 0 DET 0 -1 0 -1 -1 0 PT DGP 0 0 2 0 2 2 DEP 0 0 -1 0 -1 -1

Comparison between different PSO models

Figure 5 shows the diagnosis results with different PSO algorithm (model I-IV) of the same case (Case 6). Here, the same particle swarm simulated by annealing generator isused. It can be seen that different types of particle swarm algorithm can achieve the same accurate diagnosis results. However, the convergence processes have slightly differences. Model II has a fastest convergence speed, and model IV has the worst performance. For above models, most particles can get the area near the global optimization point after the annealing process. This will make the particle inertia be the most important factor to affcet the convergence process. Mode II is the fastest model for using the largest particle inertia and Model IV have a good speed at the begining but then it become much slower for the inertia value become smaller.

Comparison among different optimization methods

Traditional optimization method and GA method are commonly used in model based gas turbine falut diagnosis. GA algorithm is firstly used in the gas turbine fault diagnosis by Zedda and Singh for its global seaching ability.

In order to get accurate results, several problems were found such as slow diagnosis speed and difficult selection with proper parameters. In this paper, a comparison between the SA-PSO, GA and traditional method is made to test the effectivencess of the proposed method for Case 6.

Considing that both GA and SA-PSO method have global searching abilities, so here the covergence speed is taken as the key goal to be compared . Figure 6 shows the coveragence process with these two methods. It can be seen that the GA method is much slower than the SA-PSO method. The time cost by the SA-PSO method is on the range of 2~3 seconds involving generation number of 200 by using a PC with 3.3GHz i5 processor and 4GB RAM. It is much shorter than the GA method which is 12~15 seconds.

Figure 7 and Figure 8 show the diagnosis process and result for the same case. It can be seen the traditional method is much faster than the SA-PSO method. The smearing effect, which usually occures with traditional method in the diagnosis process, can be effectively eliminated by SA-PSO algorithm.

Figure 6. Fitness with GA and SA-PSO algorithm

C o n c l u s i o n s

A new fault diagnosis method using Intelligent optimization algorithm is proposed in this paper and it is tested in a diagnostic platform for LM2500 gas turbine engine. The effectiveness of the method has been verified by six cases compared with the original PSO, the GA, and the traditional method. Different PSO models are also tested and it can be a reference for selecting parameters. Some conclusions have been obtained as follows:

(1) The method can be used successfully for single and multi component degradation diagnosis. [START_REF] Fadlun | Measurement error influence on gas turbine operability for condition-based maintenance and reliabilityavailability improvement[END_REF] The effect of algorithm parameters on diagnostic accuracy and speed are studied. A big particle inertia is advantageous for PSO model. [START_REF] Zedda | Gas turbine engine and sensor fault diagnosis using optimization techniques[END_REF] The diagnostic result of this new algorithm is compared with other algorithms, such as GA, the original PSO, and traditional method. The result shows that the proposed method is much faster than the GA and original PSO. And the smearing effect which is observed by traditional method can be effectively eliminated. (4) The numbers of iterations and population for diagnosis are substantially smaller than the GAbased method while guarantee the accuracy. Diagnostic time, on a PC with i5 processor of 3.3GHz and 4GB RAM, is on the range of 2-3s, which is much shorter than GA based method (12-15s). The result shows that the method is suitable to be used in gas path diagnostics.

Figure 2 .

 2 Figure 2. Model based gas turbine fault diagnostics

Figure 1 .

 1 Figure 1. Gas tur ine fault diagnosti s approa h

Figure 3 .

 3 Figure 3. Fault diagnostic procedure of SA-PSO method

Figure 4 .

 4 Figure 4. Fault diagnostic process comparison between original PSO and SA-PSO

Figure 5 .

 5 Figure 5. Fitness with different PSO Model

Figure 7 .Figure 8 .

 78 Figure 7. Fitness with traditional and SA-PSO algorithm

Table 1 .

 1 Performance parameters and measurement parameters

	Variable type	Notation	Description
		DGC	Compressor flow rate degradation
		DEC	Compressor efficiency degradation
	Performance	DGT	High pressure turbine flow rate degradation
	parameters	DET	High pressure turbine efficiency degradation
		DGP	Power turbine flow rate degradation
		DEP	Power turbine efficiency degradation
		PL	Power load
		T1	Compressor inlet temperature
		P1	Compressor inlet pressure
		n1	Gas generator speed
	Measurement parameters	P4 PRC PRP	Power-turbine exhaust pressure Compressor pressure ratio Power turbine expansion ratio
		T2	Compressor outlet temperature
		T34	high pressure turbine exhaust temperature
		T4	Power turbine exhaust temperature
		Qf	Fuel flow

Lots of work focusing on the effect of different function has been done

[4.5]

. Equation 3 is regarded as the preferred one because it provides a robust estimation.

Table 2 .

 2 : Parameters

		for SA-PSO
	SA	T_max	100
		T_min	0.1
		a	0.93
		K	0.8
	PSO	population size	5
		Max number of	200
		iterations
		minimum global	1e-6
		error
		PSO model type	I-IV

Table 4 .

 4 Diagnosis result of different cases

	Component Symbols	Case1	Case2	Diagnosis degradation(%) Case3 Case4 Case5	Case6
	CP	DGC DEC	-2.995 -0.997	7.2e-4 -3.5e-4 2.3e-3 -1.0e-3 -3.021 -0.995	6.7e-3 -1.2e-3 -1.005 -3.009
	HT	DGT DET	-2.5e-3 1.9925 6.2e-3 -0.9984 -2.2e-3 -1.008 -2.9e-4 1.984	1.995 -0.996	1.8e-2 -2.1e-2
	PT	DGP DEP	3.5e-3 -2.9e-3 -3.5e-3 -0.995 3.8e-3 2.018	-1.9e-2 2.039 3.4e-2 -0.953	2.039 -1.050

ACKNOWLEDGMENTS

This work was supported by Program of Science and Technology Committee of Shanghai Municipal.


Shuichi Torii
email: torii@mech.kumamoto-u.ac.jp.

Radzuan Razak

EFFECT OF ALUMINA NANOFLUIDS WITH DIFFERENT VOLUME FRACTION WITH DIFFERENT FILL CHARGE RATIO ON THE HEAT TRANSFER PERFORMANCE OF SELF OSCILLATING MULTI HEAT PIPE

Keywords: nanofluids, self-oscillating multi heat pipes, effective thermal conductivity

This experiment study was performed to investigate the effect of different volume fraction of alumina (Al 2 O 3 ) with different fill charge ratio to the heat transfer performance of a selfoscillating multi heat pipe at different constant heat flux. Working fluid employed was alumina (Al 2 O 3 ) with 0.1%, 0.5%, 1.0% and 5.0% volume fraction. The fill charge ratio f cr was varied from 30% to 100% at 10% interval. Heat pipe is consisted of an evaporation section, a condenser section and an adiabatic section. Evaporation and condenser sections have the same size and are connected by four circular parallel tubes. External dimension of evaporation and condenser section were a 45mm of length, a 45mm of width and a thickness of 8mm, also the internal dimension were 42mm, 42mm, and 5mm respectively. Adiabatic section (four circular parallel tubes) with a length of 45mm long, an inner diameter of 5mm and an outer diameter of 6mm was employed. The evaporation section is heated at constant heat flux for each different volume fraction and heat flux was 5W/cm 2 to 30W/cm 2 at 5W/cm 2 interval. 5 temperature measurements, 4 at the corner and center of the evaporation section and 3 temperature measurements at along the center line of condenser section were employed. 15°C of cooling water was supplied to the condenser section (to cool) in water tank at rate of 3.5l/min. A 7400Pavacuum pressure was maintained to create self-oscillating heat transfer. The thermal performance of higher volume fraction of nanofluids was enhanced in comparison with lower volume case.

INTRODUCTION

Modern technology makes it possible to produce particle less than 100 nm in diameter and for suspending in conventional fluid such as water, engine oil, and ethylene glycol. This new class of fluid is referred as "nanofluids". Choi [1] is the first person to use and introduce nanofluids as working fluid for heat transfer. Compared to micron size particles, nano particles have much larger relative surface areas and have a great heat transfer enhancement potential. Its nano size particles reduce particle clogging and reducing pumping power to other liquids that have equivalent heat transfer capability. Based on this idea, many studies are conducted to explore the superiority of nanofluids's property, such as large surface area to volume ratio, stable suspension and no flow passage clogging which are suitable in heat transfer application that millimeter size particles fluid did not offer before. This is because larger surface areas not only improve heat transfer capability but also increase the stability of suspension. Oxide nanoparticle such as Al 2 O 3 and CuO have excellent dispersion properties in water, oil and ethylene glycol and form suspension as demonstrated by Eastma et al. [2]. Torii and Yang [3] also demonstrated the superiority of nano fluid that significant enhancement of heat transfer performance due to suspension of nanodiamond particle in the circular tube flow is observed in comparison with water. Adding nanoparticles to fluid can effectively increase the thermal conductivity ratio of the fluid as studied by Teng et al. [4].

Heat pipes introduced by Akachi [5] is "The loop type of heat pipe" which is a self-oscillating type of heat pipe that contains three sections: an evaporation section, an adiabatic section and a condenser section. This type of heat pipes usually is made by a closed circuit copper tube with calculated internal diameter and is also filling with a fixed amount of working fluid. Yoon et al. [6] also investigated the heat transfer characteristics of a self-oscillating heat pipe using pure water as a working fluid. Although excellent results of thermal conductivity were obtain, this study only applied with low heat fluxes. Closed-end oscillating heat pipe (CEOHP) has more reasonable thermal performance compared to those obtained by conventional heat sink cooler as showed by Rittidech et al. [7]. Guo and Nutter [8] had a deep study and showed the result of a type of heat pipethermosyphon: An experiment study of axial conduction through a thermosyphon pipe wall. The experiment data showed that the conduction through the pipe wall causes the wall temperature to decrease along the evaporator section. It also increases the overall heat transfer coefficient, evaporation heat transfer coefficient, and condenser heat transfer coefficient of the thermosyphon. However, the heat transfer associated with axial conduction decreased as the heat flux increased. Park et al. [9] had an investigation about heat transfer characteristic of a twophase closed thermosyphon to the fill charge ratio. The author discovered dry-out phenomenon occurs when fill charge ratio is lower than 20%.

Based on the forgoing good result of researchers about application of nanofluids to heat exchanger, in this research, different volume fraction Al 2 O 3 fluid was applied to the selfoscillating multi-heat pipe to study the heat transfer performance with different fill charge ratios and different constant heat fluxes. The particle size used in this experiment is less than 30nm in diameter.

EXPERIMENT METHOD

Nanofluids applied in this experiment were Alumina (Al 2 O 3 ) nanofluid. Pure water was mixed with Al 2 O 3 particles around 30mm in diameter. The volumes fraction of nanofluid applied in this experiment were 0.1%, 0.5%, 1.0%, and 5.0%.

The following equation is one used to determine volume fraction of Alumina Al 2 O 3 nanofluid, i.e., 0.1% volume fraction.

V nnp /V nnf =0.1% Where, V nnp : the volume of nanoparticles in solution V nnf : the volume of obtain nanofluid solution After getting accurate weight of nanoparticles corresponding with volume of pure water, pure water and nanoparticles were mixed with the aid of ultrasonic vibration machine in 4 hours.

A schematic diagram of experimental apparatus is described in Fig. 1. In this figure, we can see all the equipment that is used during the experiment. It is important to point out that the self-oscillating heat pipe is made in Laboratory of Kumamoto University -Japan and the detailed drawing of the self-oscillating heat pipe can be found in Fig. 2. The heat pipe is heated on the heating section by the heater block which is made by copper containing 5 heaters inside. Also thermo-glue is used between two surfaces to improve heat transfer (5 heaters are not shown in this figure).

Here, the working fluid in the insulated fort tubes is independentlly oscilated between the cooling and heating parts, as seen in Fig. 2.

In this experiment, the working was Alumina nanofluid 0.1% volume fraction. The fill charge ratio f cr which is ratio between fill charge volume of Alumina nanofluid 0.1% volume fraction and internal evaporator section volume, was changed from 30% to 100% at 10% interval. Here, 100% implies that the evaloration section is filled with the working fluid. The heat fluxes q applied to the evaporator section for each case of fill charge ratio were 5W/cm 2 , 10W/cm 2 , 15W/cm 2 , 20W/cm 2 , 25W/cm 2 and 30W/cm 2 , respectively. After getting result and comparing heat transfer performance between the fore going fill charge ratios, we could obtain the optimal fill charge ratio of the heat pipe corresponding to the heat flux. The temperature of each section of the heat pipe was measured with thermocouples. The inlet temperature of the cooling water was set at 15°C using the controllable thermostatic bath and was supplied to the water tank at rate of 3.5l/min. After filling a calculated amount of working fluid in the heat pipe, a 7400Pa vacuum pressure was maintained to create self-oscillating heat transfer. In order to get accurate vacuum pressure value, the vacuum pressure value was read on the power indicator transducer after 20 minutes and when changing from one fill charge ratio to another, the working fluid was filled when the evaporator section of the heat pipe was 25°C. Temperature was measured until experiment reachs a steady state and the steady time is at least 30 minutes. During experiment, the adiabatic section and evaporator section was insulated with fiber glass to prevent heat loss, and the condenser section was inserted to a water tank made by plastic containing cooling water circularly flowed from the thermostatic bath.

DATA REDUCTION

Effective thermal conductivity is calculated as following equation:

  T T NA QL k L T T k NA NAq Q c H eff c H eff      ' ( 1
)
Where Q is total heat load, q' is the heat flux from evaporator section to the condenser section, L is the length from the center of evaporator section to the condenser section, approximately the distance between T H measured point and T C measured point. N is the number of tubes of the adiabatic section (N=4).

A is the heat flux area of the inner part of the tube in the adiabatic section (m 2

). A=(πd 2 )/4. From equation (1), we obtain

    T T Q d N L T T Q NA L k c H c H eff     4 / 2  (2)
Where L=0.1m, N=4, and d=0.005m

Applying with numerical value of L, N, and d we obtain,

k eff =1273.88 Q/((T H -T C )) (3)
All data collected here are calculated using the above equation and effective thermal conductivity versus fill charge ratio are shown for different volume fractions.

RESULTS AND DISCUSSION

From the result of the temperature of the evaporator and condenser section, we obtain effective thermal conductivity of the heat pipe. Figures 3 to 6 show the variation of effective thermal conductivity k eff versus fill charge ratio f cr . As shown in Figs. 3 to 6 at all six constant heat fluxes, effective thermal conductivity of heat pipe tends to increase when fill charge ratio increases.

In other words, the effective thermal conductivity becomes higher with an increaee in the fill charge ratio. From Figs. 3 to 6, the increase in effective thermal conductivity k eff is increasing gradually similar if compare to other heat fluxes. This can be explained that when heat flux is applied to the evaporator section of the heat pipe increases, high amount of heat transfer throught the wall of the heat pipe to the condenser section was cooled properly by cooling water.

Therefore, heat transfer efficiency increases. In addition, the increasing of fill charge ratio f cr helps to tranfer large amount of heat from evaporator section to condenser section to enhance heat transfer efficiency. In other words, high temperature at higher fill charge ratio makes condensation performance increase and this corresponds to an increase in heat transfer performance. Thus, this type of heat pipe shows excellent transport characteristic even for a very simple structure. effective thermal conductivity of heat pipe tends to increase when fill charge ratio increases. In other words, the higher fill charge ratio corresponds to the higher effective thermal conductivity. From Figs. 7 to 10, the effective thermal conductivity k eff is amplifying gradually with an increase in heat flux. This can be explained that when volume fraction increases the surface area increases. Thus, heat transfer coefficient increases. An increase in surface area yields the increasing of effective thermal conductivity. The volume fraction effects heat transfer performance. From this result, we can see that with the use of alumina nanofluids with different volume fraction, the effective thermal conductivity of the heat pipe enhances effectively at heat flux q=10W/cm 2 for alumina nanofluid 0.5%, 1.0% and 5.0% volume fraction and this enhancement intensifies when heat flux increases.

This can be explained that when heat flux applied to the evaporator section of the heat pipe is increased, high amount of heat transfer from the wall of the heat pipe to the condenser section was cooled properly by cooling water. Therefore, heat transfer efficiency increases.

In addition, the increasing of fill charge ratio f cr helps to transfer large amount of heat from evaporator section to condenser section. In other words, higher fill charge ratio makes condensation performance increase and this corresponds with an increasing in heat transfer performance. Note that beyond heat flux q=15W/cm 2 , in the condenser section particles tend to aggregate and this makes the suspension of concentration decrease resulting in attenuation of heat transfer enhancement. However, we can conclude that this type of heat pipe shows excellent transport characteristic even for a very simple structure.

CONCLUSION

When heat flux was applied, the higher volume fractions (%) with higher fill charge ratio f cr the higher effective thermal conductivity k eff . 1. This multi heat pipe has capability to increase heat transfer performance corresponding with the increase of heat fluxes in all four cases of using alumina (Al 2 O 3 ) with 0.1%, 0.5%, 1.0% and 5.0% volume fraction. 2. The effective heat transfer performance is intensified with an increase in volume fraction. This trend becomes larger the higher heat flux.

3. The heat transfer performance of this multi heat pipe can be enhanced around 15% when alumina (Al 2 O 3 ) with 0.1% volume fraction is replaced by alumina (Al 2 O 3 ) 5.0%. [ 5 ] A k a c h i H . , 1 9 9 0 , " S t r u c t u r e o f H e a t P i p e , " U n i t e d S t a t e p a t e n t , P a t e n t N o . 4 9 2 1 0 4 1 .

[ 6 ] Y o o n , S . H . , O h , C . , C h o i , J . H . , 2 0 0 2 , " A S t u d y O n H e a t t r a n s f e r C h a r a c t e r i s t i c o f a S e l f -O s c i l l a t i n g H e a t P i p e " K S M E I n t e r n a t i o n a l J o u r n a l , V o l . 1 6 N o . 3 , p p . 3 5 4 -3 6 2 .

Figure 1 .Figure 2 .

 12 Figure 1. Schematic design of the experimental apparatus

Figure 3 .

 3 Figure 3. Effective thermal conductivity of the heat pipe versus fill charge ratio

Figure 7 .

 7 Figure 7. Effective thermal conductivity of the heat pipe versus fill charge ratio

  e l o p m e n t s A p p l i c a t i o n o f N o n -N e w t o n i a n F l o w s , F E D -v o l . 2 3 1 / M D -v o l . 6 6 , A S M E : 9 9 -1 0 5 , e d i t e d b y D . A . S i g i n e r a n H . P . W a n g , N e w Y o r k.

[ 2 ]

 2 E a s t m a n , J . A . , S t e p h e n U . S . , C h o i , S . L . , T h o m p s o n , L . J . a n d S h i n p y o . L . , " E n h a n c e d T h e r m a l C o n d u c t i v i t y T h r o u g h t h e d e v e l o p m e n t o f N a n o f l u i d s " , P r o c e e d i n g o f t h e s y m p o s i u m o n N a n o p h a s e a n d n a n o c o m p o s i t e M a t e r i a l I I , M a t e r i a l R e s e a r c h S o c i e t y , B o n t o n , M A , V o l . 4 5 7 , p p . 3 -1 1 .[ 3 ] T o r i i , S . a n d Y a n g . W . -J., 2 b u l e n t P i p e F l o w " , J o u r n a l o f H e a t T r a n s f e r , V o l . 1 3 1 .

[ 4 ]

 4 T e n g . , T . -P., H s u , H . -G., M o , H . -E., C h e n , C . C . , 2 0 1 0 , " T h e r m a l E f f i c i e n c y o f h e a t P i p e w i t A l u m i n a N a n o f l u i d s " , J o u r n a l o f A l l o y s a n d C o m p o u n d 5

[ 7 ]

 7 R i t t i d e c h , S . , B o o n y a e m , H . and T i p n e t , T . , 2 0 0 5 , " C P U C o o l i n g o f D e s k t o p P C b y C l o s e d -e n d O s c i l l a t i n g H e a t P i W e i G . a n d D a r i n W . N . , 2 0 0 9 , " A n E x p e r i m e n t a l S t u d y o f A x i a l C o n d u c t i o n T h r o u g h a T h e r m o s y p h o n P i p e W a l l " , A p p l i e d T h e r m a l E n g i n e e r i n g , 2 9 , p p . 3536 -3 5 4 1 . [ 9 ] P a r k , Y . J . , K a n g , H . K . a n d Kim, C . J . , 2 0 0 2 , " H e a t T r a n s f e r C h a r a c t e r i s t i c o f T w o -p h a s e C l o s e d T h e r m a l s y p h o n e t o t h e F i l l C h a r g e R a t i o " , I n t e r n a t i o n a l J o u r n a l o f H e a t a n d M a s s T r a n s f e r , 4 5 , p p . 4 6 5 5 -4661. .


Xin Liu

Zhengwei Wang

Lingjiu Zhou

Xavier Escaler

Yongyao Luo

Failure investigation of a Francis turbine under the cavitation conditions

Keywords:

Acoustic-Structural coupling methodology has been used to numerically investigate the failure of a Francis turbine runner that broke down within one day of operation. The goal of the work has been to determine the feasibility of a resonance condition under cavitating flow conditions. The obtained results prove that the presence of large vapor cavities increase the natural frequencies of the diametrical modes of vibration in comparison with the pure water condition. In particular, a possible tuning between the hydrodynamic excitation and the 4ND mode has been found.

INTRODUCTION

FEM-based acoustic fluid approach has not yet been used to estimate the added mass effects on turbine runners under cavitation conditions due to the difficulties of modelling the cavity properties and the own limitations of the acoustic fluid hypothesis. On the other hand, the previous understanding in the field of the hydroelectric is that cavitation leads to structural erosion, but would change the inherent structural dynamic characteristics because of lack of test data in situ and related research. From a design point of view, the benefit in cavitation flow analysis is obvious because of the fact that many submerged bodies suffer from this phenomenon, such as pump impellers, turbine runners, injectors, marine propellers, hydrofoils and so on.

The case reported here corresponds to a Francis turbine destroyed within one day. Majority of the traditional analysis approaches such as static/dynamic stresses analysis, modal without vapor, Von Karman excitation estimation, have failed in this case. The new hope came from there was sharply noise when the turbine was running in situ and later CFD analysis found out the existence of severer cavitation. Therefore, the research purpose of this work is to establish a reliable FSI model, and to investigate whether there is some kind of relation between cavitation and runner failure by FEM-based acoustic fluid approach under cavitation conditions.

METHODS

In FSI problems, the structural dynamics equation must be considered along with the Navier-Stokes (NS) and the continuity equations. The discretized structural dynamics equation can be formulated based on the structural elements. The NS and continuity equations are simplified to get the acoustic wave equation assuming that the fluid is compressible, which means that its density changes due to pressure variations (i.e. acoustic medium), and that there is no fluid flow [START_REF]ANSYS Inc. Ansys theory reference[END_REF].

p ρR 0 M u C0 0 p +0 C 0 0 0C p K -R -R +0 K                                                          S liquid F,liquid liquid F,gas F,gas gas u F p F 0= 00 F p K                      (1)
where MF, CF and KF are the fluid mass, damping and stiffness matrices, respectively. R is a coupling matrix that represents the effective surface area associated with each node on the fluid-structure interface. FF is the applied fluid pressure vector at the interface obtained by integrating the pressure over the area of the surface. P and U are the fluid pressure and the structure displacement vectors, respectively. MS, CS and KS are the structure mass, damping and stiffness Article Title -2 matrices, respectively, and FS is the structural load vector.

Description of Francis unit and damage

This hydropower plant has three Francis units with the same design and dimensions as outlined on Figure 1. The runners have 14 blades and a diameter of 2.72 m. The wicket gates have 24 guide vanes. The runner rotational speed is 250 rpm. The rated discharge and power per unit are 164 m 3 /s and 44 MW respectively. When units 1 and 2 were put into start test during one day, almost all blades suffered a local fracture starting close to the junction between the trailing edge and the band, as shown on Figure 2.

CFD simulation under cavitation

A steady CFD numerical simulation using a Zwart-Gerber-Belamri (ZGB) cavitation model [START_REF] Zwart | A two-phase flow model for predicting cavitation dynamics[END_REF] was carried out. According to the actual trial operating conditions, one typical output power operating condition, G1, were selected corresponding to 16MW.

Using the state-of-the-art polyhedral and trimmed cell meshing technique, the grid convergence test was carried out based on the flow features. The final mesh adopted had about 227000 elements in the spiral case, about 150000 elements in the guide vanes, 347000 elements in the turbine runner and 100000 elements in the draft tube. The Reynolds-averaged Navier-Stokes equations with RNG k-ε turbulence model were applied. The ZGB cavitation model considers the total interphase mass transfer rate per unit volume, R, as indicated:

  e 3 (1 ) 2 if 3 v nuc v v vap v l P P R F P P             3 2 if 3 v v v c cond v l P P R F P P        
where the evaporation coefficient Fvap=50, the condensation coefficient Fcond=0.01, the bubble radius To visualize large cavitation structures and to determine their locations, dimensions and shapes the void ratio, α, which is defined as:

V V V  
where VV is the vapor volume and V is the total volume of the cavity. Consequently, the cavity density, c, can now be calculated as:

  1 c V L       
where V is a vapor density of and L is the water density, which have been taken as 1.7410 -2 (saturated vapor at 20 C) and as 1000 kg/m 3 , respectively.

Runner dynamic behavior in water with cavitation

The runner material parameters for the FEM model were specified according to the actual values indicated in Table 1. To consider the effect of the water, the runner was surrounded by the fluid domain. The water properties were set as 1000 kg/m 3 for density and 1450 m/s for speed of sound. The fluid mesh was built by extending the structure mesh so that the nodes on the interface could be shared by both domains, as shown in Figure 4. The numerical modal analysis has permitted to calculate the five first natural frequencies and mode shapes of the runner which are shown in Figure 7.

To characterize the heterogeneous mixture of vapor and water inside the mapped cavity volumes various simulations have been carried out considering that the average void ratio could range from 0.1 to 0.9. However, the numerical results indicate that the central core of the cavities present a higher void ratio meanwhile the outer layers of the volumes have lower void ratios that approach the zero value for pure water. Consequently an average void ratio of 0.8 with an average density of 200 kg/m3 and a sonic velocity of 29.72 m/s has been considered to be an accurate choice [START_REF] Brennen | Homogeneous bubbly flows -in Cavitation and bubble dynamics[END_REF], as shown in Figure 6, for the mapped volumes obtained with isocontours of α = 0.1. The natural frequencies in air, fair, water, fwater, and with cavitation, fG1, and the frequency reduction ratios, FRR's, calculated in relation to the natural frequencies in air are indicated in Table 2.

To calculate the possible excitations and their diametrical nodes, the classical RSI formula [START_REF] Kubota | Vibration of rotating bladed disc excited by stationary distributed forces[END_REF]Tanaka, 2011) is given:

g r n Z k m Z    
where k is the diametrical mode, Zg is the number of guide vanes, Zr is the number of runner blades and n an m are an arbitrary integers.

Taking into account the current runner design, an excitation with a diametrical mode number of k = 4 is obtained for n = 1 and m = 2 that corresponds to a frequency of 24fn = 100 Hz from the rotating point of reference. This result predicts that there is a risk of resonance if the 24fn excitation frequency coincides with the natural frequency of the 4ND runner mode of vibration.

Considering the numerical accuracy of the acoustic FSI method, a ±5% error has been added to the frequency results of 4ND. In particular, the frequency at which the mode achieves the 100 Hz value corresponds to an average void ratio of 0.8 which is a feasible estimation.

Conclusion

The acoustic FSI method has been applied considering surrounding water, heterogeneous vapor cavities and solid boundaries to analyze the added mass effects on a Francis turbine runner with the aim of understanding the reason for an actual catastrophic failure.

The numerical results show that the presence of the cavities is indeed able to affect the dynamic characteristics of the structure. A set of feasible hydrodynamic parameters have been found that permit to hypothesize the tuning of the RSI pressure fluctuations with the 4ND mode of vibration in terms of natural frequency and rotating mode that would explain a fast blade cracking due to fatigue loading.

Nevertheless, more studies with different types of turbines compared against experimental data will be necessary to further validate the methodology used to simulate added mass effects under cavitation.

Figure 1 .

 1 Figure 1. Schematic of the unit.

Figure 2 .

 2 Figure 2. Photograph of the broken runner.

  and the nucleation site volume fraction αnuc=5×10 -4 .

Figure 3 .

 3 Figure 3. Visualization of cavitation vortexes in the turbine runner. The predicted cavities can be clearly identified in Figure 3. The largest cavities are observed in condition G1 at 16 MW. In each of the runner channels, a single inter blade vortex cavity is observed that detaches from the band and extends towards the blade suction side close to the trailing edge.

Figure 4 .

 4 Figure 4. Mesh of the Francis turbine runner with the surrounding water In order to simulate the main cavities in the acoustic FSI model, the shapes of the vapor volumes calculated with CFD have been mapped into the FEM keeping unchanged the loads and physical geometries. The nodes of the elements located on the cavity surface have been transferred to the FE model using mapping nodes to the nearest nodes, as shown in Figure 5.

Figure 5 .

 5 Figure 5. Mapping of cavity shapes into FE model (right) from CDF simulation results

Figure 6 .Figure 7 .

 67 Figure 6. Speed of sound in a bubbly air/water mixture at

Table 1 .

 1 Properties of the Francis turbine runner material.

	Properties Value	Young's modulus (GPa) 201	Density (kg/m 3 ) 7850	Poisson's ratio 0.3

Table 2 .

 2 Natural frequencies and frequency reduction ratios in air and water.

	Mode	0ND	1ND	2ND	3ND	4ND
	fair	83.67	117.9	149.81	228.83	270.52
	fwater	59.87	62.99	69.03	84.85	90.51
	FRRwater	28.4%	46.6%	53.9%	62.9%	66.5%
	fG1	61.65	64.45	68.65	87.02	98.406
	FRRG1	26.3%	45.0%	53.0%	60.7%	63.4%

ACKNOWLEDGMENTS

This work was supported by National Natural Science Foundation of China (Nos. 51279083 and 51409148).