Donghyuk Kang
email: kang@me.aoyama.ac.jp

Satoshi Yamazaki

Shusaku Kagawa

Byungjin An

Motohiko Nohmi

Kazuhiko Yokota

Effect of accumulators on cavitation surge in hydraulic systems

Keywords:

The analysis of cavitation surge was performed to investigate the effect of accumulators on cavitation surge. The accumulator was modeled by using the momentum equation with a mass and a damping and a stiffness coefficients. The mass, damping and stiffness coefficients were associated with a pipe length between an accumulator and a main pipe, a valve resistance and a compliance of fluid in an accumulator, respectively. The upstream accumulator with the valve resistance had the stability effect and caused the increase of the angular velocity of cavitation surge. The downstream accumulator had the stability effect at small mass flow gain factors and caused the increase/ decrease of the angular velocity of cavitation surge at low/ large cavitation compliances. The amplitudes of flow and pressure oscillations can be reduced by the installation of the upstream accumulator.

INTRODUCTION

Since cavitation surge is caused by the interaction of cavitation with hydraulic systems [START_REF] Young | Study of Cavitating Inducer Instabilities[END_REF], an installation of an accumulator in hydraulic systems is planned for suppression of cavitation surge. To predict the effect of accumulators on cavitation surge, the stability analysis of cavitation surge was performed for the hydraulic systems. The hydraulic systems consisted of an upstream tank, an inlet pipe, an upstream accumulator, a cavitating pump, a downstream pipe, an outlet pipe, a downstream accumulator, an outlet pipe and a downstream tank. The effect of the accumulators was qualitatively discussed in the ranges of mass flow gain factor and cavitation compliance shown in the experiments [START_REF] Brennen | The Bubbly Flow Model for the Dynamic Characteristics of Cavitating Pump[END_REF][START_REF] Rubin | An Interpretation of Transfer Function Data for a Cavitating Pump[END_REF][START_REF] Yonezawa | Numerical Evaluation of Dynamic Transfer Matrix and Unsteady Cavitation Characteristics of an Inducer[END_REF][START_REF] Kang | Estimation of the Dynamic Characteristics of a Double-Suction Centrifugal Pump in Cavitation Surge[END_REF].

METHODS

The dynamic of hydraulic systems were treated in terms of lumped-parameter models [START_REF] Cervone | Evaluation of the Dynamic Transfer Matrix of Cavitating Inducer by Means of a Simplified "Lumped-Parameter[END_REF][START_REF] Kang | Analytical study of Cavitation Surge in Hydraulic system[END_REF] which simplifies the description of the physical effects between two measuring points. The lump-parameter model is usally considered valid when the dimensions of a hydralic system are shorter than the acoustic wave length at the considered frequency.

For simplicity of the analytical model, the following assumuptions are adopted in the present analysis.

(1) The flow is one-dimensional.

(2) The flow and pressure oscillations are expressed as a temporally harmornic oscillation and infinitesimal. Thus, the second or higher order quantities are negligible.

(3) The working fluid is incompressible. The pressure loss is considered as the loss coefficient. (4) The elastic deformation of all pipes is negligible and the cross-sectional area of a pipe is constant. [START_REF] Kang | Estimation of the Dynamic Characteristics of a Double-Suction Centrifugal Pump in Cavitation Surge[END_REF] The compliance of tanks is large and thus the pressure oscillation inside tanks is negligible. (6) The cross-section area of all pipes connected to a main pipe from an accumulator is the same as the horizontal section area of accumulators. Under the above assumptions, the dimensionless pressure and flow oscillations can be written in the complex form as follows.

t j t j e t e t         ) ( , ) (     (1)
Here,

/ D U t t  
is the dimensionless time and

* * * /U D    
is the dimensionless angular frequency. A * is the cross-sectional area at the suction pipe and D * is the diameter of the impeller.  * is the angular velocity and U * is the tip velocity of the impeller.  * is the density of the working fluid and j is the imaginary unit. First, we consider the flow through the pipes shown in Fig. 1. By putting Eq.( 1) into the dimensionless unsteady continuity and momentum equations and subtracting the steady terms of them from themselves, the following equations can be introduced. 1. Default analytical parameters

u i    (2) o 2    (3)
1 1 1 1 1 i 2 2 ˆ      l j l    (4) u u u u u 1 u 2 2 ˆ       l j l    (5) d d d d d d 2 2 2 ˆ       l j l    (6) u u u u u 2 2 2 ˆ      l j l   (7)
Here, the subscripts i, u, 1, d, 2, o indicate the inlet tank, the inlet pipe, the upstream pipe, the downstream pipe, the outlet pipe and the outlet tank, respectively.  is the loss coefficient.

        j l , ( 8
) ad ad ad ad ad 2 2 2 ˆ         j l . (9)
Here,  is the valve resistance. The subscripts au and ad indicate the upstream and downstream accumulators, respectively. The continuous equations for the accumulators can be expressed as

au 1 u      , ( 10
) ad d 2      . ( 11
)
With respect to the fluid inside the accumulators, the flow coefficients au  and ad  can be obtained as following equations.

au au au au au au        j c v j     ( 12
)
ad ad ad ad ad ad        j c v j     (13)
Here,

* * * / D A v v 
is the dimensionless fluid volume of the accumulator.  is the specific heat ratio and c is called the compliance. By putting Eqs.( 12) and (13) in Eqs.( 8) and ( 9), the momentum equations with the mass and damping and stiffness coefficients were introduced as follows.

0 ) ( 1 2 2 1 au au au au au 2 au            c j l (14) 0 ) ( 1 2 2 2 ad ad ad ad ad 2 ad            c j l ( 15
)
The mass, damping and stiffness coefficients corresponding to the first, second and third terms of the left sides are associated with the pipe length, the valve resistance and the compliance of fluid, respectively. These parameters will be examined in the present study. A cavity volume vc * is assumed to be formed upstream of the cavitating pump. Then, the dimensionless continuity equation can be expressed as

dt t dv t t ) ( ) ( ) ( c u d    . ( 16
)
Here,

* * * c c / D A v v  is the dimensionless cavity volume.

The change of the cavity volume c

dv can be considered to be functions of the upstream cavitation number

2 * * * v * u u / ) ( 2 U p p     (17)
and the upstream flow coefficient

AU Q / u u   . ( 18
)
Thus, it can be written as

  u u u u c u u c c u u         Kd Md d v d v dv                    (19)
where M and K are mass flow gain factor and cavitation compliance, respectively. The pressure rise supplied from the cavitating pump can be expressed as the total pressure rise coefficient. The total pressure rise coefficient is defined as

2 tu td p / ) ( 2 U p p     . ( 20
)
Here. ptu and ptd are the total pressures at the upstream of the cavitating pump and the downstream of the cavitating pump, respectively. Using the Bernoulli equation with the pressure rise of the cavitating pump, the following equation is obtained.

  d p u d p u d         l j      (21)
Here, lp=lp * /D * is the dimensionless inertial length of the cavitating pump. lp * is the mean value of the chord length. The first term of the right side indicates the unsteady total pressure rise of the cavitating pump. The second term of the right side shows the unsteady dynamic pressure rise. The last term of the right side indicates the inertia term.

The pressure rise of the cavitating pump p ˆcan be considered to be functions of the upstream cavitation number u and the discharge flow coefficient d on the assumption that cavitation occurs at the upstream of the cavitating pump and that the total pressure rise depends on the discharge flow rate. Thus, we can represent the unsteady pressure rise as

  u p d p u u p d d p p 2 d u            S R                   . ( 22
)
Here, Rp and Sp are called flow gain and pressure gain in the present study, respectively. From above formulations, we obtained the homogeneous linear equations. The homogeneous linear equations have the complex angular frequencies expressed as

I R      . ( 23
)
Here, R and I show the angular velocity and the damping rate, respectively. For I <0, the infinitesimal amplitude grow, which means cavitation surge. The zero damping rate indicates the onset boundary of cavitation surge.

Table 1 shows the tested parameters. These parameters are used in the present paper as the default parameters. That is, the increase of cavitation compliance has the stabilizing effect and the increase of the mass flow gain factor causes cavitation surge. Point A is located in the upper region of the onset boundary of cavitation surge for cau=0. This means that the state of Point A is stable. As cau is increased up to 100, Point A is located in the lower region of the onset boundary of cavitation surge and thus its state becomes cavitation surge. Assuming that cavitation compliance and mass flow gain factor are not changed by the modification of the hydraulic systems in the present study, we can say that the upstream accumulator has the instability effect. For cau=1, the cavitation surge region is more widened than that for cau=100. for various compliances cau with au=0 and lau=0. The horizontal and vertical axes represent cavitation compliance and the angular velocity, respectively. All results show that angular velocity decreases with the increase of cavitation compliance. The upstream accumulator causes the increase of the angular velocity of cavitation surge. This is believed to be due to the decrease of the pipe length by the installation of the upstream accumulator. Namely, it says that the upstream accumulator acts as the inlet tank. Figure 3(a) shows the stability map of cavitation surge with the downstream accumulator for various compliances cad with ad=0 and lad=0. For cad=100, the stable region is widened mainly at small mass flow gain factors as compared to the result for cad=0. Thus, the downstream accumulator has the stability effect at small mass flow gain factors. For cad=1, the onset boundary of cavitation surge at small mass flow gain factors is similar to that for cad=100 and the onset boundary of cavitation surge at large mass flow gain factors approaches to that for cad=0. For cad=100 and cad=1, the angular velocities increases/ decrease at small/ large cavitation compliances. This is different from the result of the upstream accumulator that the angular velocity increases at all cavitation compliances. Figure 4(a) shows the stability map of cavitation surge with the upstream accumulator for various resistances au with cau=100 and lau=0. Forau=0.5, the onset boundary of cavitation surge is shifted to the cavitation surge. This indicates that the stable region becomes large and that the upstream accumulator with the valve resistance has the stability effect.

RESULTS AND DISCUSSION

As the valve resistance is increased up to au=10.0, the onset boundary of cavitation surge at small mass flow gain factors is moved to the right side but the onset boundary of cavitation surge at large mass flow gain factors approaches to the result without the accumulator. Figure 4(b) shows the angular velocities with the upstream accumulator for various resistances au with cau=100 and lau=0. For au=0.5, the angular velocity is slightly decreased at large cavitation compliances as compared to the result for au=0. For au=10, the angular velocity is almost the same as that without the accumulator. Figure 5(a) shows the stability map of cavitation surge with the downstream accumulator for various valve resistances ad with cad=100 and lad=0. Forad=0.5, the onset boundary of cavitation surge is slightly shifted to the left side from that for ad=0. For ad=10, the onset boundary of cavitation surge at

c au =0 c au =100 c au =1 c ad =0 c ad =100 c ad =1 c ad =0 c ad =100 c ad =1
A large mass flow gain factors approaches to the result without the accumulator. All results show that the valve resistance of the downstream accumulator has the instability effect. Figure 5(b) shows the angular velocities with the downstream accumulator for various resistances ad with cad=100 and lad=0. The tendency, the increase/ decrease of the angular velocity at low/ large mass flow gain factors, is not changed by the valve resistance. Figure 6(a) shows the stability map of cavitation surge with the upstream accumulator for various pipe lengths lau with cau=100 and au=0. For lau=3, the stable region becomes large as compared to the result for lau=0. For lau=10, the stable region is more widened. Figure 6(b) shows the angular velocity of cavitation surge with the upstream accumulator for various pipe lengths lau with cau=100 and au=0. For lau=3, the angular velocity is smaller than the result for lau=0. For lau=10, the angular velocity becomes smaller.

Figure 7(a) shows stability map of cavitation surge with the downstream accumulator for various pipe lengths lad with cad=100 and ad=0. For lad=3, the cavitation surge region increases as compared to the result for lad=0. For lad=10, the cavitation surge region at low mass flow gain factors becomes larger. All results of Figs 3(a), 5(a) and 7(a) show that the onset boundary of cavitation surge at large mass flow gain factors is not affected by the downstream accumulator. Figure 7(b) shows the angular velocities of cavitation surge with the downstream accumulator for various pipe lengths lad with cad=100 and ad=0. As shown in Figs. 3(b) and 5(b), the angular velocity increase/decrease at low/large mass flow gain factors. We observed the intersections of the lines for lad=0, lad=3 and lad=10 and the line without the accumulator. The intersections for lad=0, lad=3 and lad=10 are shown by ①, ② and ③ , respectively. As lad is increased, the cavitation compliance of the intersection increases.

Figure 8 shows the amplitude ratios and phases of the flow and pressure oscillations with the upstream accumulator. The three cases were examined. The first case is without the accumulator, the second case is with the upstream accumulator for cau=100 au=0 and lau=0 and the final case is with the upstream accumulator for cau=100 au=10 and lau=0. Without the accumulator, the amplitudes of 1 and u are larger than those of d and 2 as shown in Fig. 8(a) and the phases of d and 2 are about -35 degrees as shown in Fig. 8(b). However, the amplitudes of all pressure oscillations

 au =0  au =10  au =0.5  au =0  au =10  au =0.5

Without accumulator

Without accumulator

 ad =0  ad =10  ad =0.5  ad =0  ad =10  ad =0.5

Without accumulator

Without accumulator are almost equal and nearly 5 times higher than the amplitude of u as shown in Fig. 8(c). It can be observed that the phase differences between the pressure and flow oscillations are about -90 degrees due to the inertia effect as shown in Fig. 8(d). For cau=100 au=0 and lau=0, the amplitudes of 1, d and 2 and all pressure oscillations are very lower than the amplitude of u. This indicates that the amplitudes of the flow and pressure oscillations can be reduced by the installation of the upstream accumulator.

The phases ofd and 2 are about 25 degrees. For cad=100 ad=10 and lad=0, all magnitudes and phases are similar to the results without the accumulator except for the phase of u.

Figure 9 shows the amplitude ratios and phases of the flow and pressure oscillations with the downstream accumulator.

For cad=100 ad=0 and lad=0, the amplitudes of d and 2 are slightly larger than those without the accumulator as shown Fig. 9(a). This result shows that the flow induced from the change of the cavity volume easily flows to the downstream of the cavitating pump. The phases of d and 2 are about 84 degrees as shown in Fig. 9(b). The amplitudes of d and 2 are largely decreased as compared to those the downstream accumulator as shown Fig. 9(c). All result show that the upstream flow and pressure oscillations are not sinigificantly affected by the downstream accumulator. For cad=100 ad=10 and lad=0, the ampliutudes and phases of the pressure and flow oscillations are similar to the results without the accumulator.

CONCLUSION

The effects of the accumulators on cavitation surge based on the one dimensional linear stability analysis by the lumpparameter models were investigated. The following conculusions are obtained.

(1) The upstream accumulator with the small valve resistance has the stability effect and causes the increase of the angular velocity.



  l * :pipe length [m] lp :dimensionless chord length of impeller * * p / D l  [-] lp * :chord length of impeller [m] M :mass flow gain factor [-] :dimensionless complex angular velocity=  D  /U  [-]

Figure 1 .

 1 Figure 1. Schematic of the analytical model in the present study

Fig. 1 .

 1 Fig.1. The momentum equations of the accumulator can be expressed as au au au au au 1 2 2 ˆ     

Figure 2 (

 2 Figure 2(a) shows the stability map of cavitation surge with the upstream accumulator for various compliances cau with au=0 and lau=0. The abscissa and the ordinate show mass flow gain factor and cavitation compliance, respectively. The upper and lower regions of the onset boundary of cavitation surge indicate the stable and cavitation surge regions, respectively. All results show that the lines, the onset boundaries of cavitation surge, indicate the positive slope.That is, the increase of cavitation compliance has the stabilizing effect and the increase of the mass flow gain factor causes cavitation surge. Point A is located in the upper region of the onset boundary of cavitation surge for cau=0. This means that the state of Point A is stable. As cau is increased up to 100, Point A is located in the lower region of the onset boundary of cavitation surge and thus its state becomes cavitation surge. Assuming that cavitation compliance and mass flow gain factor are not changed by the modification of the hydraulic systems in the present study, we can say that the upstream accumulator has the instability effect. For cau=1, the cavitation surge region is more widened than that for cau=100.

Figure 2 (

 2 Figure 2(b) shows the angular velocities at the onset boundary of cavitation surge with the upstream accumulator

Figure 3 (

 3 b) shows the angular velocities at the onset boundary of cavitation surge with the downstream accumulator for various compliances cad with ad=0 and lad=0.

Figure 2 .Figure 3 .

 23 Figure 2. Onset condition of cavitation surge and angular velocity for the upstream accumulator with  au =0 and l au =0

Figure 4 .

 4 Figure 4. Onset condition of cavitation surge and angular velocity for the upstream accumulator with c au =100 and l au =0

Figure 5 .

 5 Figure 5. Onset condition of cavitation surge and angular velocity for the downstream accumulator with c ad =100 and l ad =0

( 2 )Figure 6 .Figure 7 .

 267 Figure 6. Onset condition of cavitation surge and angular velocity for the upstream accumulator with c au =100 and  au =0

Figure 8 .Figure 9 .

 89 Figure 8. Amplitude ratios and phases of the flow and pressure oscillations for the upstream accumulator


Xiaocui Wu

Yiwei Wang
email: wangyw@imech.ac.cn

Chenguang Huang

Xianqian Wu

Chao Yu

S Y Mpo

O Tatin

Experiment on the cloud cavitating flow around an axisymmetric projectile near the free surface

Keywords:

The typical experiment on the cloud cavitating flow around an axisymmetric projectile near the free surface is carried out in a Split-Hopkinson pressure bar (SHPB) launching system in this paper. An axisymmetric projectile with flat head is used. The features of cavities are obtained in conditions with different submerged depth. Development of the cavity length and effect of free surface are discussed. The free surface can delay the shed process of the cavity on the upper side, which causes the length difference between the upper and down side. With the decrement of submerged depth, the cavity length difference increases.

INTRODUCTION

Cavitation, as a critical phenomenon, has been attractive in the past decades. Especially when the interaction between free surface and cloud cavitating flow is involved, the problem becomes more complex. Relevant studies in literature are very limited. An understanding of the interactions is still inadequate [START_REF] Faltinsen | Hydrodynamics of high-speed marine vehicles[END_REF].There are a few theoretical and numerical approaches established [START_REF] Faltinsen | The effect of gravity and cavitation on a hydrofoil near the free surface[END_REF][START_REF] Bal | A bem for the prediction of free surface effects on cavitating hydrofoils[END_REF][START_REF]High-speed submerged and surface piercing cavitating hydrofoils, including tandem case[END_REF][START_REF]The effect of finite depth on 2d and 3d cavitating hydrofoils[END_REF][START_REF] Amromin | Analysis of body supercavitation in shallow water[END_REF]. Other works involving influence of the free surface are about the supercavitating flow in shallow water [START_REF] Chen | Properties of natural cavitation flows around a 2-d wedge in shallow water[END_REF]. However, experimental researches about cloud cavitation near the free surface are still not seen in literature. This is mainly because the free surface can hardly be produced in cavitation water tunnels.

In the present paper, typical experiments are performed in a scaled underwater launch system based on an SHPB device, and an axisymmetric projectle is used as the testing model. The unsteady non-axisymmetrical characteristics of cavities evolution are obtained. The effect of the free surface on reentry jets and cavities shedding are studied, while influence of the free surface on cavities is also investigated.

EXPERIMENTAL SETUP AND MODEL

The typical experiment is carried out in a Split-Hopkinson pressure bar (SHPB) launching system( as shown in Figure 1). The scaled underwater launch system mainly consists of four parts: the launching system (1, 2, 3, and 5), the water tank (4 and 10), the stain-sampling system (6, 7, and 8), and the high-speed camera (9). The projectile is transiently accelerated by the launching system with slight disturbance on the water. The signal of the stress wave propagating in the transmission bar and the projectile can be obtained by using the strain sampling system. The high-speed camera is used to capture the trajectory and cavities features. The tank dimensions are C×B×H=2000mm×1000mm×800mm, and the top of the tank is directly exposed to the open air (as shown in Figure 2). More details about the system can be refer to [START_REF] Yu | Study of characteristics of cloud cavity around axisymmetric projectile by large eddy simulation[END_REF]. An axisymmetric projectile with flat head is used in this article (as shown in Figure 3). The length of the projectile is 150 mm, and the diameter is 37mm. Lengths of the cavity can represent the transient behavior of cavitating flows. The red line in Figure 3 marks for the cavity length on the lower side (L down ) and the black line is the cavity length on the upper side (L u p ). H is the distance from the free surface and the projectile. Lengths of cavities are gotten through pixel analysis on pictures. For example, the projectile with the length of 150mm shows about 488 pixels. Then, 1 pixel stands for 0.31mm, and the deviation is 0.31mm. The development of the cavity length will be discussed in the following. The cavitation number is defined as

σ = p out -p v 0.5ρ l V 2 ∞ ( 1
)
where p out is the static pressure of outlet, 1 atm. p v is the vapour pressure, 2.97 kPa. ρ l is the water density, 998.2 kg/m3. V ∞ is the inflow velocity, 18.5 m/ s, thus the cavitation number is σ = 0.58 . The value of V ∞ is also gotten through pixel analysis on pictures. With the same method, we can find how the head position of the projectile varies with time and calculate the velocity of the projectile. The Reynolds number is defined as

Re= ρ l V ∞ l µ ( 2
)
where l is the length of the projectile-0.15m, which is treated as the characteristic length. µ is the water dynamic viscosity, 0.001003kg/ms. Thus the Reynolds number is Re = 2.7×10 6 . Experiments of five working conditions are performed with various water depths from 5mm to 40mm.

RESULTS AND DISCUSSION

Cavity evolution in a shedding cycle

The quasi-periodic development of cavity shape is obtained from the experimental results. Photographs at different moments are shown in Figure 4. The value of H is gotten through measuring the distance from the free surface and the projectile by experiment. For the first shedding cycle of cavity evolution, the time is short so the difference of H can be neglected with the length L. In this case, the distance from the free surface and the projectile is considered as a constant value as H=15mm. The variation of cavity length is shown in Figure 5.

The cavity evolution in every cycle can be divided into two stages. Firstly, the bubble is in the growth stage (as shown in Figure 4 (a) and in Figure 5 (Stage 1-1)). Then the re-entry jet is generated and developed (as shown in Figure 4 (b-c) and in Figure 5 (Stage 1-2)). The cavity sheds in stage 2-1 and collapses in the stage 2-2 (as shown in Figure 5), in which the re-entry jet is also generated for the second time. In the stages of 1-2, 2-1 and 2-2 (as shown in Figure 5), there are a significant difference between the lengths of in the up and down sides of the projectile. Longer cavities appear alternately.

Effect of free surface

The difference of cavity lengths between up and down sides is obtained to analyze the relationship between the up and down sides (as shown in Figure 6). In the picture, the t represents the time of the development of the cavity. is the difference of cavity lengths between up and down sides. It can be seen that the difference of lengths stays approximately stable at growth stage (Stage1-1), and its length varies in a very small range at Stage1-2. There is a significant difference between the lengths of in the up and down sides of the projectile at Stage2-1, 2-2. The cavity length difference reach a maximum when the down side cavity sheds.

Re-entry jet is the key factor on the shedding of cavity and vortex, which has been widely studied. Re-entry jet is a liquid stream, which is opposite to the direction of the main flow in the cavity. Re-entry jet is produced by the adverse pressure gradient near the closure the cavities. It indicates that the formation of a re-entrant jet is the main reason for the unsteadiness. The cavity will be cut off by the re-entry jet. Then the shedding cavity collapses which leads to a high pressure in local domain. Another re-entry jet will be induced accordingly. Thus the re-entry jet plays a significant function in the vapor cloud shedding process. In the present condition under the effect of the free surface, the re-entry jet on the upper side is generated earlier and faster than that on the lower side. Consequently the frontier of the re-entry jet on the surface is oblique and the re-entry jet on the upper side reaches the shoulder of the projectile firstly. The shedding of the cavity on each side is also not synchronous, which induces the differences of lengths and shapes between cavities on the upper and lower side.

Cavity evolution varied with different depth

In order to investigate the relationship between the cavity length and the distance from the free surface and the projectile, the other condition is investigated with different H. The value of H is 10mm (as shown in figure 7 and figure 8). Refer to the evolutions aforementioned and in contrast with the condition of 15mm depth, we can see the cavity on the upper side stays approximately more stable, and its length varies in a very small range after the growth stage. The free surface can delay the shed process of the cavity on the upper side. With the decrement of distance from the free surface and the projectile, the cavity length on the upper side also increases. For example, the cavity length on the upper side keeps around 50mm at H=10mm.

The maximum length difference between cavities on the upper and lower side are shown in the curve in Figure 9. The difference increases when the depth decreases. And the values are similar in the conditions with the depths of 25mm and 40mm, which is because the effect of the free surface is negligible in these conditions. The differences are induced by the randomness of cavity fluctuation.

CONCLUSIONS

Experiments on the cloud cavitating flow around an axisymmetric projectile near the free surface have been presented in this paper. Cavity evolution in a shedding cycle, effect of free surface on the cavity lengths and the cavity lengths varied with different distance from the free surface and the projectile are discussed. The results show that:

The development of cavity shape is quasi-periodic. The cavity evolution in every cycle contains the growth of the cav- The difference is larger and the cavity is more stable on the upper side when the free surfaces is closer to the free surface.

The mechanism of the free surface effect is worthy of further analysis. And in the condition that the depth is 5mm, ventilation phenomenon occurs that the air above the free surface is injected into the cavity, which is also need to be investigated in the future.

Figure 1 .

 1 Figure 1. Underwater launch system

Figure 2 .

 2 Figure 2. A picture of water tank

Figure 3 .

 3 Figure 3. The test model

Figure 4 .

 4 Figure 4. Time evolution of cavity patterns(H=15mm)

Figure 5 .

 5 Figure 5. The comparison of the variation of cavity length between up and down sides (H=15mm))

Figure 6 .

 6 Figure 6. Time evolution of cavity patterns(H=15mm)

Figure 7 .Figure 8 .

 78 Figure 7. Time evolution of cavity patterns(H=10mm)

Figure 9 .

 9 Figure 9. The variation of maximum cavity length difference with depth

ACKNOWLEDGMENTS

This research was sponsored by National Natural Science Foundation of China (Grant numbers 11332011 and 11202215) and the Youth Innovation Promotion Association of CAS (Grant number 2015015).


Ashraful Alam
email: alam@matsue-ct.ac.jp

Manabu Takao

Toshiaki Setoguchi

Heuy Dong Kim

S Y Mpo

O Tatin

Passive Enhancement of Pressure Oscillation in Cavity-Induced Supersonic Mixing

Keywords:

The present work reports of the numerical investigation on the passive technique used to enhance the supersonic cavity-induced pressure oscillations. The unsteady flows over a shallow cavity at Mach 1.73 were modified at the leading edge by changing the leading edge geometry. 3 (three) leading edge geometries namely circular-contoured, slanted and circular contour on 30 o slant edge were used to intensify the excitation of the shear layer. The density-based algorithm in ANSYS Fluent 13.0 was used to solve the time-dependent Navier-Stokes equations. The SST k-ω model was used to model the turbulence within an unstructured mesh solver. The numerical code was validated comparing the numerical results with the experimental data. The pressure fluctuations with time at different locations in and after the cavity were analyzed. The modification of cavity leading edge shown influencing on the pressure fluctuations.

INTRODUCTION

Recent interest in supersonic combustion, such as scramjets, raise recommenced research in supersonic flow mixing. The propulsion system in scramjets requires rapid mixing of airfuel to achieve the efficient combustion and at the same time to minimize the size of the combustion chamber. The difficulty of mixing enhancement originates mainly from the inherently low growth rates of supersonic free shear layers [START_REF] Chinzei | Spreading of two-stream supersonic mixing layers[END_REF][START_REF] Papamoschou | Observations of supersonic frre shear layers[END_REF][START_REF] Papamoschou | The compressible shear layer: an experimental study[END_REF]. Variety of augmentation techniques that were employed efficiently in the subsonic flow mixing failed to work in the supersonic flow regimes and some mechanisms likely increase mixing at the expense of significant total pressure losses [START_REF] Gutmark | Mixing enhancement of supersonic free shear flows[END_REF]. Therefore, the researchers experienced interest in the use of active and passive control techniques to promote mixing. The use of cavity in supersonic free shear flow is the most renowned one of mixing enhancement techniques and the well established cavity-induced oscillatory features are utilized to enhance flow mixing [START_REF] Handa | Supersonic mixing enhanced by cavity-induced threedimensional oscillatory flow[END_REF][START_REF] Handa | On the jump in the frequency of acoustic oscillations in supersonic flows over rectangular cavity[END_REF][START_REF] Yu | Cavity-actuated supersonic mixing and combustion control[END_REF][START_REF] Sato | Advanced mixing control in supersonic airstream with a wall-mounted cavity[END_REF][START_REF] Seiner | Historical survey on enhanced mixing in scramjet engine[END_REF][START_REF] Nenmeni | Cavity-induced mixing enhancement in confined supersonic flows[END_REF].

Over last few decades, many researches have been conducted to investigate the external flows over cavity [START_REF] Charwat | An investigation of separated flows -part i: The pressure field[END_REF][START_REF] Hahn | Experimental investigation of separated flow over a cavity at hypersonic speed[END_REF][START_REF] Rockwell | Review -self-sustaining oscillations of flow past cavities[END_REF][START_REF] Zhang | An investigation of supersonic oscillatory cavity flows driven by thick shear layers[END_REF]. Some researches have also been conducted for the internal flows over cavity [START_REF] Matsuo | Effect of non-equilibrium condensation of moist air on supersonic cavity flow[END_REF][START_REF] Matsuo | A new passive method of controlling cavity-induced pressure oscillations[END_REF][START_REF] Setoguchi | A numerical study on passive control of cavity induced pressure fluctuations in two-dimensional supersonic flow[END_REF] in relation to the supersonic airintake that is an essential part of the air breathing engine. The high-speed driven-cavity flowfields are often characterized by the compression waves, expansion waves, shear layer and an oscillation inside the cavity. The supersonic air-intake flow tends to be unstable due to the separation induced by the shock / boundary layer or the shock / shear layer interactions [START_REF] Seddon | Intake Aerodynamics[END_REF][START_REF] Dougherty | Time accurate navier-stokes computation of self-excited two-dimensional unsteady cavity flows[END_REF]. The supersonic flows over the cavity accompanying with the shear layer separation from the upstream lip as shown in Fig. 1 and a series of vortices travel downstream in the cavity and strikes or reattaches on the rear end with the generation of compression waves. These compression waves then propagate upstream within the cavity and excite the shear layer at the cavity upstream. As a result, the process is sustained by completing the feedback loop as first described by [START_REF] Krishnamurty | Acoustic radiation from two dimensional rectangular cutouts in aerodynamic surfaces[END_REF] and later on by [START_REF] Rossiter | Wind tunnel experiments on the flow over rectangular cavities at subsonic and trasonic speeds[END_REF]. The modification of the shear layer continues by the subsequent excitations, and consequently resonance occurs when the frequency and phase of the shear layer match with those of the compression waves. The resonance, thus produced, largely amplifies the cavity pressure fluctuations. oscillations. That is, how is it possible to make a reflection of upstream propagating compression waves so that they can excite the shear layer with higher intensity? In order to answer this question, as a passive technique, 3 (three) leading edge geometries named circular-contoured, slanted and circular contour on 30 o slant edge were used to test the shear-layer modification in the present study. The numerical analyses were performed solving this problem. The density-based algorithm in ANSYS Fluent 13.0 was used to solve the time-dependent Navier-Stokes equations. The SST k-ω model was used to model the turbulence within an unstructured mesh solver. The pressure fluctuation with time was measured at different locations in the cavity and in the main stream region.

NUMERICAL MODELING & COMPUTATIONAL METHODOLOGY 1.1 Governing equations

Computational Fluid Dynamics (CFD) investigations were performed to simulate the supersonic flows over cavity with various leading edge geometries. The flows under study was treated as compressible, viscous, unsteady and turbulent. The governing equations used in the present computations are the time-dependent Navier-Stokes equations given in vector form as shown below:

∂W ∂Q ∂ ∂t V QdV + [F -G]dA = V HdV (1)
where the vectors W, Q, F, and G are defined as,

W =            ρ ρu ρv ρE            , Q =            p u v T            , F =             ρv ρvu + p i ρvv + p j ρvE + pv             (2)
and,

G =            0 τ xi τ yi τ i j v j + q            (3)
The vector H contains the source terms of the body forces and energy sources. In the above equations, ρ, v, E, and p are the density, velocity, total energy per unit mass, and the pressure of the fluid, respectively. τ is the viscous stress tensor, and q is the heat flux. The total energy E is related to the total enthalpy H by E = Hp/ρ where, H = h + |v| 2 /2.

The Jacobian matrix ∂W/∂Q in Eqn. ( 1) is replaced with the preconditioning matrix Γ for transforming the system into the conservation form, and the matrix Γ is given by,

Γ =            Θ 0 0 ρ T Θu ρ 0 ρ T u Θv 0 ρ ρ T v ΘH -δ ρu ρv ρ T H + ρC p            (4)
where,

ρ T = ∂ ρ ∂T | p (5)
and the parameter Θ is given by,

Θ = 1 U 2 r - ρ T ρC p ( 6
)
C p is the specific heat at constant pressure. The reference velocity U r is chosen locally such that the eigenvalues of the system remain well conditioned with respect to the convective and diffusive time scales [START_REF] Weiss | Preconditioning applied to variable and constant density flows[END_REF].

To close the governing equations, the k-ω SST model which is a two equation eddy-viscosity (Shear Stress Transport) turbulence model [START_REF] Mentor | Zonal two equation k -ω turbulence models for aerodynamic flows[END_REF][START_REF] Mentor | Two equation eddy-viscosity turbulence models for engineering applications[END_REF][START_REF] Mentor | Ten years of industrial experience with the sst turbulence model[END_REF] was employed in the present computation. This turbulence model is an effective blend of the robust and accurate formulation of the Wilcox's k-ω model in the near-wall region with the free-stream independence of the k-ε model in the far field. The turbulence kinetic energy k and the specific dissipation rate ω are determined by the following transport equations.

∂ ∂t ( ρk) + ∂ ∂x i ( ρku i ) = ∂ ∂x j (Γ k ∂k ∂x j ) + Gk -Y k (7) ∂ ∂t ( ρω) + ∂ ∂x i ( ρωu i ) = ∂ ∂x j (Γ k ∂ω ∂x j ) + G ω -Y ω + D ω (8)
In these equations, Gk represents the generation of the turbulence kinetic energy due to the mean velocity gradients, G ω represents the generation of ω, and D ω is the crossdiffusion term. Γ k and Γ ω represent the effective diffusivity of k and ω, respectively. Y k and Y ω represent the dissipation of k and ω due to turbulence.

The effective diffusivities for the model are given by,

Γ k = µ + µ t σ k (9)
Γ ω = µ + µ t σ ω (10)
where, σ k and σ ω are the turbulent Prandtl numbers for k and ω, respectively, and they are given by;

σ k = 1 F 1 /σ k,1 + (1 -F 1 )/σ k,2 (11)
and

σ ω = 1 F 1 /σ ω,1 + (1 -F 1 )/σ ω,2 (12)
The turbulent viscosity µ t is computed as;

µ t = ρk ω 1 max[1/α * , SF 2 /a 1 /ω] ( 13
)
where S is the strain rate magnitude, defined as;

S ≡ √ 2S i j S i j (14)
and

S i j = 1 2 ∂u j ∂ x i + ∂u j ∂ x i (15)
The coefficient α * damps the turbulent viscosity causing a low-Reynolds-number correction. It is given by,

α * = α * ∞ α * 0 + Re t /R k 1 + Re t /R k (16)
where,

Re t = ρk µω (17)
The blending functions F 1 and F 2 are given by,

F 1 = tanh(Φ 4 1 ) (18)
Φ 1 = min       max √ k 0.09ωy , 500µ ρy 2 ω , 4ρk σ ω,2 D + ω y 2       (19)
D + ω = max 2ρ 1 σ ω,2 1 ω ∂k ∂ x j ∂ω ∂ x j , 10 -10 (20)
and,

F 2 = tanh(Φ 2 2 ) (21)
Φ 2 = max       2 √ k 0.09ωy , 500µ ρy 2 ω       (22)
where, y is the distance to the next surface and D + ω is the positive portion of the cross-diffusion term.

The cross-diffusion term D ω is defined as;

D ω = 2(1 -F 1 ) ρσ ω,2 1 ω ∂k ∂ x j ∂ω ∂ x j (23)
Model constants are given as;

σ k,1 = 1.176, σ ω,1 = 2.0, σ k,2 = 1.0, σ ω,2 = 1.168 a 1 = 0.31, α ∞ = 0.52, α * = α * ∞ = 1.0, α * 0 = β i /3 β i = 0.072, R k = 6

Numerical methods

The density-based algorithm of ANSYS Fluent 13.0 [START_REF]Fluent User Guide[END_REF] was used to solve the preconditioned governing equations. A fully implicit method was implemented on the present spatial domain. The convective fluxes were formulated using the Roe's flux difference splitting scheme [START_REF] Roe | Approximate riemann solvers, parameter vectors and difference schemes[END_REF]. The 3rd order accuracy was conceived from the original MUSCL [START_REF] Van Leer | Towards the ultimate conservative difference scheme v: A second order sequel to godunov's method[END_REF] finite volume scheme that is a blend of central differencing and second-order upwind schemes, in which the physical domain is subdivided into numerical cells and the integral equations were applied to each cell. Second-order central difference scheme was used for viscous terms. For the time derivatives, an implicit multistage time stepping scheme was used.

Computational conditions

A cavity was placed in the bottom wall of a supersonic channel in order to model the supersonic flow over cavity. Schematics of the cavity in supersonic channel is shown in Fig. 1(b). The nominal channel height is of H=24 mm. The leading edge of the cavity was modified into various geometrical shapes to enhance the flowfield oscillations. The geometric detail of various leading edge shapes is shown in Figs. 2. Three geometric configurations of the leading edge, circularcontoured (Fig. 2(a)), slanted (Fig. 2(b)) and circular contour on 30 o slant edge (Fig. 2(c)), were used in the present computational work and they are referred as CLE-R12∼24, SLE0∼30 and S30-CLE-R12∼36, respectively. In case of circular-contoured and circular contour on 30 o slant leading edges, the radius of curvature R 1 and R 2 were changed between 12 mm and 24 mm with 4 mm intervals and between 12 mm and 36 mm with 6 mm intervals, respectively; while the slant angle θ was varied as 0, 15, 30 and 45 degrees with the vertical axis. The cavity depth was kept constant at D=12 mm. The measuring points of pressure fluctuations were set at locations of a, b, c, d, e, f 1 , f 2 , f 3 , and g as indicated in Fig. 1(a). The information about the length of cavity is shown on the figures. The upstream domain was extended to 5D upstream from the cavity inlet, while the domain downstream of cavity was extended to 10D. This wide computational domain is required to ensure a domain-independent solution. The computational area was meshed with the structural mesh elements using ICEM 13.0. The typical mesh used is shown in Fig. 1(b). The grids were densely clustered near the boundary or in the shear layers and in the cavity in order to provide more reasonable predictions. The fineness of the computational grids was examined through the several preliminary tests to ensure that the obtained solutions were independent of the grids employed. According to those results, the grids used in the present work are 200×80 in the channel and 60×50 in the cavity.

The supersonic incoming flow having free stream Mach number and Reynolds numbers of M=1.78 and Re=2.3×10 6 , respectively, flows over the cavity. The inlet total pressure and temperature were p 0 =101.3kPa and T 0 =293 K, respectively. The boundary conditions used were the inflow conditions and pressure outlet conditions of ANSYS Fluent 13.0 at the inlet and outlet boundaries of the computational domain, respectively. Adiabatic and no-slip wall conditions were applied to the solid wall surfaces. A solution convergence was obtained when the residuals for each of the conserved variables were reduced below the order of magnitude 6. Another convergence criterion was to check the conserved quantities directly through the computational boundaries. The net mass flux was investigated when there was an applicable imbalance through the computational boundaries.

RESULTS AND DISCUSSION

To validate the computational code developed for the present numerical simulations, the supersonic flow (inlet Mach number M inlet =1.83) over the 2D open cavities of length to depth ratios L/D of 1.0 and 2.0 was investigated and the predicted results were compared with the experimental and theoretical results [START_REF] Heller | The physical mechanism of flow-induced pressure fluctuations in cavities and concepts for their suppression[END_REF][START_REF] Zhang | An investigation of supersonic oscillatory cavity flows driven by thick shear layers[END_REF][START_REF] Nishioka | Some thoughts on the mechanism of supersonic cavity flow oscillation, part 2: A new formula for the oscillation frequency[END_REF][START_REF] Takakura | Numerical study on supersonic internal cavity flows: What causes the pressure fluctuations[END_REF]. Figure 4 shows the comparison of Strouhal numbers St. The abscissa and ordinate in the figure denote the cavity length to depth ratio, L/D and Strouhal number, St (normalized frequency: f L/U ∞ ), respectively. The solid and broken lines are obtained from the formulae of predicting oscillation frequencies proposed by Nishioka et al. [START_REF] Nishioka | Some thoughts on the mechanism of supersonic cavity flow oscillation, part 2: A new formula for the oscillation frequency[END_REF] and by Heller and Bliss [START_REF] Heller | The physical mechanism of flow-induced pressure fluctuations in cavities and concepts for their suppression[END_REF], respectively. The open circle represents the experimental results reported by Zhang and Edward [START_REF] Zhang | An investigation of supersonic oscillatory cavity flows driven by thick shear layers[END_REF] and by Takakura et al. [START_REF] Takakura | Numerical study on supersonic internal cavity flows: What causes the pressure fluctuations[END_REF]. The closed diamond represents the predicted results from the present simulation. The comparison shows that the predicted results have a good agreement with the experimental data [START_REF] Zhang | An investigation of supersonic oscillatory cavity flows driven by thick shear layers[END_REF][START_REF] Takakura | Numerical study on supersonic internal cavity flows: What causes the pressure fluctuations[END_REF] and with the theoretical results [START_REF] Heller | The physical mechanism of flow-induced pressure fluctuations in cavities and concepts for their suppression[END_REF][START_REF] Nishioka | Some thoughts on the mechanism of supersonic cavity flow oscillation, part 2: A new formula for the oscillation frequency[END_REF].

As discussed in earlier section, 3 (three) geometric configurations of the leading edge were tested in the present study. At the first stage, circular-contoured (Fig. 2(a)), slanted (Fig. 2(b)) leading edges were tested. Analyzing the results, the circular contour edge on 30 o slant was used to test in order to clarify the best suited leading edge shape which can excite the shear-layer with higher intensity.

The pressures were measured at 9 (nine) different locations of a, b, c, d, e, f 1 , f 2 , f 3 , and g as indicated in Fig. 1(a) to Analyzing Figs. 8, 9 and 10, it can be seen that the peaks of power spectra of pressures at point d and f 1 are highest for the cavity having the CLE-R24, SLE30 and S30-CLE-R24 leading edge configurations. Here, the high peak indicates the high energy of turbulent fluctuations in the flowfield. In addition, the higher the amplitude of pressure fluctuation and peak of power spectra the higher flow mixing.

CONCLUSION

The present numerical works was performed investigating the passive technique to intensify the shear-layer excitation in the supersonic flows over cavity. The intense shear-layer excitation results in the increase of cavity-induced pressure oscillations. The leading edge of the cavity was modified and

Figure 1 .

 1 Figure 1. The basic configuration of high-speed flows over cavity

Figure 2 .

 2 Figure 2. Details of computational domain and grids

( a )Figure 3 .

 a3 Figure 3. Details of the various leading edge geometry

Figure 4 .

 4 Figure 4. Comparison of the Strouhal numbers St

Figure 5 .

 5 Figure 5. Time histories of pressure for the cavity with circular-contoured leading edge (CLE-R12∼24)

Figure 6 .

 6 Figure 6. Time histories of pressure for the cavity with slanted leading edge (SLE0∼45)

Figure 7 .

 7 Figure 7. Time histories of pressure for the cavity with circular contour to 30 o slant edge (S30-CLE-R12∼36)

  3 (three) geometric configurations namely circular-contoured, slanted and circular contour on 30 o slant edge were tested in the present study. The time dependent Navier-Stokes equations with the SST k-ω turbulence model were solve numerically to analyze the problem. The pressure fluctuations with time were measured at different locations in the cavity and in the main stream region. The pressure fluctuations were found high for the leading edges of circular-contoured with the curvature radius of 24 mm (CLE-R24), 30 o slanted (SLE30) and circular

Figure 8 .Figure 9 .

 89 Figure 8. Power spectra of pressure for the cavity with circular-contoured leading edge (CLE-R12∼24)

Figure 10 .

 10 Figure 10. Power spectra of pressure for the cavity with circular contour to 30 o slant edge (S30-CLE-R12∼36)


Ching-Yeh Hsin
email: hsin@mail.ntou.edu.tw

Shih-Yun Wang

Chiung-Chih Lin

Jiahn-Horng Chen

Forng-Chen Chiu

Hydrodynamic Analysis and Design of Marine Current Turbine Blades

Keywords:

In this paper, two design procedures of the marine current turbine blade geometries are presented. The first one is similar to the propeller design method, the lifting line method is used for the design of the optimum circulation distribution, and the lifting surface method is then adopted for the blade geometry design. The second design procedure uses Genetic Algorithm method to design the turbine blade geometry. Hydrodynamic performances of the marine current turbine are then computed and analyzed by the potential flow boundary element method and the viscous flow RANS method. Two design examples, including a 20kw floating type current turbine, are demonstrated in the paper, and the design results show the geometries designed by both methods satisfy the design goal; however, the geometry designed by Genetic Algorithm method has a better result. It is believed that both methods are applicable for the current turbine blade designs.

.

INTRODUCTION

When extracting the ocean current energy, two most common devices are the horizontal axis current turbine and the vertical axis current turbine. For the horizontal axis current turbine design, most people use the methods for wind turbine blade designs. However, the current turbines operate in the water, and their physical behaviors are more like marine propellers. In this paper, two turbine blade design procedures are presented, and the design cases are demonstrated. One design procedure is similar to the propeller designs, a lifting line method ( [START_REF] Goldstein | On the vortex theory of screw propellers[END_REF]~ [START_REF] Kerwin | Optimum Circulation Distributions for Single and Multi-Component Propulsors[END_REF]) for the current turbine is first developed to obtain the optimum circulation distribution, and the blade geometry is then designed by the lifting surface vortex lattice method modified from the propeller lifting surface design method ([5]~ [START_REF] Greeley | Numerical methods for propeller design and analysis in steady flow[END_REF]). Finally the potential flow boundary element method is used to confirm the satisfaction of designed forces [START_REF] Luo | Investigation of the Current Turbine Performance by Computational Methods[END_REF]. The other design procedure is to use Genetic Algorithm [START_REF] Wang | Improvement of Marine Current Turbine Blade Design by Using the Genetic Algorithm[END_REF] to adjust the pitch angle and camber distributions of the blade, and the objective is to find a geometry which can provide the maximum torque and the minimum thrust. After completing the designs, both the boundary element method and the viscous flow RANS method ( [START_REF] Hsin | Applying Different Computational Methods to the Simulation of Flow Field around a Horizontal Axis Marine Current Turbine[END_REF]& [START_REF] Hsin | Applying the Boundary Element Method to the Analysis of Flow Around Horizontal Axis Marine Current Turbines[END_REF]) will be applied to the analysis of the performances of the current turbines for the final check of the designs.

ANALYSIS METHODS

Two computational methods are used for the current blade performance analysis, and they are the potential flow boundary element method (BEM) and viscous flow RANS method. The commercial software STAR CCM+ is used for the viscous flow computations.

The Boundary element Method (BEM) used in this paper is a perturbation potential based boundary element method, and the governing equation is:

1 ( ) 1 2 ( ) [ ( ) ] ( : ) ( ; ) 1 ( ) ( : ) B W q q S q S q p q dS n r p q n r p q q dS n r p q                (1)
In equation [START_REF] Goldstein | On the vortex theory of screw propellers[END_REF], S B denotes the body surface, and S W denotes the wake surface. is the strength of perturbation potentials, or equivalent to the dipole strength, and / n    is the source strength. r(p:q) is the distance between the panel point q and the induced point p. The term 1/ r is the potential induced by a unit strength source, and (1/ ) / r n   is the potential induced by a unit strength dipole.   is the dipole  strength in the wake from the Kutta condition, and the source strength in the wake is zero since the wake has no thickness. The discretized form of the equation ( 1) is:

, , , , , 1 1 1 1, p W P N N N M j i j j i j i m l m l P j j m l a b w i N n                (2)
In equation ( 2  over a wake panel. A wake alignment scheme is used in this boundary element method to correct the turbine blade wake geometry based on the induced velocities downstream, and a simple viscous correction is used for including the viscous effect. The turbine performances including the axial forces, torques, powers, pressure distributions and the circulation distributions can be obtained from the computations, and so are the pressure distributions and the circulation distributions.

The commercial software STAR CCM+ is used for the viscous flow computations. It solves the RANS equations by a finite volume method. The computational domain and boundary conditions are shown in Figure 1, and the grid topology is shown in Figure 2. Approximately 4 million polyhedral grids are normally used for the computations. A very dense mesh is adopted at tip of the blades in order to capture the phenomenon of the complex tip flow. Three turbulence models have been used for investigation, and they are "Realizable k-ε two layer", "Standard k-ε" and "k-ω SST". The computational results show that using three different turbulence models are equally good, and the STAR-CCM+ default turbulence model "Realizable k-ε two layer" has been chosen for subsequent computations.

We will first verify two computational methods. For both the BEM and RANS computations, the computational conditions and the geometric specification for the verification case are listed in Table 1, and the computer depiction of this marine current turbine is shown in Figure 3. For this case, the inflow is assumed to be axis-symmetrical, and the flow is assumed to be steady. Four different inflow velocities: 1.5, 2.0, 2.5, 3.0 m/s are used to investigate the performance of the marine current turbine at different tip speed ratios (TSR). To verify the computational results of different methods, we first introduce several important coefficients related to the marine current turbine performance. These coefficients are the tip speed ratio (TSR), axial force coefficient, torque coefficient and power coefficient. They are defined as follows:

2 TSR Rn V    (3)
In equation ( 3), R is the radius of turbine in meter, n is the rotational speed of turbine in rps, and V ∞ is the inflow velocity in m/s.

x Q 2 2 P Q 3 Q C , C , 1 1 2 2 P C C 1 2 w w X AV AV R F TSR AV            (4)
In equation ( 4), F X is the axial force working on turbine in Newton, Q is the torque working on turbine in Newton-meter. P W is the power of current turbine, and

2 w P nQ   .
The case we use for verification is the designed turbine blade later described in section 3.1. The computational axial force coefficients, torque coefficient and power coefficients are shown in Figure 4. The comparisons show that both methods predict similar results near the design TSR (5.19). Torque coefficient decreases as the value of tip speed ratio increases. However, the axial force coefficients for different tip speed ratios do not vary much. Figure 5 shows the comparisons of pressure distributions computed by the BEM and RANS methods at TSR=5.19, and we can see that the pressure distributions from BEM and RANS are different near the leading edge and trailing edge. These differences may be due to the viscous effects, and thus result in the discrepancies in the computational forces shown in Figure 4. In general, the discrepancies of the computational results between the BEM and RANS are acceptable, and the trends predicted by two methods are the same.

We will then validate two computational methods by the experimental data. The experiment conducted by Bahaj et. Al. [START_REF] Bahaj | Power and thrust measurements of marine current turbines under various hydrodynamic flow conditions in a cavitation tunnel and a towing tank[END_REF] is used here. In their experiment, a three-blade turbine with 0.8 meter diameter is tested at the inflow speed of 1.73m/s. Figure 6 and Figure 7 show the comparisons of computational results from both BEM and RANS methods and the experimental data for both the axial force and power coefficients. From two figures, we can see that the RANS method shows a better prediction for the axial force, and both methods have similar accuracies for the power prediction. The errors are around 5% to 10%, and the trends are predicted correctly. From the verification and validation of two computational methods, we think both methods are reliable with acceptable accuracies.

The Lagrange Multiplier Method

The first design approach is established based on the propeller design method. That is, the lifting line method is first used to obtain the optimum circulation distribution, or, the loading distribution. In the lifting line method, we can obtain the forces of the current turbine blade by Kutta-Joukowski Law:

      U F  (5)
In equation ( 5), U  is the resultant inflow, and   is the circulation distribution. Therefore, we can obtain the following equation:

            * * * * ( ) ( ) 0 0 0 ˆˆˆˆ ˆa t r a a t t t t a a a t x a t t e e e U
V r u r wr V r u r wr V r u r e V r u r e

F e F e                           e U V r u r wr V r u r       U V r u r wr V r u r   U V r u r wr V r u r           U V r u r wr V r u r   U V r u r wr V r u r U V r u r wr V r u r       U V r u r wr V r u r U V r u r wr V r u r       U V r u r wr V r u r U V r u r wr V r u r       U V r u r wr V r u r ( 6
)
Where  is the rotational speed, a V and t V are axial and tangential inflow velocities, and * a u and * t u are axial and tangential induced velocities. a F and t F are axial and tangential forces which contribute as thrust and torque. Figure 8 shows this force-velocity diagram. Since the induced velocities depend on the magnitude of the circulation, we can express the forces as functions of discrete circulation distributions. We can then optimize the circulation distribution to obtain the best power coefficient.

The Lagrange Multiplier Method is used to obtain the optimum circulation distribution. The Lagrange Multiplier Method can transfer a constrained problem to a nonconstrained problem by introducing the Lagrange multiplier λ. For the current turbine blade design, the constrained optimization problem is to find a circulation distribution which provides the best power coefficient ( PW C ), that is, the minimum axial force coefficient ( X C ) with a given torque coefficient ( Q C ). Therefore, the constraint is

* Q Q C C  , and
* Q C is the objective torque coefficient.
The design problem thus can be stated as:

* min 0 X Q Q C C C subject to      ( 7
)
As described earlier, the torque coefficient and the axial force coefficient are functions of circulation distribution:

    12 , [ , ,... ] xx T M QQ CC CC       Γ Γ Γ (8)
We thus can define the Lagrangian of this optimization problem as:

  * , ( )
X Q Q L C C C     Γ (9)
To get the minimum value of X C , we take the gradient of equation ( 9), and L  can be expressed as:

* 0 0 XQ QQ L C C CC         (10)
We can further derive to obtain:

* 0, 1,..., 0 Q X ii QQ C C for i M CC          (11)
The optimum circulation distribution thus can be obtained by solving the equation,   0 G  X , and

* , QQ L G CC           Γ X (12)
Once the optimum circulation distribution is obtained, the blade geometry, namely, the pitch and the camber distributions can be obtained by the lifting surface vortex lattice method by distributing the radial circulation distribution to the chord-wise direction. The potential flow boundary element method is then used to confirm the satisfaction of designed forces.

Genetic Algorithm Method

The second design procedure we used is based on the Genetic Algorithm method. This optimization method was inspired by Darwin's theory of evolution. The individual whose gene fit the environment better will have better potential for survival. After eliminating unsuitable individual for several generations, the remaining individuals indicate they might be the answer that we are looking for. Our objective is to design a new geometry of the current turbine that has a bigger C Q and smaller C T than the current turbine designed by the first method described in section 2.1. In this procedure, the optimization method is the Genetic Algorithm, and the computations are carried out by the BEM since it is much more efficient than the RANS method. We considered the pitch angle distributions as the gene in the individual. We designed a procedure that uses Genetic Algorithm to adjust the gene (pitch angle) to reach our objective (bigger torque and smaller axial force). The following is the entire process of the method: 1. Generate the first generation individual randomly.

    60% 40% (50%) i Q T F C C Smooth    ( 13
)
The smoothness is defined as the differences of the designed pitch and camber values at adjacent radii. 3. Calculate the fitness function of each individual.

 

Fitness i i i F P F   (14)
Fitness Function (Pi) is the probability the individual being selected to the next process. Pick individual by human selection (pick which has higher Pi) or natural selection (pick randomly). 4. Do crossover and mutation to generate the individual of the next generation. 5. Repeat step 3, 4 until the amount of the individual of the next generation is the same as the amount of the first generation. 6. Repeat step 2, 5 until there is an individual reach the objective we decided or reach the generation number we decided. The designs from this procedure will finally be double checked by the RANS method.

DESIGN RESULTS

Two design cases are shown in this paper, and the first one is to redesign the current turbine shown in Table 1, and the second case is a 20kW current turbine.

Design Case 1

For the first design case, the design point is at inflow speed 2.0m/s, and TSR=5.19 as in Table 1, and the designed CQ=0.0754. Two different designs are listed as follows:  The first one is that we use the current turbine lifting line method with the Lagrange multiplier method to obtain the "optimum circulation distribution", and use the current turbine lifting surface method to design the pitch and camber distributions. The first design is named as "original design".  Based on the same required torque, the circulation distribution is then adjusted to tip unloading from the optimum circulation distribution. The current turbine lifting surface method is then adopted to design the pitch and camber distributions from the adjusted circulation distribution. This design is named as "new design". The circulation distributions of these two designed turbines are shown in Figure 9, and one can see that the circulation distribution of the "new design" is tip unloaded, that is, less loading near the tip. Figure 10 shows the pitch and camber distributions of two designs, and it is interesting to see that the pitch distributions of two turbines are quite close; however, the camber distributions show that the "new design" is a tip unloaded current turbine since it has a relatively low camber values near the tip.

The performances of these two turbines computed by the BEM method are shown in Table 2, and the results show that both the torque and thrust are very close for two designs. In order to investigate the differences of two designs, we then look into the pressure distributions of two designs. Figures 11a and11b show the pressure distributions computed by BEM of two designs. We can see that the "new design" has a higher loading at mid-chord near the hub (r/R=0.30), and this is due to the tip unloading. The pressure distributions of two designs are very close at r/R=0.51, and this reflects the circulation distributions in Figure 9.

The "new design" shows lower loadings at mid-chord near the tip (r/R=0.81 and 0.90), and it also shows a flatter pressure distributions near the leading edge near tip. Therefore, from the pressure distributions point of view, "new design" seems to be a better design. The GA method is then used to design this turbine, and comparisons of the calculated axial force coefficients and torque coefficients using two design methods are shown in Table 3. The comparisons show that both Lagrange method and Genetic Algorithm methods have similar design results.

In Table 3, GA(2) is to decrease the geometry smooth factor, and adjust objective function from 50% to 20%. Table 3 shows the Comparisons of three design results, and it shows that GA (2) performance is better than the designs from both GA (1) and Lagrange method. We then remove GA (1) design result since its performance is not as good as GA (2). Table 4 is the comparison of the design results from Lagrange method and Genetic Algorithm method. The computational force coefficients shown in the table are from two computational methods. The comparisons show that BEM method and RANS method predict similar trends. However, the RANS method gives a larger forces and power.

We then observed the pressure distributions of the geometry designed by two different design procedures as in Figure 12. Since the Lagrange Multiplier method use the lifting line and lifting surface methods for the geometry designs, we can see that its design results have more physical considerations. For example, the pressure distributions in Figure 12 show that the pitch angles of the Lagrange Multiplier method are designed at angles of attack closer to the ideal angles of attack than those from the Genetic Algorithm method. The pressure distributions of the geometries designed by two design methods at r/R=0.5 and 0.9

Design Case 2

The second design case demonstrated here is a 20kw floating type current turbine. There are two sets of turbines for this 20kw current turbine, and each turbine should provide at least 10kw power. The geometric parameters are shown in Table 5, and the design torque is to have at least 4070 N-m (C Q =0.07388). For this design case, we are going to set the design goal to be the maximum power since we can have larger axial forces for the floating type current turbine. Therefore, the first procedure will be used for an initial design, and the GA method will then be used for the final designs. Table 6 shows the force and power coefficients of the designed geometry computed by the BEM and RANS methods. We can see that the torque coefficients are larger than the required one which means this design fulfills the design goal. The BEM is the computational method used in the GA method, and the RANS method is the method to double check the design. Since the predicted power by RANS method is larger than the BEM, it means that both computational methods predict this design reaches the design goal. Figure 13 shows the performance curves computed by RANS method, and Figure 14 shows the power curve computed by the RANS method. In Figure 14, the line with circles is the rated power of the electricity generator at different inflow speeds. Both Figures 13 and 14 are important in the operation of the current turbine.

CONCLUSIONS

We have developed two different design procedures based on the Lagrange Multiplier method and Genetic Algorithm method respectively. Two design cases, including a 20kw floating type current turbine, are demonstrated in the paper, and the design results show the geometries designed by both methods satisfy the design goal. For the analysis, two Computational methods, RANS and BEM methods, are used to predict the performances of the current turbines, and both methods predict the similar trend. It is believed that both design procedures presented in this paper are applicable, and they can be used for the current turbine blade designs in the future. To continue this research, the blade element momentum method commonly used for the wind turbine analysis will be developed to confirm two computational methods presented in this paper. Also, the current turbine performance in the unsteady inflows will be computed and investigated.

  forms of the integrations of (1/ ) / r n   and 1/ r over a panel. w represents the discrete forms of the integration of (1/ ) / r n 

Figure 1 .

 1 Figure 1. Computational domain and boundary conditions for RANS computations

Figure 2 .

 2 Figure 2. The grid topology around current turbine

Figure 3 .

 3 Figure 3. The computer depiction of a marine current turbine

Figure 4 .

 4 Figure 4. Force and power coefficients computed by two different computational methods

Figure 5 .

 5 Figure 5. The Pressure distributions computed by different computational methods at r/R=0.5 and r/R=0.7 for TSR=5.19

Figure 6 .

 6 Figure 6. The comparisons of computational results and the experimental data of the axial force coefficients

Figure 7 .

 7 Figure 7. The comparisons of computational results and the experimental data of the power coefficients

Figure 8 .

 8 Figure 8. The force-velocity diagram of a marine current turbine blade section

Figure 9 .

 9 Figure 9. The circulation distributions of two designs

Figure 10 .

 10 Figure 10. The pitch distributions (upper) and camber distributions (lower) of two designs

Figure 11b .

 11b Figure 11b. The pressure distributions at radii r/R=0.81 and r/R=0.90 of two designed turbines

Figure 12 .Figure 13 .

 1213 Figure 12. The performance curves computed by RANS method

Table 1 .

 1 The computational conditions and geometric specifications of the current turbine

	Foil Geometry	NACA66/a=0.8meanline
	Number of blade	3
	Radius(m)	0.4
	Rotational Speed(rps)	4.13
	Inflow velocities(m/s)	1.50, 2.00, 2.50, 3.00
	Tip Speed Ratio(TSR)	6.92, 5.19, 4.15, 3.46

Table 2 .

 2 The performances of two designs computed by the BEM method

	TSR=5.19	C X	C Q	C	P w
	"Original Design"	0.7464	0.07554 0.39205
	"New Design"	0.7436	0.07538 0.39122

Table 3 .

 3 The axial force coefficients C X and torque coefficients C Q computed by the BEM method for the designed geometries at the inflow velocity V=2m/s

		BEM
		Lagrange Method	GA(1)	GA (2)
	C X	0.7436	0.7403	0.6923
	C Q	0.0754	0.0757	0.0756

Table 4 .

 4 The forces and power coefficients computed by the BEM and the RANS methods and the differences between two design methods

			BEM		RANS
		Lagrange	GA (2)	Lagrange	GA (2)
	C X	0.7436	0.6923	0.7509	0.7211
	C Q	0.0754	0.0756	0.0768	0.0779
	C PW	0.3913	0.3924	0.3986	0.4043
			Differences
	C X		-6.899%		-3.969%
	C Q		+0.265%		+1.432%
	C PW		+0.281%		+1.430%

Table 5 .

 5 The geometric parameters of the 20kw floating type current turbine

	Foil Geometry	NACA66/a=0.8meanline
	Number of blade	3
	Radius(m)	5.0
	Rotational Speed(rps)	0.5
	Inflow velocities(m/s)	1.50

Table 6 .

 6 The force and power coefficients of the designed 20kw floating type current turbine

		C X	C Q *10	C	P w
	BEM	0.7194	0.08073	0.4227
	RANS	0.7700	0.08680	0.4545

ACKNOWLEDGMENTS

This research project is funded by the Ministry of Science and Technology under National Energy Program, Phase II, Taiwan, R.O.C. The authors would like to thank for their support.