Numerical Analysis of Cavitation Instabilities In Inducer Blade Cascade

The cavitation behavior of a four-blade rocket engine turbopump inducer was simulated by the CFD code Fine/TurboTM . The code was modified to take into account a cavitation model based on a homogeneous approach of cavitation, coupled with a barotropic state law for the liquid/vapor mixture [1–4]. In the present study, the numerical model of unsteady cavitation was applied to a four-blade cascade drawn from the inducer geometry. Unsteady behavior of cavitation sheets attached to the inducer blade suction side depends on the flow rate and cavitation number σ. Numerical simulations of the transient evolution of cavitation on the blade cascade were performed for nominal flow rate and different cavitation numbers, taking into account simultaneously the four blade-to-blade channels. Depending on the flow parameters, steady or unsteady behaviors spontaneously take place. In unsteady cases, sub synchronous or super synchronous regimes were observed. Some mechanisms responsible for the development of these instabilities are proposed and discussed.Copyright © 2005 by ASME


Introduction
One of the most prejudicial consequences of the cavitation in the rocket engine turbopump inducers is the generation of system and machinery instabilities. Machinery instabilities may appear for a particular range of the cavitation number, near to the performance drop. These unsteady phenomena are associated with different rotating nonsymmetrical cavitation patterns, which are characterized by sub-or supersynchronous frequencies.
Cavitation instabilities lead to downstream and upstream pressure fluctuations, which may be at the origin of system instabilities. They induce also axial and radial forces on the blades and shaft, which can disturb the bearings working and endanger the turbopump running.
For these reasons, inducer cavitation instabilities have been studied for many years in the aim of observing, better understanding, quantifying, and predicting such phenomena. Extensive studies have been carried out mainly from the viewpoint of system instability ͓1-3͔. Several experimental works have been performed to visualize global cavitation patterns and to measure the radial loads on the shaft due to cavitation instabilities ͓4-8͔. More recently, theoretical and numerical models have been developed to predict these instabilities and to analyze their origin ͓9-15͔.
In spite of those relevant works, more detailed cavitating flow analyses are needed for the understanding of local cavitation instabilities.
Through collaborations between the laboratory LEGI-Laboratoire des Ecoulements Géophysiques et Industriels, the Rocket Engine Division of Snecma, the French Space Agency CNES, and NUMECA International, a cavitation physical model was implemented in the commercial code FINE/TURBO™. The code allows steady and unsteady cavitating flows calculations in twodimensional ͑2D͒ or three-dimensional ͑3D͒ geometries, as presented in Refs. ͓16-20͔, and it can be used as a tool to analyze inducer instabilities.
In this context, and in complement to previous work ͓13,14͔, the main target of the present study is to propose a qualitative analysis of cavitation instabilities. The presented analyses are based on some numerical results obtained with the unsteady numerical code in a blade cascade geometry, derived from a 3D inducer. Subsynchronous and supersynchronous cavitation regimes have been simulated and analyzed. This paper illustrates numerical results obtained and proposes qualitative analyses of mechanisms responsible for these instability phenomena.

Models
The applied numerical code solves the conservative unsteady Reynolds-averaged Navier-Stokes ͑RANS͒ equations of a homogeneous fluid ͓16-21͔. In the present study, the Spalart-Allmaras model was used to simulate turbulent behavior ͓22͔.
Associated with these basic equations, the barotropic state law proposed by Delannoy and Kueny ͓23͔ was implemented in the code to model mass transfer in cavitating flows. This simple model involves in the hypotheses of instantaneous vaporization and condensation phenomena, as well as the no-slip condition between liquid and vapor phases. It does not take into account thermal effects, i.e., the energy equation is not considered in the calculation. The liquid-vapor mixture is characterized by a mixture specific mass given as a function of void ratio ␣, by the This barotropic law depends on the minimum speed of sound, c min , in the mixture and on the value of the density ratio, v / l , as illustrated in Figs. 1 and 2.
In the present study, the Reynolds number is in the order of 10 6 . The value of the minimum speed of sound c min is taken equal to 4.7 m / s, which is higher than the values used in our numerical previous works concerning cavitating flows in a Venturi apparatus ͓24͔. The density ratio v / l is imposed to be 0.1. These choices were made to improve the robustness of these first feasibility computations with FINE/TURBO™ code.
For the same blade cascade geometry considered in the present work, some numerical tests were carried out in a previous work to evaluate the influence of the speed of sound c min , the Reynolds number, and the used time step and mesh on the cavitation instability pattern. Most important results obtained from these influence analyses are described in Ref. ͓13͔.
Regardless of considered hypotheses, previous studies ͓24-27͔ pointed out that the barotropic approach can succeed to quantita-tive and qualitative predictions of cavitating flow global parameters ͑i.e., characteristic frequencies, vapor structure size, and pump performance͒ as well as mean local parameters ͑void ratio and velocity distribution͒.
These physical models have been implemented in the code FINE/TURBO™ developed by NUMECA International. This threedimensional structured mesh code solves the time dependent Reynolds-averaged Navier-Stokes equations. Time accurate resolutions use the dual time stepping approach proposed by Jameson ͓28͔. Pseudotime derivative terms are added to the equations. They march the solution toward convergence at each physical time step.
This kind of resolution is devoted to highly compressible flows. In the case of low-compressible or incompressible flows, its efficiency decreases dramatically. This well-known problem has been addressed by many authors and solved by introducing a preconditioner ͓21͔. This one is based on the studies presented in Refs. ͓29,30͔, and has been modified to take into account the cavitation model.
The discretization is based on a finite volume approach. We use a second order central scheme that must be associated with two artificial dissipation terms, respectively, of second and fourth orders as initially proposed by Jameson ͓28͔. The first one is activated in the strong pressure and density gradient areas. The other one is used in the whole domain. The applied form of the artificial dissipation term leads to a central second order accurate convection scheme. The pseudotime integration is made by a four-step Runge-Kutta procedure.
The physical time derivative terms are discretized with a second order backward difference scheme that ensures a second order accuracy in time.
A meaningful numerical work has been performed recently by Pouffary ͓17͔ to improve the preconditioner and the stability of numerical code for calculations of cavitating flows. A more detailed description of the code is given in Refs. ͓19-21͔. Some influence tests of numerical parameters, including the effects of second and fourth order dissipation terms, are presented in Ref. ͓19͔.
An example of studied inducer geometry is given in Fig. 3. Because 3D unsteady calculations in inducer geometries are very time consuming, a 2D approached geometry was adopted ͑even if it cannot represent completely the original three-dimensional case͒.
The studied inducer was designed to operate, even at a nominal flow rate, with nonzero incidence angle. Nevertheless, in the range of running conditions, the angle of attack remains very small for the considered inducer geometry, because of the variable thread and the associated curvature of the blades. In these conditions, the inlet backflow is not relevant and a 2D geometry approach can be used for qualitative analyses. It is important to note that tip clearance and backflow effects are neglected and, consequently, the applicability of the approach proposed in this paper might be limited to high inducer flow rates.
Hence, computations are performed in a four-blade cascade derived from an entire inducer. The transformation from the 3D geometry to the 2D blade cascade leads us to neglect the peripheral cavitation in the inducer: Only cavitation sheets attached on the blades will be considered. To obtain the 2D geometry, the 3D inducer ͑Fig. 3͒ is cut at a constant radius R c equal to 70% of the tip radius R. The resulting shape of the blade cut in the plane ͑R c , Z͒ is given in Fig. 4.
We use a 215ϫ 31 structured mesh per blade-to-blade channel, giving 26,660 internal nodes when calculating the four-blade cascade. The boundary conditions are imposed velocity at the mesh inlet, imposed static pressure at the outlet, and periodicity or connection conditions between the different channels of the blade cascade ͑Fig. 4͒. The distributions of the inflow velocity and the downstream pressure are uniform.
The time step, mesh, and turbulence model are chosen to put the attention on the low frequency fluctuations of the attached cavity, more than to the local unsteadiness in the cavitation sheet wake ͑cloud shedding͒.
v Several four-channel computations are performed at the nominal flow rate, for different cavitation numbers ͑based on outlet reference pressure͒ varying from low cavitating conditions down to the final performance drop of the cascade. The corresponding head drop chart is drawn in Fig. 5 at a nominal flow rate.
Stable configurations correspond to cavitation numbers lower than 0.65 or higher than 0.8 ͑square dots in the figure͒. For Circular dots in Fig. 5 correspond to unsteady regimes and define the range of cavitation instabilities of the cascade. Three kinds of instabilities have been observed from our calculations: a supersynchronous behavior, a subsynchronous configuration, and an unstable alternate cavitation.
Calculations are done in relative referential: Sub-or supersynchronous configurations are determined by the unsteady visualization of the cavitation sheets. For example, if the nonsymmetrical cavitation sheet pattern rotates in the same direction as the inducer, the characteristic frequency of the instability is added to the inducer rotation frequency and a supersynchronous configuration is established. In another way, subsynchronous regimes are characterized by cavitation patterns that rotate in the opposite direction as the inducer. v For downstream ϳ 0.75, two kinds of cavitation configurations have been observed. The first one is the alternate blade cavitation ͑Fig. 7͒. The system here is symmetrical, characterized by two radial opposite large vapor sheets and two opposite small ones.
For this cavitation number, the calculations pointed out also the appearance of another cavitation configuration, associated with a supersynchronous behavior. Figure 8 illustrates this kind of instability. It is characterized by a nonsymmetrical rotating pattern of vapor sheets: We observe four different sizes of cavitation sheets in the four channels.
Calculations indicate that this cavitation pattern rotates in the blade cascade in the same direction as the inducer rotation, but more quickly. v For downstream ϳ 0.7, another type of alternate blade configuration has been observed. As illustrated in Fig. 9, for this cavitation number, the cavitation pattern is initially characterized by two opposite large sheets in Channels 1 and 3, separated by two channels without cavitation. During two periods of inducer rotation, this configuration shifts and the large vapor structures change to Channels 2 and 4. This kind of configuration was studied in a previous work presented in Ref. ͓13͔.
For smaller cavitation numbers ͑ϳ0.65͒, a subsynchronous configuration appears ͑Fig. 10͒. In this case, the cavitation pattern is not symmetrical: We observe four sheets of different sizes. The largest sheet occurring in Channel 3 obstructs the flow in this channel and leads to the increase of the cavitation sheet in the upper channel. This configuration rotates in the opposite direction as the inducer rotation.
v Based on the results presented here above, we propose in this paper a qualitative analysis of the mechanisms responsible for the instability phenomena observed.
According to our calculations, and by making an analogy with rotating stall phenomenon observed in compressor inlet ͓31͔, the subsynchronous configuration seems to be associated with the obstruction of a channel by the largest sheet, as schematized in Fig. 11. This obstruction deflects the flow and modifies the angle of attack of the neighboring upper channel, which leads to a pressure decrease and to the increase of the cavitation sheet size in this channel. Otherwise, as flow rate is imposed as a constant in the inducer, the flow angle of attack in the channel below the obstructed channel should decrease, which leads to the reduction of vapor structure size in this channel. Further works are needed to improve the comparison between cavitating subsynchronous configuration and rotating stall phenomena.
To identify and analyze the mechanism responsible for the supersynchronous instabilities, we will adopt the nomenclature introduced in Fig. 11.
First, we will study the time evolution of the inlet and outlet flow rates in each channel during an inducer rotation period ͑Fig. 12͒, as well as the corresponding transient evolution of the total pressure variation ⌬P tot ͑Fig. 13͒. The difference between outlet and inlet flow rates in a channel ͑Fig. 12͒ represents the variation of volume of the vapor structure in the channel. These results concern the unsteady cavitation pattern illustrated in Fig. 8.
These results point out the following.
-The rise of the cavitation sheet size in Channel 2 ͑between Blades 1 and 2͒ leads to a partial obstruction of this channel, which induces the mass flow rate decrease. -As a consequence of the cavity length increase, we can observe the augmentation of total pressure variation in Channel 2. As a matter of fact, the cavitation sheet in Blade 1 does not reach the channel throat ͑i.e., the channel is not completely obstructed͒. The ⌬P tot ͑or blade load͒ increase is due to the rise of the angle of attack at Blade 1, which is associated with the flow rate decrease in Channel 2. -During the rotation of the small sheet from Channel 2 to Channel 1, the flow rates and ⌬P tot in Channels 3 and 4 are slightly modified.
Contrarily to the subsynchronous configuration, the propagation of the supersynchronous instability seems not to be related only to the inlet flow velocity fields. Indeed, the instability moves in the opposite direction as the velocity fields, as schematized in Fig. 14.
Hence, we think that the mechanism responsible for this kind of instability should be related mainly to the pressure field in the blade cascade. To study pressure fields, Fig. 15 illustrates the time evolution of the blade loads during the propagation of supersynchronous instability for a complete inducer rotation period ͑corresponding to Fig. 8͒.
In agreement with the flow rates and ⌬P tot time evolution analyses, we note that the loads of Blades 2 and 3 are slightly modified. Concerning Blade 2 ͑Fig. 15͑a͒͒ the pressure along the blade pressure side decreases during the development of the cavi- tation sheet in Channel 2. The load of Blade 3 ͑Fig. 15͑b͒͒ is slightly modified at the blade suction side ͑probably due to the decrease of the cavitation sheet in Channel 1͒. The loads concerning Blades 1 ͑Fig. 15͑c͒͒ and 4 ͑Fig. 15͑d͒͒ are strongly modified by the cavitation instability. The increase of the vapor sheet length at the suction side of Blade 1 is associated with a strong pressure rise, on the pressure side, near to the leading edge. This peak pressure could be the mechanism responsible for the decrease of the sheet in Channel 1.
To summarize, from the analyses of the whole of these unsteady results, we can propose the following scenario in order to explain the cavitation supersynchronous configuration ͑Figs. 8 and 16͒.
-The growth of the cavitation sheet on the Blade 1 suction side leads to the increase of the angle of attack on this blade. -A consequence of this phenomenon is the rise of the pressure at the Blade 1 pressure side, near the leading edge. -The pressure field in Channel 1, located below Blade 1, is modified by this pressure rise, which leads to the decrease of the cavitation sheet attached to Blade 4. -The presence of the cavitation sheet in Channel 4 will change the Blade 4 angle of attack, inducing the progressive rise of the attached sheet cavitation, and the continuation of the instability propagation.
The Laboratory LEGI, in collaboration with the French Space Agency CNES and the rocket engine division of the Snecma, develops, for many years, numerical codes to simulate and analyze cavitation phenomena, mainly occurring in turbopump inducer geometries. In the scope of collaboration with NUMECA International, a cavitation model was implemented in the commercial code FINE/TURBO™. It has been successfully applied to perform calculations of 3D steady cavitating flows ͓18͔, as well as unsteady 2D simulations of cavitating flows ͓16͔.
In this paper, we presented first results obtained by 2D unsteady calculation in a four-blade cascade drawn from one inducer geometry. Computations were performed at the nominal flow rate for several cavitation numbers. The resulting performance chart was presented. A good agreement was found with previous work performed for the same geometry, by using another numerical code based on a different numerical scheme ͓13͔. characterized by four identical cavitation sheets stable alternate blade cavitation ͑for ϳ 0.75͒ supersynchronous configuration ͑also for ϳ 0.75 by in-termittency͒ unstable alternate cavitation ͑for ϳ 0.7͒ subsynchronous configuration ͑for ϳ 0.65͒ Based on the qualitative analyses of numerical results ͑i.e., density, velocity, pressure fields, and time evolution of the mass flow rates and of the blade loads͒, the mechanisms responsible for cavitation instabilities were analyzed and described. An original scenario concerning supersynchronous instability origin was more particularly proposed and discussed.
The duration of the calculations presented in this paper does not allow us to perform reliable quantitative evaluation of the characteristic frequencies. In previous work ͓13͔, by using the same physical model and a different computational method, many timeconsuming calculations have been carried out to evaluate the characteristic frequencies of instabilities in this same geometry. The main contribution of the present paper, in complement of the previous one, is to propose physical scenarios to describe cavitation instabilities. This study was also a test case used to evaluate the capability of the new 3D numerical code to perform unsteady calculations and to simulate cavitation instability. This step of the work is very important to prepare 3D unsteady calculation in the entire inducer geometries.
Further work is now needed to assess the prediction capability of the model. Presently, a detailed study is being performed to improve the applied physical models by taking into account thermodynamic effects ͓32͔. Numerical model aspects are also in progress to succeed 3D unsteady cavitating calculations in inducer real geometries ͓33͔.
This research was supported by the French Space Agency CNES ͑Centre National d'Etudes Spatiales͒ and the Snecma company. The authors wish also to express their gratitude to NU-MECA International for its cooperation to the development of the numerical code. c min ϭ minimum sound speed in the medium ͑ms −1 ͒ P ϭ local static pressure ͑Pa͒ P ref ϭ reference pressure ͑Pa͒ P ! ϭ vapor pressure ͑Pa͒ R c ϭ cut radius for the passage from three dimensions to two dimensions ͑m͒ T 0 ϭ inducer rotation period ͑s͒ T ref ϭ reference time= blade passage time= T 0 / 4 ͑s͒ V ref ϭ reference velocity ͑=R c ⍀͒͑ ms −1 ͒ y+ ϭ nondimensional distance to the boundary ␣ ϭ void ratio ⍀ ϭ angular rotation speed of the inducer ͑rd s −1 ͒ ϭ mixture density ͑kg/ m 3 ͒ l ͑= ref ͒, ! ϭ liquid ͑=ref͒, vapor density ͑kg/ m 3 ͒ ϭ ͑P ref − P ! ͒ / ͑V ref 2 / 2͒ cavitation number ϭ static pressure coefficient= ͑P − P ref ͒ / ͑ L V ref 2 / 2͒ ⌬P tot ϭ total pressure variation in a channel; ⌬P tot = P downstream − P upstream ͑Pa͒